ADA1: Example cumulative project file
Erik’s NESARC Project
Erik Erhardt
October 01, 2019
- 1 Document overview
- 2 Research Questions
- 3 Data Management
- 3.1 Week03: Data Subset, Univariate Summaries And Plots
- 3.2 Data is complete
- 4 Graphing and Tabulating
- 4.1 Thursday ———
- 4.2 Categorical variables
- 4.3 Tables for categorical variables
- 4.4 Numeric variables
- 4.5 Graphing numeric variables
- 4.6 Week04: Bivariate graphs
- 4.6.1 Tuesday ———
- 4.6.2 Scatter plot (for regression): x = numerical, y = numerical
- 4.6.3 Box plots (for ANOVA): x = categorical, y = numerical
- 4.6.4 Thursday ———
- 4.6.5 Mosaic plot or bivariate bar plots (for contingency tables): x = categorical, y = categorical
- 4.6.6 Logistic scatter plot (for logistic regression): x = numerical, y = categorical (binary)
- 5 Statistical methods: inference and testing
- 5.1 Week05: Simple linear regression, logarithm transformation
- 5.1.1 Tuesday ———
- 5.1.2 Thursday ———
- 5.1.3 1. Scatter plot, add a regression line.
- 5.1.4 2. Fit the linear regression, interpret slope and \(R^2\) (R-squared) values
- 5.1.5 3. Interpret the intercept. Does it make sense?
- 5.1.6 4. Try plotting on log scale (\(x\)-only, \(y\)-only, both).
- 5.1.7 5. Does log transformation help?
- 5.2 Week06: Correlation and Categorical contingency tables
- 5.3 Week07: Inference and Parameter estimation (one-sample)
- 5.4 Week08: Hypothesis testing (one- and two-sample)
- 5.5 Week09: ANOVA and Assessing Assumptions
- 5.6 Week10: Nonparametric methods and Binomial and multinomial proportion tests
- 5.7 Week11: Two-way categorical tables and Simple linear regression, inference
- 5.7.1 Tuesday ———
- 5.7.2 Two-way categorical analysis.
- 5.7.3 Thursday ———
- 5.7.4 Simple linear regression.
- 5.7.4.1 (1 p) Plot the data and, if required, transform the variables so a roughly linear relationship is observed. All interpretations will be done on this scale of the variables.
- 5.7.4.2 (0 p) Fit the simple linear regression model.
- 5.7.4.3 (1 p) Assess the residuals for lack of fit (interpret plots of residuals vs fitted and \(x\)-value).
- 5.7.4.4 (1 p) Assess the residuals for normality (interpret QQ-plot and histogram).
- 5.7.4.5 (1 p) Assess the relative influence of points.
- 5.7.4.6 (1 p) Test whether the slope is different from zero, \(H_A: \beta_1 \ne 0\).
- 5.7.4.7 (1 p) Interpret the \(R^2\) value.
- 5.8 Week12: Logistic regression and Experiments vs Observational studies
- 5.1 Week05: Simple linear regression, logarithm transformation
- 6 Poster presentation
- 6.1 Week13: Complete poster in HW document
- 6.1.1 Tuesday ———
- 6.1.2 Thursday ———
- 6.1.2.1 Title
- 6.1.2.2 1. (1 p) Introduction
- 6.1.2.3 2. (1 p) Research Questions
- 6.1.2.4 3. (1 p) Methods
- 6.1.2.5 4. (1 p) Discussion (while this would follow your results, let’s put it here so you have a full column to show the results of the analysis of both research questions)
- 6.1.2.6 5. (1 p) Further directions or Future work or Next steps or something else that indicates there more to do and you’ve thought about it.
- 6.1.2.7 6. (1 p) References
- 6.1.2.8 7. (2 p) Results for your first research question.
- 6.1.2.9 8. (2 p) Results for your second research question.
- 6.2 Week14: Posters, finishing up
- 6.1 Week13: Complete poster in HW document
- 7 Additional items
- References (from Week02)
1 Document overview
This document is organized by Week and Class number. The in-class assignments are indicated by the Tuesday and Thursday Class numbers. Each week’s homework is often a combination of Tuesday and Thursday, with a small extension. Therefore, “fleshing out” the Tuesday and Thursday sections with a little addition is often sufficient for your homework assignment; that is, you won’t need a separate “homework” section for a week but just extend the in-class assignments. Rarely, the homework assignment is different from the in-class assignments and requires it’s own section in this document.
Consider your readers (graders):
- organize the document clearly (use this document as an example)
- label minor sections under each day (use this document as an example)
- For each thing you do, always have these three parts:
- Say what you’re going to do and why.
- Do it with code, and document your code.
- Interpret the results.
1.1 Global code options
# I set some GLOBAL R chunk options here.
# (to hide this message add "echo=FALSE" to the code chunk options)
# In particular, see the fig.height and fig.width (in inches)
# and notes about the cache option.
knitr::opts_chunk$set(comment = NA, message = FALSE, warning = FALSE, width = 100)
knitr::opts_chunk$set(fig.align = "center", fig.height = 4, fig.width = 6)
# Note: The "cache=TRUE" option will save the computations of code chunks
# so R it doesn't recompute everything every time you recompile.
# This can save _tons of time_ if you're working on a small section of code
# at the bottom of the document.
# Code chunks will be recompiled only if they are edited.
# The autodep=TRUE will also update dependent code chunks.
# A folder is created with the cache in it -- if you delete the folder, then
# all the code chunks will recompute again.
# ** If things are working as expected, or I want to freshly compute everything,
# I delete the *_cache folder.
knitr::opts_chunk$set(cache = TRUE) #, autodep=TRUE)1.2 Document
1.2.1 Naming
Note: Each class save this file with a new name, updating the last two digits to the class number. Then, you’ll have a record of your progress, as well as which files you turned in for grading.
ADA1_HW_ALL_NESARC_Project_05.RmdADA1_HW_ALL_NESARC_Project_06.RmdADA1_HW_ALL_NESARC_Project_07.Rmd…
A version that I prefer is to use a date using Year-Month-Day, YYYYMMDD:
ADA1_HW_ALL_NESARC_Project_20190903.RmdADA1_HW_ALL_NESARC_Project_20190905.RmdADA1_HW_ALL_NESARC_Project_20190910.Rmd…
1.2.2 Structure
Starting in Week03, we will concatenate all our Homework assignments together to retain the relevant information needed for subsequent classes. You will also have an opportunity to revisit previous parts to make changes or improvements, such as updating your codebook, modifying your research questions, improving tables and plots. I’ve provided an initial predicted organization of our sections and subsections using the # and ## symbols. A table of contents is automatically generated using the “toc: true” in the yaml and can headings in the table of contents are clickable to jump down to each (sub)section.
2 Research Questions
2.1 Week01: Personal Codebook
2.1.1 Class 02 Rmd, codebook
Dataset: National Epidemiologic Survey on Alcohol and Related Conditions (NESARC), with codebook wv1codebook.pdf.
Initial thinking: While nicotine dependence is a good starting point, I need to determine what it is about nicotine dependence that I am interested in. It strikes me that friends and acquaintances that I have known through the years that became hooked on cigarettes did so across very different periods of time. Some seemed to be dependent soon after their first few experiences with smoking and others after many years of generally irregular smoking behavior. I decide that I am most interested in exploring the association between level of smoking and nicotine dependence. I add to my codebook variables reflecting smoking levels (e.g., smoking quantity and frequency).
Topic of interest: Is there an association between nicotine dependence and the frequency and quantity of smoking on people up to 25 years old? The association may differ by ethnicity, age, gender, and other factors (though we won’t be able to test these additional associations until next semester in ADA2).
How I did it: I look through the codebook wv1codebook.pdf and find some variables of interest. I searched the text with “Ctrl-F” (find) to find these variables. For each variable, I copy/paste the description here, then formatted so it’s organized. You can choose to use a table or an outline format. I found this text format to be very easy to format. I retained the “frequency” of each response because it’s interesting to know, and because it was already in the codebook — this value is not required for your codebook.
2.1.2 Codebook
Dataset: NESARC
Primary association: nicotine dependence vs frequency and quantity of smoking
Key:
RenamedVarName
VarName original in dataset
Variable description
Data type (Continuous, Discrete, Nominal, Ordinal)
Frequency ItemValue Description
ID
IDNUM
UNIQUE ID NUMBER WITH NO ALPHABETICS
Nominal
43093 1-43093. Unique Identification number
Sex
SEX
SEX
Nominal
18518 1. Male
24575 2. Female
Age
AGE
AGE
Continuous
43079 18-97. Age in years
14 98. 98 years or older
SmokingStatus
CHECK321
CIGARETTE SMOKING STATUS
Nominal
9913 1. Smoked cigarettes in the past 12 months
8078 2. Smoked cigarettes prior to the last 12 months
22 9. Unknown
25080 BL. NA, never or unknown if ever smoked 100+ cigarettes
TobaccoDependence
TAB12MDX
NICOTINE DEPENDENCE IN THE LAST 12 MONTHS
Nominal
38131 0. No nicotine dependence
4962 1. Nicotine dependence
SmokingFreq
S3AQ3B1
USUAL FREQUENCY WHEN SMOKED CIGARETTES
Ordinal
14836 1. Every day
460 2. 5 to 6 Day(s) a week
687 3. 3 to 4 Day(s) a week
747 4. 1 to 2 Day(s) a week
409 5. 2 to 3 Day(s) a month
772 6. Once a month or less
102 9. Unknown
25080 BL. NA, never or unknown if ever smoked 100+ cigarettes
DailyCigsSmoked
S3AQ3C1
USUAL QUANTITY WHEN SMOKED CIGARETTES
Discrete
17751 1-98. Cigarette(s)
262 99. Unknown
25080 BL. NA, never or unknown if ever smoked 100+ cigarettes
Ethnicity
ETHRACE2A
IMPUTED RACE/ETHNICITY
Nominal
24507 1. White, Not Hispanic or Latino
8245 2. Black, Not Hispanic or Latino
701 3. American Indian/Alaska Native, Not Hispanic or Latino
1332 4. Asian/Native Hawaiian/Pacific Islander, Not Hispanic or Latino
8308 5. Hispanic or Latino
Depression
MAJORDEPLIFE
MAJOR DEPRESSION - LIFETIME (NON-HIERARCHICAL)
Nominal
35254 0. No
7839 1. Yes
IncomePersonal
S1Q10A
TOTAL PERSONAL INCOME IN LAST 12 MONTHS
Continuous
43093 0-3000000. Income in dollars
Height_ft
S1Q24FT
HEIGHT: FEET
42363 4-7. Feet
730 99. Unknown
* change 99 to NA
Height_in
S1Q24IN
HEIGHT: INCHES
Continuous
3572 0. None
38760 1-11. Inches
761 99. Unknown
* change 99 to NA
Weight_lb
S1Q24LB
WEIGHT: POUNDS
Continuous
41717 62-500. Pounds
1376 999. Unknown
* change 999 to NAAdditional variables were created from the original variables:
CREATED VARIABLES
DaySmoke
estimated number of days per month a subject smokes
Continuous (due to way estimated), assumes 30 days per month using SmokingFreq)
1-30.
TotalCigsSmoked
estimated number of cigarettes per month a subject smokes (DaysSmoke * DailyCigsSmoked)
Continuous
0-large.
CigsSmokedFac
quartiles of TotalCigsSmoked
Ordinal
ranges for each of the four quarters
SmokingFreq3
three levels of smoking frequency
Ordinal
from SmokingFreq
1. daily = 1. Every day, 2. 5 to 6 Day(s) a week, 3. 3 to 4 Day(s) a week
2. weekly = 4. 1 to 2 Day(s) a week, 5. 2 to 3 Day(s) a month
3. monthly = 6. Once a month or less
4. Never = 7. Never
Height_inches
Total height in inches
Height_ft * 12 + Height_in2.2 Week02: Literature Review
2.2.1 Tuesday ———
2.2.2 Class 03 Research questions
Dataset: National Epidemiologic Survey on Alcohol and Related Conditions (NESARC), with codebook NESARC_W1_CodeBook.pdf.
Research question:
- Is nicotine dependence [S3AQ10D] associated with smoking frequency [S3AQ3B1] and quantity [S3AQ3C1]?
- Is nicotine dependence [S3AQ10D] associated with depression [S4AQ1]?
- Is the associated between nicotine dependence [S3AQ10D] and depression [S4AQ1] different by demographics, such as Education or Sex?
2.2.3 Thursday ———
2.2.4 Class 04 Citations and Literature review
2.2.5 Citations
- Beckschäfer et al. (2014) said an interesting thing about X and Y and I will try that relationship with my variables A and B.
- Richert (2013) said something similar, but it was with adults instead of adolescents.
- There were other citations that are related and I want to read them more carefully (Murphy 2012; Dean and ebrary, Inc 2014; Beckschäfer et al. 2014)
For more documentation on bibliographies and citations with R Markdown, see http://rmarkdown.rstudio.com/authoring_bibliographies_and_citations.html. For general help with R Markdown, see https://www.rstudio.com/wp-content/uploads/2015/03/rmarkdown-reference.pdf.
Citations to appear in parentheses are cited in brackets like this (Beckschäfer et al. 2014; Richert 2013; Murphy 2012; Dean and ebrary, Inc 2014). A citation you want to refer to in text, such as Beckschäfer et al. (2014), doesn’t have the brackets.
2.2.6 Week 02 Homework Literature review
Dataset: National Epidemiologic Survey on Alcohol and Related Conditions (NESARC), with codebook NESARC_W1_CodeBook.pdf.
Research question:
- Is nicotine dependence [S3AQ10D] associated with smoking frequency [S3AQ3B1] and quantity [S3AQ3C1]?
- Google scholar search: “nicotine dependence smoking frequency”
- Citation: pdf file available for Dierker et al. (2007)
- Interesting points: Figures 2 and 3, quantity and frequency both positively related to probability of dependence.
- Others: Kandel and Chen (2000) and Caraballo, Novak, and Asman (2009)
You don’t need to include images in your literature review. I’m providing these figures and tables to illustrate what they look like.
- Is nicotine dependence [S3AQ10D] associated with depression [S4AQ1]?
- Google scholar search: “nicotine dependence depression”
- Citation: pdf file available for Breslau (1995)
- Interesting points: Table 2, Smoking and Nicotine Dependence both associated with Education. Table 3, Major depression associated with being nicotine dependent and Sex.
- Others: Breslau et al. (1998)
- Is the associated between nicotine dependence [S3AQ10D] and depression [S4AQ1] different by demographics, such as Education or Sex?
- In Question 2 we see differences by Education and Sex.
I have decided to further focus my question by examining whether the association between nicotine dependence and depression differs based on how much a person smokes. I am wondering if at low levels of smoking compared to high levels, nicotine dependence is more common among individuals with major depression than those without major depression. I add relevant depression questions/items/variables to my personal codebook as well as several demographic measures (age, gender, ethnicity, education, etc.) and any other variables I may wish to consider.
All required variables have been found and added to my personal codebook (by expanding Class 03).
Citations are at the bottom of this document in the “References” section.
3 Data Management
3.1 Week03: Data Subset, Univariate Summaries And Plots
3.1.1 Background
3.1.1.1 Purpose of study: smoking and depression in young adults
(Note: This is the last paragraph of the research proposal used to remind the reader of the research objectives.)
The present study will examine young adults from the National Epidemiologic Survey of Alcohol and Related Conditions (NESARC). The goals of the analysis will include 1) establishing the relationship between major depression and nicotine dependence; and 2) determining whether or not the relationship between major depression and nicotine dependence exists above and beyond smoking quantity.
3.1.1.2 Variables
Variables from NESARC that will be used include (already in codebook, above):
IDNUM(unique ID number with no alphabetics),MAJORDEPLIFE(has subject experienced a major depression?),CHECK321(has subject smoked cigarettes in past 12 months?),AGE(age of subject),TAB12MDX(tobacco dependence past 12 months),S3AQ3B1(usual frequency when cigarettes smoked),ETHRACE2A(ethnicity of subject),SEX(gender of subject), andS3AQ3C1(usual smoking quantity of cigarettes when smoked).
3.1.2 Tuesday ———
3.1.3 Data subset
First, the data is placed on the search path.
# data analysis packages
library(tidyverse) # Data manipulation and visualization suite
library(forcats) # Factor variables
library(lubridate) # Dates
## 1. Download the ".RData" file for your dataset into your ADA Folder.
## 2. Use the load() statement for the dataset you want to use.
##
## load("AddHealth.RData")
## load("addhealth_public4.RData")
## load("NESARC.RData")
# read data
load("NESARC.RData")
dim(NESARC)[1] 43093 3008The variables of interest are specified and a subset of the dataset including only these variables is created and stored in the data frame nesarc_sub.
# variables to include in our data subset
nesarc_sub <-
NESARC %>%
select(
IDNUM
, SEX
, AGE
, CHECK321
, TAB12MDX
, S3AQ3B1
, S3AQ3C1
, ETHRACE2A
, MAJORDEPLIFE
, S1Q10A
, S1Q24FT
, S1Q24IN
, S1Q24LB
)Check size and structure of data subset. Note that categorical variables are already Factor variables, but the levels do not have meaningful labels, yet.
[1] 43093 13'data.frame': 43093 obs. of 13 variables:
$ IDNUM : Factor w/ 43093 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ SEX : Factor w/ 2 levels "1","2": 1 2 2 1 1 2 1 2 2 2 ...
$ AGE : num 23 28 81 18 36 34 19 84 29 18 ...
$ CHECK321 : Factor w/ 3 levels "1","2","9": NA NA NA NA NA NA NA NA NA NA ...
$ TAB12MDX : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ S3AQ3B1 : Factor w/ 7 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ...
$ S3AQ3C1 : num NA NA NA NA NA NA NA NA NA NA ...
$ ETHRACE2A : Factor w/ 5 levels "1","2","3","4",..: 5 5 5 5 2 2 2 1 1 5 ...
$ MAJORDEPLIFE: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 2 1 ...
$ S1Q10A : num 17500 11000 6000 27000 42000 500 2500 24000 65000 4000 ...
$ S1Q24FT : num 5 5 5 5 6 5 5 5 5 5 ...
$ S1Q24IN : num 7 1 4 7 1 5 9 4 4 10 ...
$ S1Q24LB : num 195 127 185 180 220 140 160 120 140 150 ...3.1.3.1 EXTRA: AddHealth, joining waves 1 and 2 together
Joining AddHealth waves 1 and 2 together into a single dataset can be done if you want to use variables from when the participants were both adolecents and adults.
load("AddHealth.RData")
load("addhealth_public4.RData")
# size of datasets, to verify the data is available
dim(AddHealth)[1] 6504 2829[1] 6504 976[1] "AID"[1] "aid"# Rename W4's unique ID variable "AID" and "aid" to "ID"
# so both W1 and W4 have same variable name (lower/uppercase matters)
AddHealth <-
AddHealth %>%
dplyr::rename(ID = AID)
addhealth_public4 <-
addhealth_public4 %>%
dplyr::rename(ID = aid)
# note that the first column name has be updated
colnames(AddHealth)[1][1] "ID"[1] "ID"# Join two datasets by unique ID
AH14 <-
full_join(AddHealth, addhealth_public4)
# now we have twice as many columns
dim(AH14)[1] 6504 38043.1.4 Renaming Variables
Rename variables NewName = OriginalName.
nesarc_sub <-
nesarc_sub %>%
dplyr::rename(
ID = IDNUM
, Sex = SEX
, Age = AGE
, SmokingStatus = CHECK321
, TobaccoDependence = TAB12MDX
, SmokingFreq = S3AQ3B1
, DailyCigsSmoked = S3AQ3C1
, Ethnicity = ETHRACE2A
, Depression = MAJORDEPLIFE
, IncomePersonal = S1Q10A
, Height_ft = S1Q24FT
, Height_in = S1Q24IN
, Weight_lb = S1Q24LB
)
str(nesarc_sub)'data.frame': 43093 obs. of 13 variables:
$ ID : Factor w/ 43093 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Sex : Factor w/ 2 levels "1","2": 1 2 2 1 1 2 1 2 2 2 ...
$ Age : num 23 28 81 18 36 34 19 84 29 18 ...
$ SmokingStatus : Factor w/ 3 levels "1","2","9": NA NA NA NA NA NA NA NA NA NA ...
$ TobaccoDependence: Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
$ SmokingFreq : Factor w/ 7 levels "1","2","3","4",..: NA NA NA NA NA NA NA NA NA NA ...
$ DailyCigsSmoked : num NA NA NA NA NA NA NA NA NA NA ...
$ Ethnicity : Factor w/ 5 levels "1","2","3","4",..: 5 5 5 5 2 2 2 1 1 5 ...
$ Depression : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 2 1 ...
$ IncomePersonal : num 17500 11000 6000 27000 42000 500 2500 24000 65000 4000 ...
$ Height_ft : num 5 5 5 5 6 5 5 5 5 5 ...
$ Height_in : num 7 1 4 7 1 5 9 4 4 10 ...
$ Weight_lb : num 195 127 185 180 220 140 160 120 140 150 ...The first 6 and last 6 values are given by head() and tail().
ID Sex Age SmokingStatus TobaccoDependence SmokingFreq DailyCigsSmoked
1 1 1 23 <NA> 0 <NA> NA
2 2 2 28 <NA> 0 <NA> NA
3 3 2 81 <NA> 0 <NA> NA
4 4 1 18 <NA> 0 <NA> NA
5 5 1 36 <NA> 0 <NA> NA
6 6 2 34 <NA> 0 <NA> NA
Ethnicity Depression IncomePersonal Height_ft Height_in Weight_lb
1 5 0 17500 5 7 195
2 5 0 11000 5 1 127
3 5 0 6000 5 4 185
4 5 0 27000 5 7 180
5 2 0 42000 6 1 220
6 2 0 500 5 5 140 ID Sex Age SmokingStatus TobaccoDependence SmokingFreq
43088 43088 2 18 1 0 1
43089 43089 2 18 <NA> 0 <NA>
43090 43090 1 19 <NA> 0 <NA>
43091 43091 1 18 1 1 6
43092 43092 1 29 1 1 1
43093 43093 1 18 <NA> 0 <NA>
DailyCigsSmoked Ethnicity Depression IncomePersonal Height_ft
43088 5 1 0 4000 5
43089 NA 1 0 4000 5
43090 NA 1 0 2500 6
43091 3 1 0 3000 5
43092 20 1 0 6500 5
43093 NA 1 1 6000 5
Height_in Weight_lb
43088 3 120
43089 1 160
43090 1 175
43091 8 150
43092 10 140
43093 11 145A quick summary of the variables indicates many missing values.
ID Sex Age SmokingStatus TobaccoDependence
1 : 1 1:18518 Min. :18.0 1 : 9913 0:38131
2 : 1 2:24575 1st Qu.:32.0 2 : 8078 1: 4962
3 : 1 Median :44.0 9 : 22
4 : 1 Mean :46.4 NA's:25080
5 : 1 3rd Qu.:59.0
6 : 1 Max. :98.0
(Other):43087
SmokingFreq DailyCigsSmoked Ethnicity Depression IncomePersonal
1 :14836 Min. : 1.00 1:24507 0:35254 Min. : 0
6 : 772 1st Qu.: 6.00 2: 8245 1: 7839 1st Qu.: 8800
4 : 747 Median :15.00 3: 701 Median : 20000
3 : 687 Mean :16.81 4: 1332 Mean : 28185
2 : 460 3rd Qu.:20.00 5: 8308 3rd Qu.: 36000
(Other): 511 Max. :99.00 Max. :3000000
NA's :25080 NA's :25080
Height_ft Height_in Weight_lb
Min. : 4.0 Min. : 0.000 Min. : 62.0
1st Qu.: 5.0 1st Qu.: 3.000 1st Qu.:140.0
Median : 5.0 Median : 5.000 Median :168.0
Mean : 6.7 Mean : 6.934 Mean :197.1
3rd Qu.: 5.0 3rd Qu.: 8.000 3rd Qu.:197.0
Max. :99.0 Max. :99.000 Max. :999.0
3.1.5 Coding missing values
There are two steps. The first step is to recode any existing NAs to actual values, if necessary. The method for doing this differs for numeric and categorical variables. The second step is to recode any coded missing values, such as 9s or 99s, as actual NA.
3.1.5.1 Coding NAs as meaningful “missing”
First step: the existing blank values with NA mean “never”, and “never” has a meaning different from “missing”. For each variable we need to decide what “never” means and code it appropriately.
3.1.5.1.1 NAs recoded as numeric
The easiest is DailyCigsSmoked where NA means 0 cigarettes. We replace NA with 0. The function is.na() identifies each observation that “is an NA”. I’ve updated the codebook to indicate that the original NA values were changed.
S3AQ3C1 (DailyCigsSmoked)
USUAL QUANTITY WHEN SMOKED CIGARETTES
Discrete
17751 1-98. Cigarette(s)
262 99. Unknown
25080 BL. NA, never or unknown if ever smoked 100+ cigarettes
* change to 0 Cigarett(s)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
934 884 923 573 1070 463 269 299 49 3077 23 230 34 25 851
16 17 18 19 20 21 22 23 24 25 27 28 29 30 33
40 22 59 5 5366 1 10 2 7 155 2 3 3 909 1
34 35 37 39 40 45 50 55 57 60 66 70 75 80 98
1 30 2 1 993 8 106 2 1 241 1 12 2 47 15
99
262 nesarc_sub <-
nesarc_sub %>%
replace_na(
# replace NA with new numeric value
list(
DailyCigsSmoked = 0
)
)
table(nesarc_sub$DailyCigsSmoked)
0 1 2 3 4 5 6 7 8 9 10 11
25080 934 884 923 573 1070 463 269 299 49 3077 23
12 13 14 15 16 17 18 19 20 21 22 23
230 34 25 851 40 22 59 5 5366 1 10 2
24 25 27 28 29 30 33 34 35 37 39 40
7 155 2 3 3 909 1 1 30 2 1 993
45 50 55 57 60 66 70 75 80 98 99
8 106 2 1 241 1 12 2 47 15 262 3.1.5.1.2 NAs recoded as categorical
The next two variables require defining a new category level. Both SmokingStatus and SmokingFreq with a value of NA means “never”, so we need to create a new “Never” category and put in order with the rest of the levels. I’ve updated the codebook to indicate that the original NA values were changed.
CHECK321 (SmokingStatus)
CIGARETTE SMOKING STATUS
Nominal
9913 1. Smoked cigarettes in the past 12 months
8078 2. Smoked cigarettes prior to the last 12 months
22 9. Unknown
25080 BL. NA, never or unknown if ever smoked 100+ cigarettes
* change to 3. Never
S3AQ3B1 (SmokingFreq)
USUAL FREQUENCY WHEN SMOKED CIGARETTES
Ordinal
14836 1. Every day
460 2. 5 to 6 Day(s) a week
687 3. 3 to 4 Day(s) a week
747 4. 1 to 2 Day(s) a week
409 5. 2 to 3 Day(s) a month
772 6. Once a month or less
102 9. Unknown
25080 BL. NA, never or unknown if ever smoked 100+ cigarettes
* change to 7. Never
1 2 9
9913 8078 22
1 2 3 4 5 6 9
14836 460 687 747 409 772 102 nesarc_sub <-
nesarc_sub %>%
mutate(
# first change the "factor" variables to "character" to add a new level
SmokingStatus = as.character(SmokingStatus)
, SmokingFreq = as.character(SmokingFreq )
) %>%
replace_na(
# replace NA with new category values
list(
SmokingStatus = "3"
, SmokingFreq = "7"
)
)
table(nesarc_sub$SmokingStatus )
1 2 3 9
9913 8078 25080 22
1 2 3 4 5 6 7 9
14836 460 687 747 409 772 25080 102 3.1.5.2 Coding 9s and 99s as NAs
Second step: replace the coded missing values with NA, then remove the remaining records with missing values. Not all of these variables need this operation, since the previous subset commands removed missing values by chance. However, I’m going to run each as an extra precaution. Afterall, later, if we were to make different previous data decisions, it’s possible we could end up with missing values at this point that are incorrectly coded without this step!
1 2 3 9
9913 8078 25080 22
1 2 3 4 5 6 7 9
14836 460 687 747 409 772 25080 102
0 1 2 3 4 5 6 7 8 9 10 11
25080 934 884 923 573 1070 463 269 299 49 3077 23
12 13 14 15 16 17 18 19 20 21 22 23
230 34 25 851 40 22 59 5 5366 1 10 2
24 25 27 28 29 30 33 34 35 37 39 40
7 155 2 3 3 909 1 1 30 2 1 993
45 50 55 57 60 66 70 75 80 98 99
8 106 2 1 241 1 12 2 47 15 262 nesarc_sub <-
nesarc_sub %>%
mutate(
# replace missing values
SmokingStatus = replace(SmokingStatus , SmokingStatus %in% c("9"), NA)
, SmokingFreq = replace(SmokingFreq , SmokingFreq %in% c("9"), NA)
, DailyCigsSmoked = replace(DailyCigsSmoked, DailyCigsSmoked %in% c( 99), NA)
, Height_ft = replace(Height_ft , Height_ft %in% c( 99), NA)
, Height_in = replace(Height_in , Height_in %in% c( 99), NA)
, Weight_lb = replace(Weight_lb , Weight_lb %in% c(999), NA)
# drop unused factor levels
, SmokingStatus = fct_drop(SmokingStatus)
, SmokingFreq = fct_drop(SmokingFreq )
)
table(nesarc_sub$SmokingStatus )
1 2 3
9913 8078 25080
1 2 3 4 5 6 7
14836 460 687 747 409 772 25080
0 1 2 3 4 5 6 7 8 9 10 11
25080 934 884 923 573 1070 463 269 299 49 3077 23
12 13 14 15 16 17 18 19 20 21 22 23
230 34 25 851 40 22 59 5 5366 1 10 2
24 25 27 28 29 30 33 34 35 37 39 40
7 155 2 3 3 909 1 1 30 2 1 993
45 50 55 57 60 66 70 75 80 98
8 106 2 1 241 1 12 2 47 15 ID Sex Age SmokingStatus TobaccoDependence
1 : 1 1:18518 Min. :18.0 1 : 9913 0:38131
2 : 1 2:24575 1st Qu.:32.0 2 : 8078 1: 4962
3 : 1 Median :44.0 3 :25080
4 : 1 Mean :46.4 NA's: 22
5 : 1 3rd Qu.:59.0
6 : 1 Max. :98.0
(Other):43087
SmokingFreq DailyCigsSmoked Ethnicity Depression IncomePersonal
7 :25080 Min. : 0.000 1:24507 0:35254 Min. : 0
1 :14836 1st Qu.: 0.000 2: 8245 1: 7839 1st Qu.: 8800
6 : 772 Median : 0.000 3: 701 Median : 20000
4 : 747 Mean : 6.462 4: 1332 Mean : 28185
3 : 687 3rd Qu.:10.000 5: 8308 3rd Qu.: 36000
(Other): 869 Max. :98.000 Max. :3000000
NA's : 102 NA's :262
Height_ft Height_in Weight_lb
Min. :4.00 Min. : 0.000 Min. : 62.0
1st Qu.:5.00 1st Qu.: 3.000 1st Qu.:140.0
Median :5.00 Median : 5.000 Median :165.0
Mean :5.11 Mean : 5.279 Mean :170.6
3rd Qu.:5.00 3rd Qu.: 8.000 3rd Qu.:192.0
Max. :7.00 Max. :11.000 Max. :500.0
NA's :730 NA's :761 NA's :1376 Finally, consider using na.omit() to remove any records with missing values if it won’t cause issues with analysis. The command na.omit() removes a row if any columns have an NA in that row. This can have drastic consequences. For example, imagine that you have 7 variables and 6 variables have complete data, but the 7th variable is about three-quarters NAs! Then na.omit() will drop three-quarters of your data, even though most of your variables are complete!
There’s only one advantage I can think of for running na.omit(), and that’s if you have a large dataset and only a few observations have a few NAs. You can’t analyze data with NAs, since any calculation involving NA results in an NA. Therefore imagine two scenarios: (1) An analysis with Variables 1 and 2, each with a few NA distributed throughout; first the observations with NA for either Variable 1 or 2 will be removed, then the analysis will proceed. (2) A similar analysis, but now with Variables 1 and 3, but Variable 3 has NA for different observations than Variable 2 — therefore, different observations will be removed before analysis. In summary, these two analyses will include different observations from the dataset. Therefore, by running na.omit() you will start with a dataset with no NAs, thus every analysis will be of exactly the same set of data (at the cost of not using all the possible data for each individual analysis).
Below, I do not run na.omit() since I’m concerned about having lots of NAs for some variables.
[1] 43093 13[1] 43093 13 ID Sex Age SmokingStatus TobaccoDependence
1 : 1 1:18518 Min. :18.0 1 : 9913 0:38131
2 : 1 2:24575 1st Qu.:32.0 2 : 8078 1: 4962
3 : 1 Median :44.0 3 :25080
4 : 1 Mean :46.4 NA's: 22
5 : 1 3rd Qu.:59.0
6 : 1 Max. :98.0
(Other):43087
SmokingFreq DailyCigsSmoked Ethnicity Depression IncomePersonal
7 :25080 Min. : 0.000 1:24507 0:35254 Min. : 0
1 :14836 1st Qu.: 0.000 2: 8245 1: 7839 1st Qu.: 8800
6 : 772 Median : 0.000 3: 701 Median : 20000
4 : 747 Mean : 6.462 4: 1332 Mean : 28185
3 : 687 3rd Qu.:10.000 5: 8308 3rd Qu.: 36000
(Other): 869 Max. :98.000 Max. :3000000
NA's : 102 NA's :262
Height_ft Height_in Weight_lb
Min. :4.00 Min. : 0.000 Min. : 62.0
1st Qu.:5.00 1st Qu.: 3.000 1st Qu.:140.0
Median :5.00 Median : 5.000 Median :165.0
Mean :5.11 Mean : 5.279 Mean :170.6
3rd Qu.:5.00 3rd Qu.: 8.000 3rd Qu.:192.0
Max. :7.00 Max. :11.000 Max. :500.0
NA's :730 NA's :761 NA's :1376 Now we have a dataset with no missing values.
3.1.6 Creating new variables
Additional variables that we created from the original variables are listed in this ammendment to the original codebook.
CREATED VARIABLES
DaySmoke
estimated number of days per month a subject smokes
Continuous (due to way estimated), assumes 30 days per month using SmokingFreq)
1-30.
TotalCigsSmoked
estimated number of cigarettes per month a subject smokes (DaysSmoke * DailyCigsSmoked)
Continuous
0-large.
CigsSmokedFac
quartiles of TotalCigsSmoked
Ordinal
ranges for each of the four quarters
SmokingFreq3
three levels of smoking frequency
Ordinal
from SmokingFreq
1. Daily = 1. Every day, 2. 5 to 6 Day(s) a week, 3. 3 to 4 Day(s) a week
2. Weekly = 4. 1 to 2 Day(s) a week, 5. 2 to 3 Day(s) a month
3. Monthly = 6. Once a month or less
4. Never = 7. Never
Height_inches
Total height in inches
Height_ft * 12 + Height_in3.1.6.1 From categories to numeric
DaysSmoke estimates the days per month a subject smokes by converting SmokingFreq (a factor with 6 levels) to a numeric variable using as.numeric() and multiplying by the midpoint of the range of SmokingFreq.
nesarc_sub <-
nesarc_sub %>%
mutate(
DaysSmoke =
case_when(
SmokingFreq == 1 ~ 30 # 1. Every day
, SmokingFreq == 2 ~ 4 * 5.5 # 2. 5 to 6 Day(s) a week
, SmokingFreq == 3 ~ 4 * 3.5 # 3. 3 to 4 Day(s) a week
, SmokingFreq == 4 ~ 4 * 1.5 # 4. 1 to 2 Day(s) a week
, SmokingFreq == 5 ~ 2.5 # 5. 2 to 3 Day(s) a month
, SmokingFreq == 6 ~ 1 # 6. Once a month or less
, SmokingFreq == 7 ~ 0 # 7. Never
)
)
summary(nesarc_sub$DaysSmoke) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.00 0.00 0.00 10.96 30.00 30.00 102 3.1.6.1.1 AddHealth4 Income example
Use the AddHealth Wave 4 income variable to convert from categorical to numerical.
Note that there’s also variable H4EC2, which is already a numeric income variable.
# create "Income4" from the Section 12: Economics section.
AH14 <-
AH14 %>%
mutate(
Income4 = h4ec1
)
table(AH14$Income4)
1 2 3 4 5 6 7 8 9 10 11 12 96 98
137 120 173 167 237 269 505 578 1154 698 476 247 51 302 # create a numeric version of this categorical variable
AH14 <-
AH14 %>%
mutate(
Income4_num =
case_when(
Income4 == 1 ~ 2500 # less than $5,000
, Income4 == 2 ~ 7500 # $5,000 to $9,999
, Income4 == 3 ~ 12500 # $10,000 to $14,999
, Income4 == 4 ~ 17500 # $15,000 to $19,999
, Income4 == 5 ~ 22500 # $20,000 to $24,999
, Income4 == 6 ~ 27500 # $25,000 to $29,999
, Income4 == 7 ~ 35000 # $30,000 to $39,999
, Income4 == 8 ~ 45000 # $40,000 to $49,999
, Income4 == 9 ~ 62500 # $50,000 to $74,999
, Income4 == 10 ~ 87500 # $75,000 to $99,999
, Income4 == 11 ~ 125000 # $100,000 to $149,999
, Income4 == 12 ~ 175000 # $150,000 or more
## These two would have been coded as NA previously
## You can't code them as NA in this case_when() function.
#, Income4 == 96 ~ # refused
#, Income4 == 98 ~ # don't know
)
)
table(AH14$Income4_num)
2500 7500 12500 17500 22500 27500 35000 45000 62500 87500
137 120 173 167 237 269 505 578 1154 698
125000 175000
476 247 3.1.6.2 From numeric to numeric
Height_inches is the height of a person in inches. Since the data was collected in two variables (Height_ft and Height_in) we need to add them together for our new variable.
nesarc_sub <-
nesarc_sub %>%
mutate(
Height_inches = Height_ft * 12 + Height_in
)
summary(nesarc_sub$Height_inches) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
48.0 64.0 66.0 66.6 70.0 84.0 761 The variable TotalCigsSmoked estimates the monthly number of cigarettes a subject smokes per month by multiplying DaysSmoke times DailyCigsSmoked (the usual quantity smoked per day).
Log scale. I like log2 since its interpretation is that every unit is a doubling. Notice that I added 1 to the TotalCigsSmoked before computing the log. This is because log(0) is undefined, thus I add 1 so that log(0 + 1) = 0 and all values are defined.
nesarc_sub <-
nesarc_sub %>%
mutate(
TotalCigsSmokedLog2 = log(TotalCigsSmoked + 1, base = 2)
)
summary(nesarc_sub$TotalCigsSmokedLog2) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 0.000 0.000 3.284 8.234 11.522 297 Square root scale. Square root is another transformation that spreads out values between 0 and 1 and compresses values from 1 to infinity. It can be preferred to a log transformation in some cases.
nesarc_sub <-
nesarc_sub %>%
mutate(
TotalCigsSmokedSqrt = sqrt(TotalCigsSmoked)
)
summary(nesarc_sub$TotalCigsSmokedSqrt) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.000 0.000 0.000 7.887 17.321 54.222 297 3.1.6.3 From numeric to categories based on quantiles
The numeric variable TotalCigsSmoked is converted into a factor with four roughly equivalent numbers (quartiles) stored in each level of the factor CigsSmokedFac using the cut function. The label "Q1" indicates that the original value was in the first quarter (25%) of the ordered data, that is, they are among the smallest 25% observations; while "Q4" indicates the original value was in the last quarter (25%) of the ordered data, that is, they are among the largest 25% observations.
0% 25% 50% 75% 100%
0 0 0 300 2940 nesarc_sub <-
nesarc_sub %>%
mutate(
CigsSmokedFac =
.bincode(TotalCigsSmoked
, breaks = quantile(TotalCigsSmoked, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
, include.lowest = TRUE
)
, CigsSmokedFac = factor( CigsSmokedFac
, levels = c(1, 2, 3, 4)
, labels = c("Q1", "Q2", "Q3", "Q4")
)
)In my data, more than half of the values are 0s which is the minimum value, thus they all get grouped in “Q1”; this leaves no values for “Q2” – oh well.
Q1 Q2 Q3 Q4 NA's
25080 0 8630 9086 297 3.1.6.4 From many categories to a few
The variable SmokingFreq3 collapses the SmokingFreq from 6 down to 3 categories. This is an example to show you how to do this with your own variables. Some of you will have categorical variables that have a dozen or so categories; having so many categories makes interpretation very difficulty — fewer is easier.
## SmokingFreq3 variable
# notice that I'm doing this before I've assigned the labels
# if you do this after assigning labels (with the factor() command), then you'll want to use "as.numeric()" as below
# inspect original variable
table(nesarc_sub$SmokingFreq)
1 2 3 4 5 6 7
14836 460 687 747 409 772 25080
1 2 3 4 5 6 7
14836 460 687 747 409 772 25080 # Note: it is sometimes easier to use the numeric "levels" of a factor variable, than the character "labels"
# initialize new variable with NAs
nesarc_sub <-
nesarc_sub %>%
mutate(
SmokingFreq3 =
fct_collapse(
SmokingFreq
, "1" = c("1", "2", "3") # 1. daily
, "2" = c("4", "5") # 2. weekly
, "3" = "6" # 3. monthly
, "4" = "7" # 7. never
)
)
# new variable
table(nesarc_sub$SmokingFreq3)
1 2 3 4
15983 1156 772 25080 3.1.6.5 Review results of new variables
Look at the first 10 observations in the dataset, with special attention to the new variables.
ID Sex Age SmokingStatus TobaccoDependence SmokingFreq DailyCigsSmoked
1 1 1 23 3 0 7 0
2 2 2 28 3 0 7 0
3 3 2 81 3 0 7 0
4 4 1 18 3 0 7 0
5 5 1 36 3 0 7 0
6 6 2 34 3 0 7 0
7 7 1 19 3 0 7 0
8 8 2 84 3 0 7 0
9 9 2 29 3 0 7 0
10 10 2 18 3 0 7 0
Ethnicity Depression IncomePersonal Height_ft Height_in Weight_lb
1 5 0 17500 5 7 195
2 5 0 11000 5 1 127
3 5 0 6000 5 4 185
4 5 0 27000 5 7 180
5 2 0 42000 6 1 220
6 2 0 500 5 5 140
7 2 1 2500 5 9 160
8 1 0 24000 5 4 120
9 1 1 65000 5 4 140
10 5 0 4000 5 10 150
DaysSmoke Height_inches TotalCigsSmoked TotalCigsSmokedLog2
1 0 67 0 0
2 0 61 0 0
3 0 64 0 0
4 0 67 0 0
5 0 73 0 0
6 0 65 0 0
7 0 69 0 0
8 0 64 0 0
9 0 64 0 0
10 0 70 0 0
TotalCigsSmokedSqrt CigsSmokedFac SmokingFreq3
1 0 Q1 4
2 0 Q1 4
3 0 Q1 4
4 0 Q1 4
5 0 Q1 4
6 0 Q1 4
7 0 Q1 4
8 0 Q1 4
9 0 Q1 4
10 0 Q1 43.1.7 Labeling Categorical variable levels
Look at summary of dataset to see the state of the labelling before we make changes.
ID Sex Age SmokingStatus TobaccoDependence
1 : 1 1:18518 Min. :18.0 1 : 9913 0:38131
2 : 1 2:24575 1st Qu.:32.0 2 : 8078 1: 4962
3 : 1 Median :44.0 3 :25080
4 : 1 Mean :46.4 NA's: 22
5 : 1 3rd Qu.:59.0
6 : 1 Max. :98.0
(Other):43087
SmokingFreq DailyCigsSmoked Ethnicity Depression IncomePersonal
7 :25080 Min. : 0.000 1:24507 0:35254 Min. : 0
1 :14836 1st Qu.: 0.000 2: 8245 1: 7839 1st Qu.: 8800
6 : 772 Median : 0.000 3: 701 Median : 20000
4 : 747 Mean : 6.462 4: 1332 Mean : 28185
3 : 687 3rd Qu.:10.000 5: 8308 3rd Qu.: 36000
(Other): 869 Max. :98.000 Max. :3000000
NA's : 102 NA's :262
Height_ft Height_in Weight_lb DaysSmoke
Min. :4.00 Min. : 0.000 Min. : 62.0 Min. : 0.00
1st Qu.:5.00 1st Qu.: 3.000 1st Qu.:140.0 1st Qu.: 0.00
Median :5.00 Median : 5.000 Median :165.0 Median : 0.00
Mean :5.11 Mean : 5.279 Mean :170.6 Mean :10.96
3rd Qu.:5.00 3rd Qu.: 8.000 3rd Qu.:192.0 3rd Qu.:30.00
Max. :7.00 Max. :11.000 Max. :500.0 Max. :30.00
NA's :730 NA's :761 NA's :1376 NA's :102
Height_inches TotalCigsSmoked TotalCigsSmokedLog2 TotalCigsSmokedSqrt
Min. :48.0 Min. : 0.0 Min. : 0.000 Min. : 0.000
1st Qu.:64.0 1st Qu.: 0.0 1st Qu.: 0.000 1st Qu.: 0.000
Median :66.0 Median : 0.0 Median : 0.000 Median : 0.000
Mean :66.6 Mean : 187.5 Mean : 3.284 Mean : 7.887
3rd Qu.:70.0 3rd Qu.: 300.0 3rd Qu.: 8.234 3rd Qu.:17.321
Max. :84.0 Max. :2940.0 Max. :11.522 Max. :54.222
NA's :761 NA's :297 NA's :297 NA's :297
CigsSmokedFac SmokingFreq3
Q1 :25080 1 :15983
Q2 : 0 2 : 1156
Q3 : 8630 3 : 772
Q4 : 9086 4 :25080
NA's: 297 NA's: 102
Informative labels are given to the factor levels. The order of the levels is also rearranged for the variables SmokingFreq, TobaccoDependence, and Sex.
# Label factors with meaningful names in the same order as the codebook
nesarc_sub$Depression <-
factor(
nesarc_sub$Depression
, labels = c("No Depression"
, "Yes Depression"
)
)
nesarc_sub$SmokingStatus <-
factor(
nesarc_sub$SmokingStatus
, labels = c("Smoked past 12 months"
, "Smoked prior to 12 months"
, "Never smoked"
)
)
# Specify the order of the levels to be different from codebook order
nesarc_sub$Sex <-
factor(nesarc_sub$Sex
, levels = c( 2
, 1
)
, labels = c( "Female"
, "Male"
)
)
# check ordering with a frequency table
table(nesarc_sub$Sex)
Female Male
24575 18518 nesarc_sub$TobaccoDependence <-
factor(nesarc_sub$TobaccoDependence
, levels = c( 1
, 0
)
, labels = c( "Nicotine Dependence"
, "No Nicotine Dependence"
)
)
# Remember to include the new category level from the original `NA` values.
nesarc_sub$SmokingFreq <-
factor(nesarc_sub$SmokingFreq
, levels = c( 7
, 6
, 5
, 4
, 3
, 2
, 1
)
, labels = c( "Never"
, "Once a month or less"
, "2 to 3 Days/month"
, "1 to 2 Days/week"
, "3 to 4 Days/week"
, "5 to 6 Days/week"
, "Every Day"
)
)
nesarc_sub$Ethnicity <-
factor(nesarc_sub$Ethnicity
# shorter labels are helpful for plots
, labels = c( "Cauc"
, "AfAm"
, "NaAm"
, "Asia"
, "Hisp"
)
# , labels = c("Caucasian"
# , "African American"
# , "Native American"
# , "Asian"
# , "Hispanic"
# )
)
nesarc_sub$SmokingFreq3 <-
factor(nesarc_sub$SmokingFreq3
, levels = c( 1
, 2
, 3
, 4
)
, labels = c( "Daily"
, "Weekly"
, "Monthly"
, "Never"
)
)
summary(nesarc_sub) ID Sex Age
1 : 1 Female:24575 Min. :18.0
2 : 1 Male :18518 1st Qu.:32.0
3 : 1 Median :44.0
4 : 1 Mean :46.4
5 : 1 3rd Qu.:59.0
6 : 1 Max. :98.0
(Other):43087
SmokingStatus TobaccoDependence
Smoked past 12 months : 9913 Nicotine Dependence : 4962
Smoked prior to 12 months: 8078 No Nicotine Dependence:38131
Never smoked :25080
NA's : 22
SmokingFreq DailyCigsSmoked Ethnicity
Never :25080 Min. : 0.000 Cauc:24507
Every Day :14836 1st Qu.: 0.000 AfAm: 8245
Once a month or less: 772 Median : 0.000 NaAm: 701
1 to 2 Days/week : 747 Mean : 6.462 Asia: 1332
3 to 4 Days/week : 687 3rd Qu.:10.000 Hisp: 8308
(Other) : 869 Max. :98.000
NA's : 102 NA's :262
Depression IncomePersonal Height_ft Height_in
No Depression :35254 Min. : 0 Min. :4.00 Min. : 0.000
Yes Depression: 7839 1st Qu.: 8800 1st Qu.:5.00 1st Qu.: 3.000
Median : 20000 Median :5.00 Median : 5.000
Mean : 28185 Mean :5.11 Mean : 5.279
3rd Qu.: 36000 3rd Qu.:5.00 3rd Qu.: 8.000
Max. :3000000 Max. :7.00 Max. :11.000
NA's :730 NA's :761
Weight_lb DaysSmoke Height_inches TotalCigsSmoked
Min. : 62.0 Min. : 0.00 Min. :48.0 Min. : 0.0
1st Qu.:140.0 1st Qu.: 0.00 1st Qu.:64.0 1st Qu.: 0.0
Median :165.0 Median : 0.00 Median :66.0 Median : 0.0
Mean :170.6 Mean :10.96 Mean :66.6 Mean : 187.5
3rd Qu.:192.0 3rd Qu.:30.00 3rd Qu.:70.0 3rd Qu.: 300.0
Max. :500.0 Max. :30.00 Max. :84.0 Max. :2940.0
NA's :1376 NA's :102 NA's :761 NA's :297
TotalCigsSmokedLog2 TotalCigsSmokedSqrt CigsSmokedFac SmokingFreq3
Min. : 0.000 Min. : 0.000 Q1 :25080 Daily :15983
1st Qu.: 0.000 1st Qu.: 0.000 Q2 : 0 Weekly : 1156
Median : 0.000 Median : 0.000 Q3 : 8630 Monthly: 772
Mean : 3.284 Mean : 7.887 Q4 : 9086 Never :25080
3rd Qu.: 8.234 3rd Qu.:17.321 NA's: 297 NA's : 102
Max. :11.522 Max. :54.222
NA's :297 NA's :297 3.1.8 Data subset rows to young smokers
In practice, I would probably subset the data much earlier. In this example I waited until the end to subset so that I could illustrate all of the other steps with complete values.
In our example, we are interested in the subset of people 25 or younger who smoked in the last 12 months and this is created using the filter() function. Note that SmokingStatus == "Smoked past 12 months" is used to choose people who smoked in the past 12 months.
Notice, the missing values have been removed (since NAs are not less than or equal to numbers).
# the subset command below will not include a row that evalutes to NA for the
# specified subset. For example:
(NA == 1)[1] NA# age and smoking subset
nesarc_sub <-
nesarc_sub %>%
filter(
Age <= 25 # people 25 or younger
, SmokingStatus == "Smoked past 12 months" # smoked in the past 12 months
)
dim(nesarc_sub)[1] 1706 20 ID Sex Age
21 : 1 Female:850 Min. :18.00
77 : 1 Male :856 1st Qu.:20.00
103 : 1 Median :22.00
122 : 1 Mean :21.61
136 : 1 3rd Qu.:23.00
150 : 1 Max. :25.00
(Other):1700
SmokingStatus TobaccoDependence
Smoked past 12 months :1706 Nicotine Dependence :896
Smoked prior to 12 months: 0 No Nicotine Dependence:810
Never smoked : 0
SmokingFreq DailyCigsSmoked Ethnicity
Every Day :1320 Min. : 1.00 Cauc:1060
3 to 4 Days/week : 91 1st Qu.: 5.00 AfAm: 211
1 to 2 Days/week : 88 Median :10.00 NaAm: 42
Once a month or less: 71 Mean :11.34 Asia: 58
5 to 6 Days/week : 68 3rd Qu.:20.00 Hisp: 335
(Other) : 65 Max. :98.00
NA's : 3 NA's :9
Depression IncomePersonal Height_ft Height_in
No Depression :1260 Min. : 0 Min. :4.000 Min. : 0.000
Yes Depression: 446 1st Qu.: 4000 1st Qu.:5.000 1st Qu.: 2.000
Median : 11000 Median :5.000 Median : 5.000
Mean : 13298 Mean :5.196 Mean : 5.425
3rd Qu.: 20000 3rd Qu.:5.000 3rd Qu.: 8.000
Max. :130000 Max. :6.000 Max. :11.000
NA's :14 NA's :14
Weight_lb DaysSmoke Height_inches TotalCigsSmoked
Min. : 85.0 Min. : 1.00 Min. :58.00 Min. : 1.0
1st Qu.:135.5 1st Qu.:30.00 1st Qu.:65.00 1st Qu.: 90.0
Median :160.0 Median :30.00 Median :68.00 Median : 300.0
Mean :165.8 Mean :25.33 Mean :67.78 Mean : 320.5
3rd Qu.:185.0 3rd Qu.:30.00 3rd Qu.:71.00 3rd Qu.: 600.0
Max. :390.0 Max. :30.00 Max. :81.00 Max. :2940.0
NA's :27 NA's :3 NA's :14 NA's :9
TotalCigsSmokedLog2 TotalCigsSmokedSqrt CigsSmokedFac SmokingFreq3
Min. : 1.000 Min. : 1.000 Q1 : 0 Daily :1479
1st Qu.: 6.508 1st Qu.: 9.487 Q2 : 0 Weekly : 153
Median : 8.234 Median :17.321 Q3 :1107 Monthly: 71
Mean : 7.448 Mean :15.987 Q4 : 590 Never : 0
3rd Qu.: 9.231 3rd Qu.:24.495 NA's: 9 NA's : 3
Max. :11.522 Max. :54.222
NA's :9 NA's :9 3.1.9 EXTRA: Working with Dates
3.1.9.1 NESARC age example
Dates can be tricky in R. The lubridate package makes them much easier.
Below I’m going to create a variable for the person’s age in years based on the difference between the person’s date of birth and the interview date. There is an AGE variable which I’m already using, but if there wasn’t, or I wanted to be more precise about their age (such as, 24.34 years old), then this is a way to do it.
I’m going to create a separate data frame just for this example to show how to do it; if you needed these variables, they should be part of your original subset above. Below, the “C” variables (CDAY) are the interview day, month, and year, and the “DOB” variables (DOBD) are the date of birth day, month, and year.
'data.frame': 43093 obs. of 6 variables:
$ CDAY : num 14 12 23 9 18 7 1 18 9 16 ...
$ CMON : num 8 1 11 9 10 9 9 9 1 8 ...
$ CYEAR: num 2001 2002 2001 2001 2001 ...
$ DOBD : num 27 15 28 27 11 28 16 10 7 25 ...
$ DOBM : num 9 12 12 6 9 12 3 11 4 3 ...
$ DOBY : num 1977 1973 1919 1983 1965 ...Combine the date fields together, then interpret them into a date object. I’ll do this for both the interview date and the date of birth. I print the head (first 6 observations) after each operation.
## interview date
# combine day, month, and year into one text string
cdate_text <- paste(dat_sub$CDAY, dat_sub$CMON, dat_sub$CYEAR)
head(cdate_text)[1] "14 8 2001" "12 1 2002" "23 11 2001" "9 9 2001" "18 10 2001"
[6] "7 9 2001" # create date object by interpretting the day, month, year with dmy()
cdate_date <- dmy(cdate_text)
head(cdate_date)[1] "2001-08-14" "2002-01-12" "2001-11-23" "2001-09-09" "2001-10-18"
[6] "2001-09-07"## date of birth
# combine day, month, and year into one text string
dob_text <- paste(dat_sub$DOBD, dat_sub$DOBM, dat_sub$DOBY)
head(dob_text)[1] "27 9 1977" "15 12 1973" "28 12 1919" "27 6 1983" "11 9 1965"
[6] "28 12 1966"# create date object by interpretting the day, month, year with dmy()
# if any failed to parse, it's because a day, month, or year is out of range
# such as "99" for missing values or refused to answer
# in these cases, the date is NA (missing), which is what we want
dob_date <- dmy(dob_text)
head(dob_date)[1] "1977-09-27" "1973-12-15" "1919-12-28" "1983-06-27" "1965-09-11"
[6] "1966-12-28"Now that we have the interview date and date of birth as date objects, we calculate the difference to give the person’s age at the time of the interview. Because the difference of dates is in the unit “days”, we convert to years.
# difference in days / 365.25 = difference in years
age_days <- cdate_date - dob_date
head(age_days)Time differences in days
[1] 8722 10255 29916 6649 13186 12672[1] 23.87953 28.07666 81.90554 18.20397 36.10130 34.69405# difference between the age we calculated and the age in the original dataset
head(age_years - NESARC$AGE)[1] 0.87953457 0.07665982 0.90554415 0.20396988 0.10130048 0.69404517Add the new age variable to your data frame, drop the ones you don’t need anymore, and you’re ready to go!
'data.frame': 43093 obs. of 7 variables:
$ CDAY : num 14 12 23 9 18 7 1 18 9 16 ...
$ CMON : num 8 1 11 9 10 9 9 9 1 8 ...
$ CYEAR : num 2001 2002 2001 2001 2001 ...
$ DOBD : num 27 15 28 27 11 28 16 10 7 25 ...
$ DOBM : num 9 12 12 6 9 12 3 11 4 3 ...
$ DOBY : num 1977 1973 1919 1983 1965 ...
$ age_years: num 23.9 28.1 81.9 18.2 36.1 ...# drop the original ones you don't need anymore
# by using subset() and select = -VAR, it excludes the VAR column
dat_sub <-
dat_sub %>%
select(-c(CDAY, CMON, CYEAR, DOBD, DOBM, DOBY))
str(dat_sub)'data.frame': 43093 obs. of 1 variable:
$ age_years: num 23.9 28.1 81.9 18.2 36.1 ...Now you’ve got an age (or other date object) you can use in your project.
3.1.9.2 AddHealth4 age example
Dates can be tricky in R. The lubridate package makes them much easier.
Below I’m going to create a variable for the person’s age in years based on the difference between the person’s date of birth and the interview date. There is an AGE variable which I’m already using, but if there wasn’t, or I wanted to be more precise about their age (such as, 24.34 years old), then this is a way to do it.
I’m going to create a separate data frame just for this example to show how to do it; if you needed these variables, they should be part of your original subset above. Below, the “i” variables (iday4) are the interview day, month, and year, and the date of birth variables (H4OD1MandH4OD1Y`) are the date of birth month and year (there is no day, so we’ll use day = 15 for the middle of the month).
Combine the date fields together, then interpret them into a date object. I’ll do this for both the interview date and the date of birth. I print the head (first 6 observations) after each operation.
AH14 <-
AH14 %>%
mutate(
idate = paste(iday4, imonth4, iyear4)
, dob = paste(15, h4od1m, h4od1y)
)
head(select(AH14, idate, dob)) idate dob
1 NA NA NA 15 NA NA
2 6 5 2008 15 11 1976
3 NA NA NA 15 NA NA
4 22 5 2008 15 1 1976
5 NA NA NA 15 NA NA
6 NA NA NA 15 NA NA# create date object by interpretting the day, month, year with dmy()
# if any failed to parse, it's because a day, month, or year is out of range
# such as "99" for missing values or refused to answer
# in these cases, the date is NA (missing), which is what we want
AH14 <-
AH14 %>%
mutate(
idate = dmy(idate)
, dob = dmy(dob)
)
head(select(AH14, idate, dob)) idate dob
1 <NA> <NA>
2 2008-05-06 1976-11-15
3 <NA> <NA>
4 2008-05-22 1976-01-15
5 <NA> <NA>
6 <NA> <NA>Now that we have the interview date and date of birth as date objects, we calculate the difference to give the person’s age at the time of the interview. Because the difference of dates is in the unit “days”, we convert to years.
# difference in days / 365.25 = difference in years
AH14 <-
AH14 %>%
mutate(
age_days = idate - dob
, age_years = as.numeric(age_days / 365.25) # change from "difftime" to "numeric"
)
head(select(AH14, idate, dob, age_days, age_years)) idate dob age_days age_years
1 <NA> <NA> NA days NA
2 2008-05-06 1976-11-15 11495 days 31.47159
3 <NA> <NA> NA days NA
4 2008-05-22 1976-01-15 11816 days 32.35044
5 <NA> <NA> NA days NA
6 <NA> <NA> NA days NA Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
24.43 27.48 28.92 28.89 30.36 33.94 1390 Now you’ve got an age (or other date object) you can use in your project.
3.1.9.3 Converting 12-hour to 24-hour time
Let’s use a time from Add Health “Wave IV Section 7: Sleep Patterns” about what time people wake up. There are three fields for hour, minute, and am/pm, and we want to convert these three fields into a single time number.
Now I'm going to ask you about your sleep patterns.
1H. On the days you go to work, school or similar activities,
what time do you usually wake up?
H4SP1H num [Hour] 96=refused, 98=don't know
H4SP1M num [Minute] 96=refused, 98=don't know
H4SP1T cat [am/pm] 1=am, 2=pm, 6=refused, 8=don't knowThere are two steps to address.
- First, “PM” adds 12 hours.
- Second, 12 PM stays as 12 and 12 AM becomes 00.
AH14 <-
AH14 %>%
mutate(
# name the hour, minute, and am_pm variables
wake_hour = h4sp1h %>% as.numeric()
, wake_min = h4sp1m %>% as.numeric()
, wake_ampm = h4sp1t %>% as.numeric()
# code missing values labels as NAs
, wake_hour = replace(wake_hour, wake_hour %in% c(96, 98), NA)
, wake_min = replace(wake_min , wake_min %in% c(96, 98), NA)
, wake_ampm = replace(wake_ampm, wake_ampm %in% c( 6, 8), NA)
# 1=am, 2=pm
, wake_ampm = case_when(
wake_ampm == 1 ~ "AM"
, wake_ampm == 2 ~ "PM"
)
# convert hour to 24-hour time
, wake_hour24 = case_when(
# hour is 1-11 and AM, then hour is 1-11
(wake_hour %in% 1:11) & (wake_ampm == "AM") ~ wake_hour
# hour is 1-11 and PM, then hour adds 12 for 13-23
, (wake_hour %in% 1:11) & (wake_ampm == "PM") ~ wake_hour + 12
# hour is 12 and AM, then hour is 0 (midnight)
, (wake_hour == 12 ) & (wake_ampm == "AM") ~ 0
# hour is 12 and PM, then hour is 12 (noon)
, (wake_hour == 12 ) & (wake_ampm == "PM") ~ 12
)
# convert minutes from 0-59 to fractions of an hour
, wake_min_frac = wake_min / 60
# add the 24-hour hour and the fraction minute to get a value from 0 to 23.99
, wake_time = wake_hour24 + wake_min_frac
)
# Here's an example of the new variables
# This looks correct :)
AH14 %>%
select(wake_hour, wake_min, wake_ampm, wake_time) %>%
na.omit() %>%
arrange(wake_time) %>%
slice(seq(1, n(), length = 10) ) wake_hour wake_min wake_ampm wake_time
1 12 0 AM 0.00
2 5 0 AM 5.00
3 5 30 AM 5.50
4 6 0 AM 6.00
5 6 30 AM 6.50
6 6 45 AM 6.75
7 7 0 AM 7.00
8 7 30 AM 7.50
9 9 0 AM 9.00
10 11 30 PM 23.50Below we see what we expect, that most people wake between 5 and 8 AM.
library(ggplot2)
p <- ggplot(data = AH14, aes(x = wake_time))
p <- p + geom_histogram(boundary = 0, binwidth = 0.5)
p <- p + geom_rug()
p <- p + labs(title = "AddHealth Wave 4 Wake time")
p <- p + theme_bw()
print(p)
3.2 Data is complete
At this point, your dataset should be complete, coded, formatted, and ready for analysis.
3.2.1 Plot entire dataset, show missing values
Plot the data subset to show the extent of missing values.
ggplot_missing <- function(dat) {
# https://github.com/njtierney/neato/blob/master/R/ggplot_missing.R
# A function that plots missing data in ggplot. For a more updated version,
# see `vis_miss` from visdat - github.com/njtierney/visdat
# dat is a data frame
# return a ggplot of the missing data.
dat2 <-
dat %>%
is.na()
# create a column indicating which rows all have data (no missing)
NO_MISSING <-
!(rowSums(!dat2) == ncol(dat2))
dat2 <-
cbind(dat2, NO_MISSING) %>%
reshape2::melt()
p <- ggplot2::ggplot(data = dat2, aes(x = Var2, y = Var1))
p <- p + ggplot2::geom_raster(aes(fill = value)) # , alpha = 0.6
p <- p + ggplot2::scale_fill_grey(name = "",
labels = c("Present","Missing"))
p <- p + ggplot2::theme_minimal()
p <- p + ggplot2::theme(axis.text.x = element_text(angle=90, vjust=0, hjust=0))
p <- p + ggplot2::labs(x = "Variables in Dataset",
y = "Rows / observations")
#ggplot2::scale_y_continuous(expand = c(0,0))
p <- p + ggplot2::scale_y_reverse(expand = c(0,0), breaks = c(1, seq(0, 10000, by=20)))
return(p)
}
4 Graphing and Tabulating
4.1 Thursday ———
4.2 Categorical variables
4.3 Tables for categorical variables
Basic tables can be created using the functions table or xtabs. Frequency tables are created for the variables TobaccoDependence, CigsSmokedFac, and SmokingFreq. When only looking at one variable, table is the same as xtabs, but xtabs will make looking at two-way, three-way, and higher tables possible.
Nicotine Dependence No Nicotine Dependence
896 810 # proportions available by passing these frequencies to prop.table()
table(nesarc_sub$TobaccoDependence) %>% prop.table()
Nicotine Dependence No Nicotine Dependence
0.5252052 0.4747948 # xtabs() can be used for one or more variables
T1 <- xtabs( ~ TobaccoDependence, data = nesarc_sub)
T1TobaccoDependence
Nicotine Dependence No Nicotine Dependence
896 810 TobaccoDependence
Nicotine Dependence No Nicotine Dependence
0.5252052 0.4747948 SmokingFreq
Never Once a month or less 2 to 3 Days/month
0 71 65
1 to 2 Days/week 3 to 4 Days/week 5 to 6 Days/week
88 91 68
Every Day
1320 (Look at the .Rmd code to see that these numbers in the next paragraph are calculated directly from the numbers summarized in the table! This way, if the data changes, the text will be correctly updated.)
In the data frame nesarc_sub, there are 896 nicotine dependent subjects and 810 subjects that are not nicotine dependent. Most of the subjects in nesarc_sub are daily smokers (68) with the remainder distributed uniformly across the first five levels of SmokingFreq.
4.3.1 Graphing frequency tables
The barplots are all created with the package ggplot2. The barplots start with the defaults for the geom_bar and add more detail to the plot with each graph.
library(ggplot2)
p1 <- ggplot(data = nesarc_sub, aes(x = TobaccoDependence))
p1 <- p1 + geom_bar()
p1
p2 <- ggplot(data = nesarc_sub, aes(x = SmokingFreq))
p2 <- p2 + geom_bar()
p2 <- p2 + labs(x = "", y = "Total Number", title = "Smoking Frequency for Young Smoking Adults")
p2 <- p2 + theme_bw()
p2 <- p2 + theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust = 1))
p2
4.3.2 Removing NAs from factor axes
When a factor variable has NAs, then NA appears as a level with frequences (see bar plots above). As of writing this, ggplot() does not have a way to remove the NAs before plotting. Therefore, we have to do it manually.
The easiest and most transparent way is to use the drop_na() function below where we only plot the non-NA values. Compare the two plots below of the same data, but the second plot had the NAs removed before plotting.
library(ggplot2)
# Original plot with NAs
p1 <- ggplot(data = nesarc_sub, aes(x = SmokingFreq))
p1 <- p1 + geom_bar()
p1 <- p1 + labs(title = "Basic plot")
p1 <- p1 + theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust = 1))
#p1
# subset excluded the NAs for the variable being plotted
# "subset() specifies the values to keep, and "!is.na()" means "keep the non-NA values"
# p2 <- ggplot(data = subset(nesarc_sub, !is.na(CigsSmokedFac)), aes(x = CigsSmokedFac))
# tidyverse syntax uses drop_na(Var) to drop observations where Var is NA
p2 <- ggplot(data = nesarc_sub %>% drop_na(SmokingFreq), aes(x = SmokingFreq))
p2 <- p2 + geom_bar()
p2 <- p2 + labs(title = "Using drop_na()")
p2 <- p2 + theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust = 1))
#p2
# grid.arrange() is a way to arrange several ggplot objects
library(gridExtra)
grid.arrange(grobs = list(p1, p2), nrow=1, top = "Excluding NAs from factor variable in ggplot")
4.3.3 Final versions
In your project, you only need to present your “Final version”. I’m illustrating the process above.
library(ggplot2)
p1 <- ggplot(data = nesarc_sub, aes(x = TobaccoDependence))
p1 <- p1 + geom_bar()
p1 <- p1 + labs(x = ""
, y = "Total Number"
, title = "Nicotine Dependence for Young Smoking Adults"
)
print(p1)
Interpretation: Among smokers 25 years and younger, there are slightly more (5% more) people who are nicotine dependent (52.5%) compared to those who are not nicotine dependent (47.5%).
p2 <- ggplot(data = nesarc_sub %>% drop_na(SmokingFreq), aes(x = SmokingFreq))
p2 <- p2 + geom_bar()
p2 <- p2 + labs(x = ""
, y = "Total Number"
, title = "Smoking Frequency for Young Smoking Adults"
)
p2 <- p2 + theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust = 1))
print(p2)
Interpretation: Among smokers 25 years and younger, most smoke every day (77.5%) while the less frequent categories are roughly evenly distributed around 4% to 5%.
4.4 Numeric variables
4.5 Graphing numeric variables
One popular way to graph a continuous variable is to use a histogram. R has the base function hist() which can be used for this purpose. We will also use the package ggplot2 to create histograms. We start with a basic histogram of the variable TotalCigsSmoked.
Experiment with the code in the next code chunk to set the binwidth = (extremely important!), decide if you need boundary = 0, add a title, and if needed the labels on the \(x\)- and \(y\)-axes.
p <- ggplot(data = nesarc_sub, aes(x = TotalCigsSmoked))
p <- p + geom_histogram(binwidth = 200)
p <- p + geom_rug()
p <- p + labs(x = "Estimated Cigarettes Smoked per Month")
p <- p + theme_bw()
p
# Use "boundary = 0" for values that can't be negative
p <- ggplot(data = nesarc_sub, aes(x = TotalCigsSmoked))
p <- p + geom_histogram(boundary = 0, binwidth = 100)
p <- p + geom_rug()
p <- p + labs(x = "Estimated Cigarettes Smoked per Month", title = "boundary at 0")
p <- p + theme_bw()
p
4.5.1 Creating Density Plots
A smoothed density curve can be added with the geom_density() function, where adjust = 2.0 is a fairly smooth curve. Experiment with adjust values such as 1.0, 1.5, 2.0, and 3.0 to decrease or increase the smoothness and decide what level best captures the general shape of your distribution.
p <- ggplot(data = nesarc_sub, aes(x = TotalCigsSmoked))
p <- p + geom_histogram(aes(y = ..density..), boundary = 0, binwidth = 100)
p <- p + geom_rug()
p <- p + geom_density(alpha = 0.2, fill = "gray50", adjust = 1.5)
p <- p + labs(x = "Estimated Cigarettes Smoked per Month")
p <- p + theme_bw()
p
4.5.2 Shape, center, and spread
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1.0 90.0 300.0 320.5 600.0 2940.0 9 [1] 1 90 300 600 2940[1] 300[1] 510The “shape” of the distribution for the estimated cigarettes smoked per month is skewed to the right. The “center” (median) of the distribution is 300 and the “spread” (interquartile range, middle 50%) for the distribution is 510.
4.5.3 Final versions
In your project, you only need to present your “Final version”. I’m illustrating the process above.
p <- ggplot(data = nesarc_sub, aes(x = TotalCigsSmoked))
p <- p + geom_histogram(aes(y = ..density..), boundary = 0, binwidth = 100)
p <- p + geom_rug()
p <- p + geom_density(alpha = 0.2, fill = "gray50", adjust = 1.5)
p <- p + labs(x = "Estimated Cigarettes Smoked per Month"
, y = ""
, title = "Total cigaretts smoked for Young Smoking Adults"
)
p <- p + theme_bw()
p
Interpretation: Among smokers 25 years and younger, the number of cigarettes smoked is right skewed with a median (center) of 300 and the IQR (spread) (interquartile range, middle 50%) for the distribution is 510. There are three distinct modes (peaks) in the distribution because of how we constructed this variable from two others, and because of the common numbers of cigarettes being relative to a pack size (20). There are several outlying values greater than 1200 (representing 40 cigarettes per day) that are possibly overestimates from the people responding.
p <- ggplot(data = nesarc_sub, aes(x = DailyCigsSmoked))
p <- p + geom_histogram(aes(y = ..density..), boundary = 0, binwidth = 2)
p <- p + geom_rug()
p <- p + geom_density(alpha = 0.2, fill = "gray50", adjust = 1.5)
p <- p + labs(x = "Estimated Cigarettes Smoked per Month"
, y = ""
, title = "Monthly cigaretts smoked for Young Smoking Adults"
)
p <- p + theme_bw()
p
Interpretation: Among smokers 25 years and younger, most smoke less than 10 cigarettes per day with another peak at 20 (1 pack per day), another small peak at 40 (2 packs per day), and some extreme outlying values at 60, 80, and 100 (3, 4, and 5 packs per day).
4.6 Week04: Bivariate graphs
4.6.1 Tuesday ———
4.6.2 Scatter plot (for regression): x = numerical, y = numerical
Research question: With my questions, I don’t have a meaningful question to ask that will relate two numeric variable.
An exploratory question whether there is a relationship between a person’s age (Age) and total number of cigaretts smoked (TotalCigsSmoked). Because age is an integer varible, I jitter the points, and I reduce the opacity (increase the transparancy) with alpha = 1/4 to see the density of points better (1/4 means that it takes 4 overlayed points to be fully opaque).
Interpretation: I added a linear regression smoother to see whether there was any linear trend with age. The horizontal line suggests that there isn’t a (linear) relationship here.
The first plot has all the points in their original locations, but they end up stacking on top of each other so you can’t tell how many points are there.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p <- p + geom_point(alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + labs(title = "Total Cigarettes Smoked vs Age")
print(p)
This is the same data as above, but I have added some horizontal jitter (displacement) to each point which doesn’t affect their value on the y-axis but allows us to see how many points are at each discrete age year.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p <- p + geom_jitter(position = position_jitter(width = 0.2), alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + labs(title = "Total Cigarettes Smoked vs Age")
print(p)
Here’s an example where there’s a strong positive relationship between the two numeric variables. Most variables in our dataset won’t exhibit such a nice relationship.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Height_inches, y = Weight_lb))
p <- p + geom_point(alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + labs(title = "Weight vs Height")
print(p)
If you have two groups you want to compare, you can create facets (small multiples) and show the relationship by each group. Below, I compare by sex. The slopes appear similar; later in the semester we’ll learn how to formally compare these lines.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Height_inches, y = Weight_lb))
p <- p + geom_point(alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + labs(title = "Weight vs Height")
p <- p + facet_wrap(~ Sex)
print(p)
4.6.3 Box plots (for ANOVA): x = categorical, y = numerical
Research question: 1. Is there a relationship between depression (Depression) and total number of cigarettes smoked (TotalCigsSmoked)?
The \(x\) variable should be a factor variable, and \(y\) should be numeric.
Interpretation: The boxplots indicate that the five number summary is about the same for both depression groups. Also, the means are about the same, though slightly higher for the people with depression.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Depression, y = TotalCigsSmoked))
p <- p + geom_boxplot(width = 0.5, alpha = 0.5)
p <- p + geom_jitter(position = position_jitter(width = 0.1), alpha = 1/4)
# diamond at mean for each group
p <- p + stat_summary(fun.y = mean, geom = "point", shape = 18, size = 6,
colour = "red", alpha = 0.8)
p <- p + labs(title = "Depression does not predict Total Cigarettes Smoked")
print(p)
4.6.4 Thursday ———
4.6.5 Mosaic plot or bivariate bar plots (for contingency tables): x = categorical, y = categorical
Research question: Is there a relationship between depression status (Depression) and nicotine dependence (TobaccoDependence)?
This strategy uses the base graphics by plotting the tables directly as bar plots. While I’m illustrating how frequency and proportion tables can be plotted, I do not recommend using the barplot() for your work.
Depression
TobaccoDependence No Depression Yes Depression
Nicotine Dependence 574 322
No Nicotine Dependence 686 124
Depression
TobaccoDependence No Depression Yes Depression
Nicotine Dependence 0.4555556 0.7219731
No Nicotine Dependence 0.5444444 0.2780269
The last graph needs a legend as well as a title. While it is possible to construct a legend and title for the last graph, it is much better to use an approach that will generate the legend automatically.
Interpretation: If a person has depression they are more likely to have nicotine dependence than if they do not have depression.
4.6.5.1 Categorical plots in ggplot
This set of plots using ggplot() are preferred.
Is there a relationship between depression status (Depression) and nicotine dependence (TobaccoDependence)?
Interpretation: If a person has depression they are more likely to have nicotine dependence than if they do not have depression.
library(ggplot2)
p <- ggplot(data = nesarc_sub, aes(x = Depression, fill = TobaccoDependence))
p <- p + geom_bar(position = "fill")
p <- p + theme_bw()
p <- p + labs(x = "Depression status"
, y = "Proportion"
, title = "Fraction of young adult smokers with and without\n nicotine dependence by depression status"
)
# the legend title can be modified, if desired (try this line)
p <- p + scale_fill_discrete(name="Tobacco Addiction\nStatus")
print(p)
Reversing color order: It is sometimes preferable to reverse the order of your color-filled variable to read the more interesting proportions from the bottom of the plot up the proportion scale. Use the function forcats::fct_rev() on your fill= variable.
library(ggplot2)
p <- ggplot(data = nesarc_sub, aes(x = Depression, fill = forcats::fct_rev(TobaccoDependence)))
p <- p + geom_bar(position = "fill")
p <- p + theme_bw()
p <- p + labs(x = "Depression status"
, y = "Proportion"
, title = "Fraction of young adult smokers with and without\n nicotine dependence by depression status"
)
# the legend title can be modified, if desired (try this line)
p <- p + scale_fill_discrete(name="Tobacco Addiction\nStatus")
print(p)
Does the a relationship between depression status (Depression) and nicotine dependence (TobaccoDependence) vary by sex?
Interpretation: It appears that nicotine dependence is slightly higher for women with depression than for men with depression.
library(ggplot2)
p <- ggplot(data = nesarc_sub, aes(x = Depression, fill = TobaccoDependence))
p <- p + geom_bar(position = "fill")
p <- p + theme_bw()
p <- p + labs(x = "Depression status"
, y = "Proportion"
, title = "Fraction of young adult smokers with and without\n nicotine dependence by depression status"
)
# the legend title can be modified, if desired (try this line)
p <- p + scale_fill_discrete(name="Tobacco Addiction\nStatus")
# added a row facets by Sex
p <- p + facet_grid(Sex ~ .)
# tilted the x-axis labels
p <- p + theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
print(p)
Research question: Is there a relationship between smoking frequency (SmokingFreq) and nicotine dependence (TobaccoDependence)?
Also compare by sex.
library(ggplot2)
p <- ggplot(data = subset(nesarc_sub, !is.na(SmokingFreq)), aes(x = SmokingFreq, fill = TobaccoDependence))
p <- p + geom_bar(position = "fill")
p <- p + theme_bw()
# tilt axis labels
p <- p + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
p <- p + labs(x = "Smoking Frequency", y = "Percent")
print(p)
4.6.5.2 Using Mosaic Plots
Mosaic plots will be used later in the course when we calculate statistics to test whether a pair of categorical variables are associated. These are helpful because they color-code each cell based on whether they provide evidence of an association. These plots are difficult to customize, so the ggplot() versions above will generally be preferred to these.
Interpretation: If a person has depression they are more likely to have nicotine dependence than if they do not have depression.
It’s easy to look at such a fine resolution (too many variables) that it’s difficult to understand what’s happening.
4.6.5.2.1 Customizing Mosaic Plots
There are several solutions to dealing with the messy labels that overlap. The help line below will show the many options you have control over. Read through the comments of each code chunk, see the new options, and see the result.
library(vcd)
# help(labeling_border)
# Initial plot
mosaic( ~ TobaccoDependence + Depression, data = nesarc_sub)
# Abbreviate labels, specify the number of characters (try a few)
mosaic( ~ TobaccoDependence + Depression, data = nesarc_sub
, abbreviate_labs = 4)
# Rotate and align labels
mosaic( ~ TobaccoDependence + Depression, data = nesarc_sub
, rot_labels = c(0, 90, 0, 0)
, just_labels = "right" # right-justifies all the labels (top and side)
)
# Rotate and align labels only y-axis label
mosaic( ~ TobaccoDependence + Depression, data = nesarc_sub
, rot_labels = c(0, 90, 0, 0)
, just_labels = c("center", "right") # order is opposite as in first argument
)
# Rotate and align labels only y-axis label, and offset the variable name
mosaic( ~ TobaccoDependence + Depression, data = nesarc_sub
, rot_labels = c(0, 90, 0, 0)
, just_labels = c("center", "right") # order is opposite as in first argument
, offset_varnames = c(0,0,0,9) # offset the y-axis var name
)
Final version
Labeling this plot well is, in my opinion, a bit too hard to be worth it. Notice in my code chunk options (you’ll need to look at the .Rmd file), that I’ve included fig.height = 8, fig.width = 8 which makes the size of the plot in the output 8 inches high by 8 inches wide.
library(vcd)
mosaic( ~ TobaccoDependence + Depression, data = nesarc_sub
, shade = TRUE
, rot_labels = c(0, 90, 0, 0) # order is: top, right, bottom, left
, just_labels = c("center", "right") # order is opposite as in first argument
, offset_varnames = c(0, 0, 0, 9) # offset the y-axis var name
# a big left margin makes space for the plot labels
, margins = c(3, 3, 3, 14) # order is: top, right, bottom, left
, main = "Fraction of young adult smokers with and without\n nicotine dependence by depression status"
, labeling_args = list(set_varnames = c(TobaccoDependence = "Tobacco Addiction Status"
, Depression = "Depression status")
)
)
Interpretation: If a person has nicotine dependence they are more likely to have depression than if they do not have nicotine dependence.
4.6.6 Logistic scatter plot (for logistic regression): x = numerical, y = categorical (binary)
Research question: 1. Frequency and quantity of smoking (TotalCigsSmoked) are markedly imperfect indices for determining an individual’s probability of exhibiting nicotine dependence (TobaccoDependence) (this is true for other drugs as well).
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = TotalCigsSmoked, y = TobaccoDependence))
p <- p + geom_jitter(position = position_jitter(height = 0.1), alpha = 1/4)
print(p)
Interpretation: These two groups are quite balanced on the probability of tobacco dependence given the total number of cigarettes smoked, though they increase together slightly.
To overlay a curve, we need to convert our response variable \(y\) to be binary 0 or 1, where we let 1 indicate a “success” in the direction of interest.
4.6.6.1 Collapsing a variable to binary
Whether you have a numeric or categorical variable that you’d like to represent with two levels, you’ll need to convert either to a numeric binary variable (values of 0 or 1, only). The value of 1 indicates “success” or the feature you’re interested in (below, 1 = someone with tobacco dependence) and 0 is “failure” (not tobacco dependent). Note that “success” does not necessarily mean “good”.
4.6.6.1.1 Numeric to binary
It is convenient to collapse a continuous variable into a categorical variable (as I did in Week03 with the cut() command) or further collapse a categorical variable into a binary variable.
One strategy to create a binary variable is to use the ifelse() function.
Let’s create a binary Age variable; if 21 years and older, than legal to purchase alcohol, otherwise it is illegal. See how the two columns below correspond as we wanted them to. I do not recommend overwriting the original column, but create a new variable in the data.frame.
nesarc_sub$Alcohol_binary <- ifelse(nesarc_sub$Age >= 21, "Legal", "Illegal")
nesarc_sub$Alcohol_binary01 <- ifelse(nesarc_sub$Age >= 21, 1, 0)
str(nesarc_sub$Alcohol_binary) chr [1:1706] "Legal" "Legal" "Legal" "Legal" "Legal" "Legal" "Legal" ... Age Alcohol_binary Alcohol_binary01
1 25 Legal 1
2 21 Legal 1
3 24 Legal 1
4 23 Legal 1
5 21 Legal 1
6 23 Legal 1
7 21 Legal 1
8 19 Illegal 0
9 19 Illegal 0
10 19 Illegal 0
11 18 Illegal 0
12 24 Legal 14.6.6.1.2 Categorical to binary
Two-level categorical to binary.
Nicotine Dependence No Nicotine Dependence
896 810 # if "Nicotine Dependence", then code it as a 1, otherwise code as a 0
nesarc_sub$TobaccoDependence01 <- ifelse(nesarc_sub$TobaccoDependence == "Nicotine Dependence", 1, 0)
# check that the frequencies are correct before and after
table(nesarc_sub$TobaccoDependence01)
0 1
810 896 # Note, sometimes it is necessary to make your variable a character type first.
# If the code above doesn't produce both 0s and 1s, then try this:
nesarc_sub$TobaccoDependence01 <- ifelse(as.character(nesarc_sub$TobaccoDependence) == "Nicotine Dependence", 1, 0)Multi-level categorical to binary.
# Smoking frequency has 7 categories
# if at least 3 Days/week, then code as a 1, otherwise code as a 0
table(nesarc_sub$SmokingFreq)
Never Once a month or less 2 to 3 Days/month
0 71 65
1 to 2 Days/week 3 to 4 Days/week 5 to 6 Days/week
88 91 68
Every Day
1320 nesarc_sub$SmokingFreq01 <-
ifelse(
nesarc_sub$SmokingFreq
%in%
c("3 to 4 Days/week"
, "5 to 6 Days/week"
, "Every Day"
)
, 1 # success
, 0 # failure
)
# check that the frequencies are correct before and after
table(nesarc_sub$SmokingFreq01)
0 1
227 1479 4.6.6.2 Plot with logistic curve
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = TotalCigsSmoked, y = TobaccoDependence01))
p <- p + geom_jitter(position = position_jitter(height = 0.1), alpha = 1/4)
# overlay a sigmodal curve to indicate direction of relationship
# As of ggplot2 2.1.0, need to use this binomial_smooth() function:
binomial_smooth <- function(...) {
geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}
p <- p + binomial_smooth() # old way: stat_smooth(method="glm", family="binomial", se=FALSE)
p <- p + labs(x = "Estimated Cigarettes Smoked per Month"
, y = "Tabacco Dependence (1=Yes, 0=No)"
, title = "Likelihood of Tabacco Dependence increases with number of cigs smoked"
)
print(p)
Interpretation: The probability of someone being tobacco dependent increases as the number of cigarettes smoked per month increases. The logistic curve model probably doesn’t apply at (or near) 0 cigarettes.
5 Statistical methods: inference and testing
5.1 Week05: Simple linear regression, logarithm transformation
5.1.1 Tuesday ———
In-class SLR intuition-building exercise and interpretation. Do not add to this document. Turn it in as its own document.
5.1.2 Thursday ———
5.1.3 1. Scatter plot, add a regression line.
Reproduce scatterplot from before, with a regression line.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p <- p + geom_jitter(position = position_jitter(width = 0.1), alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + labs(title = "Total Cigarettes Smoked vs Age")
print(p)
5.1.4 2. Fit the linear regression, interpret slope and \(R^2\) (R-squared) values
# fit the simple linear regression model
lm_TCS_Age <- lm(TotalCigsSmoked ~ Age, data = nesarc_sub)
# use summary() to parameters estimates (slope, intercept) and other summaries
summary(lm_TCS_Age)
Call:
lm(formula = TotalCigsSmoked ~ Age, data = nesarc_sub)
Residuals:
Min 1Q Median 3Q Max
-324.74 -224.96 -22.66 274.26 2623.50
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 287.263 66.339 4.330 1.58e-05 ***
Age 1.539 3.054 0.504 0.614
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.3 on 1695 degrees of freedom
(9 observations deleted due to missingness)
Multiple R-squared: 0.0001498, Adjusted R-squared: -0.0004401
F-statistic: 0.254 on 1 and 1695 DF, p-value: 0.6144The slope is 1.54. Thus, for every year increase in age, we expect an increase of 1.54 in total cigarettes smoked.
The \(R^2\) value is 1.498065210^{-4}, which is effectively 0. Compared to a model which uses only the mean for total cigarettes smoked, the model with age accounts for 1.510^{-4} additional variance.
When the \(x\) and \(y\) variables are correlated (\(r \ne 0\)), \(R^2\) will be larger (since \(R^2 = r^2\), the capitalization doesn’t matter in this case).
5.1.5 3. Interpret the intercept. Does it make sense?
The intercept of 287 is the predicted number of total cigarettes smoked for a person whose age is 0. Not only is that a gross extrapolation, but it doesn’t make sense.
5.1.6 4. Try plotting on log scale (\(x\)-only, \(y\)-only, both).
What happened to your data when you transformed it?
The plot with \((x, \log_{10}(y))\) looks much better. That’s because total cigarettes smoked (TotalCigsSmoked) has a lower bound of 0 and is right skewed (most values are low, some are extremely high).
library(ggplot2)
p1 <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p1 <- p1+ geom_jitter(position = position_jitter(width = 0.1), alpha = 1/4)
p1 <- p1+ stat_smooth(method = lm)
p1 <- p1+ labs(title = "Total Cigarettes Smoked vs Age")
#print(p1)
library(ggplot2)
p2 <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p2 <- p2+ geom_jitter(position = position_jitter(width = 0.1), alpha = 1/4)
p2 <- p2+ stat_smooth(method = lm)
p2 <- p2+ scale_y_log10()
p2 <- p2+ labs(title = "log10(Total Cigarettes Smoked) vs Age")
p2 <- p2+ labs(y = "log10(TotalCigsSmoked)")
#print(p2)
library(gridExtra)
grid.arrange(grobs = list(p1, p2), nrow=1, top = "Note: log10(Total Cigarettes Smoked) improves the symmetry about the regression line")
The Age range is limited and not an interesting variable to transform in this case. Below I provide those two plots: \((\log_{10}(x), y)\) and \((\log_{10}(x), \log_{10}(y))\). You should try them on your data. They will be most helpful for a right-skewed \(x\) variable, such as income.
library(ggplot2)
p1 <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p1 <- p1 + geom_jitter(position = position_jitter(width = 0.002), alpha = 1/4)
p1 <- p1 + stat_smooth(method = lm)
p1 <- p1 + scale_x_log10()
p1 <- p1 + labs(title = "Total Cigarettes Smoked vs log10(Age)")
#print(p1)
library(ggplot2)
p2 <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p2 <- p2 + geom_jitter(position = position_jitter(width = 0.002), alpha = 1/4)
p2 <- p2 + stat_smooth(method = lm)
p2 <- p2 + scale_x_log10()
p2 <- p2 + scale_y_log10()
p2 <- p2 + labs(title = "log10(Total Cigarettes Smoked) vs log10(Age)")
#print(p2)
library(gridExtra)
#grid.arrange(p1, p2, nrow=1, main = "log10(Age) isn't a good idea")
grid.arrange(grobs = list(p1, p2), nrow=1, top = "Note: log10(Age) isn't a good idea")
5.1.7 5. Does log transformation help?
The plot with \((x, \log_{10}(y))\) is an improvement over \((x,y)\) since the points are more symmetric around the regression line (they are not perfectly symmetric, but they are better). There might be a better transformation (such as square-root), or we might have to live with a non-symmetric response.
nesarc_sub$TotalCigsSmokedLog10 <- log10(nesarc_sub$TotalCigsSmoked)
# fit the simple linear regression model
lm_TCSL10_Age <- lm(TotalCigsSmokedLog10 ~ Age, data = nesarc_sub)
# use summary() to parameters estimates (slope, intercept) and other summaries
summary(lm_TCSL10_Age)
Call:
lm(formula = TotalCigsSmokedLog10 ~ Age, data = nesarc_sub)
Residuals:
Min 1Q Median 3Q Max
-2.2346 -0.2592 0.2426 0.5352 1.2296
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.319009 0.160521 14.447 <2e-16 ***
Age -0.004223 0.007389 -0.571 0.568
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6637 on 1695 degrees of freedom
(9 observations deleted due to missingness)
Multiple R-squared: 0.0001926, Adjusted R-squared: -0.0003972
F-statistic: 0.3266 on 1 and 1695 DF, p-value: 0.5678The slope is -0.00422. Thus, for every year increase in age, we expect an increase of -0.00422 in log 10 of total cigarettes smoked; that is, we expect a small decrease on the log scale.
Notice that the interpretation is on the scale of the regression. If the \(x\) or \(y\) variable is on a different scale (such as \(\log_{10}\)), then the interpretation is on that scale, too.
5.2 Week06: Correlation and Categorical contingency tables
5.2.1 Tuesday ———
5.2.2 Correlation
From a previous scatterplot, calculate the correlation and interpret it. The plot from the previous section and correlation are calculated below. (Unfortunately, my example is not very interesting.)
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Age, y = TotalCigsSmoked))
p <- p + geom_jitter(position = position_jitter(width = 0.1), alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + scale_y_log10()
p <- p + labs(title = "log10(Total Cigarettes Smoked) vs Age")
p <- p+ labs(y = "log10(TotalCigsSmoked)")
print(p)
# Calculate correlation
# Because of possible missing values, the "use" command excludes NAs before calculation
cor_TL_A <- cor(nesarc_sub$TotalCigsSmokedLog10, nesarc_sub$Age, use="pairwise.complete.obs")
cor_TL_A[1] -0.01387897Interpretation of correlation
The correlation is -0.0139, which is negative (but effectively 0). This suggests that there effectively no linear relationship between the \(\log_{10}\)(Total Cigarettes Smoked) and Age.
(See the Rmd file for how I included the correlation value in the text.)
Height and weight example
As a second example with a beautiful correlation, here is height and weight. We should all be so lucky to be analyzing data that looks like this!
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Height_inches, y = Weight_lb))
p <- p + geom_point(alpha = 1/4)
p <- p + stat_smooth(method = lm)
p <- p + labs(title = "Weight vs Height")
print(p)
# Calculate correlation
# Because of possible missing values, the "use" command excludes NAs before calculation
cor_H_W <- cor(nesarc_sub$Height_inches, nesarc_sub$Weight_lb, use="pairwise.complete.obs")
cor_H_W[1] 0.4996264Interpretation of correlation
The correlation is 0.5, which is positive and moderately strong. This suggests that there is a decent linear relationship between height and weight.
5.2.3 Thursday ———
5.2.4 Contingency table
From a previous multiple-variable categorical table, calculate conditional proportions and interpret them.
Research question: Is there a relationship between depression status (Depression) and nicotine dependence (TobaccoDependence), and does it differ by Sex?
(We’ll perform the statistical test for this in a couple weeks; right now we’re calculating, plotting, and interpretting the proportions.)
I’ll modify a plot from before using code from Worksheet 12.
5.2.5 TobaccoDependence by Depression and Sex
# proportions are of the last group_by() variable within combinations of the other variables
tab_long <-
nesarc_sub %>%
group_by(Sex, Depression, TobaccoDependence) %>%
summarise(Frequency = n()) %>%
mutate(Proportion = Frequency / sum(Frequency)) %>%
ungroup()
tab_long# A tibble: 8 x 5
Sex Depression TobaccoDependence Frequency Proportion
<fct> <fct> <fct> <int> <dbl>
1 Female No Depression Nicotine Dependence 263 0.464
2 Female No Depression No Nicotine Dependence 304 0.536
3 Female Yes Depression Nicotine Dependence 194 0.686
4 Female Yes Depression No Nicotine Dependence 89 0.314
5 Male No Depression Nicotine Dependence 311 0.449
6 Male No Depression No Nicotine Dependence 382 0.551
7 Male Yes Depression Nicotine Dependence 128 0.785
8 Male Yes Depression No Nicotine Dependence 35 0.2155.2.6 Interpretation of conditional proportions
The proportions above show that Females without depression have less nicotine dependence (0.464) than those with depression (0.686), and this disparity is even greater for Males (0.449 vs 0.785).
(See the Rmd file for how I included the values from the table in the text.)
# lines are sometimes easier, especially when many categories along the x-axis
library(ggplot2)
p <- ggplot(data = tab_long, aes(x = Depression, y = Proportion, colour = TobaccoDependence))
p <- p + geom_hline(yintercept = c(0, 1), alpha = 1/4)
p <- p + geom_point(aes(shape = TobaccoDependence), size = 3)
p <- p + geom_line(aes(linetype = TobaccoDependence, group = TobaccoDependence), size = 1)
p <- p + scale_y_continuous(limits = c(0, 1))
p <- p + facet_grid(. ~ Sex)
p <- p + theme_bw()
p <- p + theme(legend.position = "bottom")
p <- p + labs(title = "Proportion of Tobacco Dependence by Depression and Sex")
p <- p + labs(y = "Proportion of Tobacco Dependence")
print(p)
5.3 Week07: Inference and Parameter estimation (one-sample)
5.3.1 Tuesday ———
5.3.2 Numeric variable confidence interval for mean \(\mu\)
Using a numerical variable, calculate and interpret a confidence interval for the population mean.
# calcuate mean and CI
t_result <- t.test(nesarc_sub$TotalCigsSmokedLog10, conf.level = 0.95)
t_result
One Sample t-test
data: nesarc_sub$TotalCigsSmokedLog10
t = 138.31, df = 1696, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
2.196149 2.259333
sample estimates:
mean of x
2.227741 The sample mean of the \(\log_{10}(\)Total Cigs Smoked\()\) is 2.23. We are 95% confident that the population mean is between 2.2 and 2.26.
Yes, the interval can be very narrow when the sample size is large.
(Note that this does not mean that there’s 95% probability that the population mean is in that interval. Rather, for 95% of samples, the true population mean will be in the interval. We do not know whether the true mean is in our interval or not.)
# plot data and create table
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = TotalCigsSmokedLog10))
p <- p + geom_histogram(binwidth = 0.2, alpha = 1)
p <- p + geom_rug(alpha = 1/8)
# est
p <- p + geom_vline(aes(xintercept = as.numeric(t_result$estimate)), colour = "blue", size = 1)
p <- p + geom_rect(aes(xmin = t_result$conf.int[1], xmax = t_result$conf.int[2], ymin = -20, ymax = 0)
, fill = "blue", alpha = 1)
p <- p + labs(title = "Log 10 Total Cigs Smoked\nblue = est +- 95% CI")
print(p)
# checking model assumptions, some slight skewness observed
bs_one_samp_dist(nesarc_sub %>% drop_na(TotalCigsSmokedLog10) %>% pull(TotalCigsSmokedLog10))
5.3.3 Thursday ———
5.3.4 Categorical variable confidence interval for proportion \(p\)
Using a two-level categorical variable, calculate and interpret a confidence interval for the population proportion.
Let’s estimate the population proportion of people with Depression for the subset of people we’ve sampled from.
(Because I have at least 5 “successes” and “failures”, and a large sample size (\(n=1706\)) I know the prop.test() and binom.test() will give the same result. Try it for your data. They will only differ with quite small samples, especially if the proportions are close to 0 or 1.)
# proportions are of the last group_by() variable within combinations of the other variables
tab_dep_long <-
nesarc_sub %>%
group_by(Depression) %>%
summarise(Frequency = n()) %>%
mutate(Proportion = Frequency / sum(Frequency)) %>%
ungroup()
tab_dep_long# A tibble: 2 x 3
Depression Frequency Proportion
<fct> <int> <dbl>
1 No Depression 1260 0.739
2 Yes Depression 446 0.261# prop.test() is an asymptotic (approximate) test for a binomial random variable
p_summary <-
prop.test(
x = tab_dep_long %>% filter(Depression == "Yes Depression") %>% pull(Frequency)
, n = sum(tab_dep_long$Frequency)
, conf.level = 0.95
)
p_summary
1-sample proportions test with continuity correction
data: tab_dep_long %>% filter(Depression == "Yes Depression") %>% pull(Frequency) out of sum(tab_dep_long$Frequency), null probability 0.5
X-squared = 387.44, df = 1, p-value < 2.2e-16
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.2408466 0.2831009
sample estimates:
p
0.2614302 Note how to add error bars using geom_errorbar(). First determine the CI bounds from the prop.test() previously, then set those as limits.
[1] "statistic" "parameter" "p.value" "estimate" "null.value"
[6] "conf.int" "alternative" "method" "data.name" # here's the confidence interval bounds (the attribute tells us this is a 95% interval)
p_summary$conf.int[1] 0.2408466 0.2831009
attr(,"conf.level")
[1] 0.95[1] 0.2408466[1] 0.2831009library(ggplot2)
p <- ggplot(data = tab_dep_long %>% filter(Depression == "Yes Depression"), aes(x = Depression, y = Proportion))
p <- p + geom_hline(yintercept = c(0, 1), alpha = 1/4)
p <- p + geom_bar(stat = "identity")
p <- p + geom_errorbar(aes(min = p_summary$conf.int[1], max = p_summary$conf.int[2]), width=0.25)
p <- p + scale_y_continuous(limits = c(0, 1))
# swap axes to take less vertical space
p <- p + coord_flip()
print(p)
The sample proportion of the true population proportion of people depressed is 0.261. We are 95% confident that the population proportion is between 0.241 and 0.283.
Yes, the interval can be very narrow when the sample size is large.
(Note that this does not mean that there’s 95% probability that the population proportion is in that interval. Rather, for 95% of samples, the true population proportion will be in the interval. We do not know whether the true mean is in our interval or not.)
5.4 Week08: Hypothesis testing (one- and two-sample)
5.4.1 Mechanics of a hypothesis test (review)
Set up the null and alternative hypotheses in words and notation.
- In words: ``The population mean for [what is being studied] is different from [value of \(\mu_0\)].’’ (Note that the statement in words is in terms of the alternative hypothesis.)
- In notation: \(H_0: \mu=\mu_0\) versus \(H_A: \mu \ne \mu_0\) (where \(\mu_0\) is specified by the context of the problem).
Choose the significance level of the test, such as \(\alpha=0.05\).
Compute the test statistic, such as \(t_{s} = \frac{\bar{Y}-\mu_0}{SE_{\bar{Y}}}\), where \(SE_{\bar{Y}}=s/\sqrt{n}\) is the standard error.
Determine the tail(s) of the sampling distribution where the \(p\)-value from the test statistic will be calculated (for example, both tails, right tail, or left tail). (Historically, we would compare the observed test statistic, \(t_{s}\), with the critical value \(t_{\textrm{crit}}=t_{\alpha/2}\) in the direction of the alternative hypothesis from the \(t\)-distribution table with degrees of freedom \(df = n-1\).)
State the conclusion in terms of the problem.
- Reject \(H_0\) in favor of \(H_A\) if \(p\textrm{-value} < \alpha\).
- Fail to reject \(H_0\) if \(p\textrm{-value} \ge \alpha\). (Note: We DO NOT accept \(H_0\).)
Check assumptions of the test (for now we skip this).
5.4.2 What do we do about “significance”?
Adapted from Significance Magazine.
Recent calls have been made to abandon the term “statistical significance”. The American Statistical Association (ASA) issued its statement and recommendation on p-values (see the special issue of p-values for more).
In summary, the problem of “significance” is one of misuse, misunderstanding, and misinterpretation. The recommendation in this class is that it is no longer sufficient to say that a result is “statistically significant” or “non-significant” depending on whether a p-value is less than a threshold. Instead, we will be looking for wording as in the following paragraph.
“The difference between the two groups turns out to be small (8%), while the probability (\(p\)) of observing a result at least as extreme as this under the null hypothesis of no difference between the two groups is \(p = 0.003\) (that is, 0.3%). This p-value is statistically significant as it is below our pre-defined threshold (\(p < 0.05\)). However, the p-value tells us only that the 8% difference between the two groups is somewhat unlikely given our hypothesis and null model’s assumptions. More research is required, or other considerations may be needed, to conclude that the difference is of practical importance and reproducible.”
5.4.3 Tuesday ———
5.4.4 Two-sample \(t\)-test: TotalCigsSmoked by Depression
- Using a numerical response variable and a two-level categorical variable (or a categorical variable you can reduce to two levels), specify a two-sample \(t\)-test associated with your research questions.
- (2 p) Specify the hypotheses in words and notation (either one- or two-sided test),
- (0 p) use
t.test()to calculate the mean, test statistic, and p-value, - (3 p) state the significance level, test statistic, and p-value, and
- (2 p) state the conclusion in the context of the problem.
- (1 p) Given your conclusion, could you have committed at Type-I or Type-II error?
- (2 p) Provide an appropriate plot of the data and sample estimates in a well-labelled plot.
## If we create a summary data.frame with a similar structure as our data, then we
## can annotate our plot with those summaries.
est_mean_TCS_D <-
nesarc_sub %>%
group_by(Depression) %>%
summarise(TotalCigsSmoked = mean(TotalCigsSmoked, na.rm = TRUE)) %>%
ungroup()
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Depression, y = TotalCigsSmoked))
p <- p + geom_boxplot(width = 0.5, alpha = 0.5)
p <- p + geom_jitter(position = position_jitter(width = 0.1), alpha = 1/4)
# diamond at mean for each group
p <- p + stat_summary(fun.y = mean, geom = "point", shape = 18, size = 6,
colour = "red", alpha = 0.8)
p <- p + labs(title = "Total Cigarettes Smoked by Depression")
print(p)
# histogram using ggplot
p <- ggplot(nesarc_sub, aes(x = TotalCigsSmoked))
p <- p + geom_histogram(binwidth = 40)
p <- p + geom_rug()
p <- p + geom_vline(data = est_mean_TCS_D, aes(xintercept = TotalCigsSmoked), colour = "red")
p <- p + facet_grid(Depression ~ .)
print(p)
Welch Two Sample t-test
data: TotalCigsSmoked by Depression
t = -1.8191, df = 733.22, p-value = 0.06931
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-59.186888 2.255475
sample estimates:
mean in group No Depression mean in group Yes Depression
313.0782 341.5439 “Is the population mean total cigarettes smoked different for those with depression or not?”
- \(H_0: \mu_{D} = \mu_{ND}\) versus \(H_A: \mu_{D} \ne \mu_{ND}\)
Let \(\alpha=0.05\), the significance level of the test and the Type-I error probability if the null hypothesis is true.
\(t_{s} = -1.819\).
\(p=0.0693\), this is the observed significance of the test.
Because \(p=0.0693 > 0.05\), we have insufficient evidence to reject \(H_0\), concluding that the total cigarettes smoked does not differ by depression status.
Because we failed to reject \(H_0\), I could have committed a Type-II error (that there is a difference but we didn’t detect it).
5.4.5 Thursday ——– Fall Break
Enjoy your Fall Break!
5.5 Week09: ANOVA and Assessing Assumptions
5.5.1 Tuesday ———
Do not add to this document. Turn it in as its own document.
5.5.2 Thursday ———
5.5.3 ANOVA: Total cigarettes smoked by Ethnicity
5.5.3.1 Transform the response variable to satisfy assumptions
I’ve decided to try to find a better transformation of the total cigarettes smoked, TotalCigsSmoked, to improves the symmetry (and normality) of the distributions (the original scale was right skewed, and the \(\log_{10}\) was left skewed).
The Box-Cox transformation is a power transformation where we find a best power \(\lambda\) for \(y^{\lambda}\) such that the residuals from the model are most normal. We’re fitting ANOVA, so we first fit the model, then have the boxcox() function estimate the power that gives us residuals that are most “normal”. Finally, we’ll create a new variable \(y^{\lambda}\) for a sensible value of \(\lambda\).
# Box-Cox transformation finds the best power transformation of the response
# to make the residuals most "normal".
# A transformation close to this will help us meet our normality assumption.
# It is often preferred to choose a transformation that is a "simple" exponent,
# such as 0.5, rather than the exact estimate of the best, such as 0.4242424.
library(car)
aov_summary <- aov(TotalCigsSmoked ~ Ethnicity, data = nesarc_sub)
bc_summary <- boxCox(aov_summary
, lambda = seq(-2, 2, 0.1)
, plotit = TRUE)
[1] 0.4242424The best power transformation value is \(\lambda = 0.424\), so I will choose 0.5 (which is a square root).
When I reassess the power after taking the square root, the power \(\lambda = 1\) is close to the 95% interval, so I will use this as my transformation to ease interpretation and because the plots we see below of the data on the square-root scale are close to normal.
## Choose a square root transformation.
nesarc_sub$TotalCigsSmokedSqrt <- sqrt(nesarc_sub$TotalCigsSmoked)
library(car)
aov_summary <- aov(TotalCigsSmokedSqrt ~ Ethnicity, data = nesarc_sub)
bc_summary <- boxCox(aov_summary
, lambda = seq(-2, 2, 0.1)
, plotit = TRUE
)
[1] 0.82828285.5.3.2 ANOVA Hypothesis test
A table of the ethnic groups to compare is below.
Cauc AfAm NaAm Asia Hisp
1060 211 42 58 335 Let \(\mu_j\) = pop mean number of the square-root total cigarettes smoked for the five ethnicities identified in the dataset: Caucasian, African American, Native American, Asian, and Hispanic, numbered \((j = 1, 2, 3, 4, 5)\), respectively. We wish to test \(H_0: \mu_1=\cdots=\mu_5\) against \(H_A: \textrm{not } H_0\) (at least one pair of means differ.
Plot the data in a way that compares the means. Error bars are 95% confidence intervals of the mean.
# Plot the data using ggplot
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = Ethnicity, y = TotalCigsSmokedSqrt))
# plot a reference line for the global mean (assuming no groups)
p <- p + geom_hline(yintercept = mean(nesarc_sub$TotalCigsSmokedSqrt),
colour = "black", linetype = "dashed", size = 0.3, alpha = 0.5)
# boxplot, size=.75 to stand out behind CI
p <- p + geom_violin(width = 0.5, alpha = 0.25)
p <- p + geom_boxplot(width = 0.25, alpha = 0.25)
# points for observed data
p <- p + geom_point(position = position_jitter(w = 0.05, h = 0), alpha = 0.2)
# diamond at mean for each group
p <- p + stat_summary(fun.y = mean, geom = "point", shape = 18, size = 4,
colour = "red", alpha = 0.8)
# confidence limits based on normal distribution
p <- p + stat_summary(fun.data = "mean_cl_normal", geom = "errorbar",
width = .2, colour = "red", alpha = 0.8)
p <- p + labs(title = "Total cigarettes smoked by Ethnicity")
p <- p + ylab("Square-root of Total Cigs Smoked")
print(p)
Hypothesis test
Set up the null and alternative hypotheses in words and notation.
- In words: “The population mean square-root total cigarettes smoked is different between ethnicities.”
- In notation: \(H_0: \mu_1=\cdots=\mu_5\) versus \(H_A: \textrm{not } H_0\) (at least one pair of means differ).
Let the significance level of the test be \(\alpha=0.05\).
Compute the test statistic.
Df Sum Sq Mean Sq F value Pr(>F)
Ethnicity 4 6188 1547.0 25.17 <2e-16 ***
Residuals 1692 104009 61.5
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
9 observations deleted due to missingnessThe \(F\)-statistic for the ANOVA is \(F = 25.2\).
- Compute the \(p\)-value from the test statistic.
The p-value for testing the null hypothesis is \(p = 2.83\times 10^{-20}\).
- State the conclusion in terms of the problem.
Because \(p = 2.83\times 10^{-20} < 0.05\), we reject \(H_0\) in favor of \(H_1\) concluding that the mean square-root total cigarettes smoked is different between ethnicities.
5.5.3.3 Check assumptions
- Check assumptions of the test.
- Residuals are normal
- Populations have equal variances.
Check whether residuals are normal.
Plot the residuals and assess whether they appear normal.
# Plot the data using ggplot
df_res <- data.frame(res = aov_summary$residuals)
library(ggplot2)
p <- ggplot(df_res, aes(x = res))
p <- p + geom_histogram(aes(y = ..density..), binwidth = 4)
p <- p + geom_density(colour = "blue")
p <- p + geom_rug()
p <- p + stat_function(fun = dnorm, colour = "red", args = list(mean = mean(df_res$res), sd = sd(df_res$res)))
p <- p + labs(title = "ANOVA Residuals")
print(p)
There are a few extreme observations, and a large peak of points in the center. They are not symmetric, but they are unimodal.
- Plot the residuals versus the normal quantiles. If the residuals are normal, then the will fall on the center line and very few will be outside the error bands.
# QQ plot
par(mfrow=c(1,1))
library(car)
qqPlot(aov_summary$residuals, las = 1, id = list(n = 0, cex = 1), lwd = 1, main="QQ Plot")
There are deviations from normality, and I think this largely has to do with rounding the number of cigarettes smoked. Those lowest values are a group of very small numbers (1 cig - 10 cig in a month), and because those are bounded at 0, we’re getting that the residuals are larger than expected under the assumption that they are normal. The points in the middle under the red line are a group of points of regular smokers who smoke a moderate amount.
- A formal test of normality on the residuals tests the hypothesis \(H_0:\) The distribution is Normal vs \(H_1:\) The distribution is not Normal. We can test the distribution of the residuals.
Three formal tests for normality are reported below.
Shapiro-Wilk normality test
data: aov_summary$residuals
W = 0.97776, p-value = 1.44e-15
Anderson-Darling normality test
data: aov_summary$residuals
A = 12.22, p-value < 2.2e-16
Cramer-von Mises normality test
data: aov_summary$residuals
W = 1.9661, p-value = 7.37e-10All three tests reject normality.
In summary, both the plots and tests for normality suggest these residuals are not normal. The original plots do indicate that for each ethnic group the square-root transformation has made the data roughly symmetric, but not normal. We have a few outliers we may want to worry about (remove and see if the results change). Ultimately, we will probably want to use a nonparametric method to test the ANOVA hypothesis, a method that doesn’t require the normality assumption.
Check whether populations have equal variances.
Look at the numerical summaries below.
# calculate summaries
nesarc_sub_summary <-
nesarc_sub %>%
select(Ethnicity, TotalCigsSmokedSqrt) %>%
group_by(Ethnicity) %>%
summarise(
m = mean(TotalCigsSmokedSqrt, na.rm = TRUE)
, s = sd(TotalCigsSmokedSqrt, na.rm = TRUE)
, n = length(TotalCigsSmokedSqrt)
) %>%
ungroup()
nesarc_sub_summary# A tibble: 5 x 4
Ethnicity m s n
<fct> <dbl> <dbl> <int>
1 Cauc 17.4 8.16 1060
2 AfAm 14.5 6.94 211
3 NaAm 15.9 7.64 42
4 Asia 14.5 5.98 58
5 Hisp 12.7 7.64 335The standard deviations are between 5.98 and 8.16, which is not a wide range.
- Formal tests for equal variances. We can test whether the variances are equal between our three groups. This is similar to the ANOVA hypothesis, but instead of testing means we’re tesing variances. \(H_0: \sigma^2_1=\sigma^2_2=\sigma^2_3\) versus \(H_A: \textrm{not } H_0\) (at least one pair of variances differ).
## Test equal variance
# assumes populations are normal
bartlett.test(TotalCigsSmokedSqrt ~ Ethnicity, data = nesarc_sub)
Bartlett test of homogeneity of variances
data: TotalCigsSmokedSqrt by Ethnicity
Bartlett's K-squared = 16.623, df = 4, p-value = 0.002288# does not assume normality, requires car package
library(car)
leveneTest(TotalCigsSmokedSqrt ~ Ethnicity, data = nesarc_sub)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 4.7731 0.0007871 ***
1692
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Fligner-Killeen test of homogeneity of variances
data: TotalCigsSmokedSqrt by Ethnicity
Fligner-Killeen:med chi-squared = 20.064, df = 4, p-value =
0.0004852Interpret the result of the appropriate test. If normality was reasonable then use Bartlett, otherwise use Levene.
Both the Levene and Fligner tests reject equal variance. However, we have a large sample size, so we are “over powered” to identify small deviations from equality of variance. While this difference probably wouldn’t violate our assumption so badly as to invalidate the test, we will still want to use the nonparametric test next week, instead.
5.5.3.4 Post Hoc pairwise comparison tests
- If the ANOVA null hypothesis was rejected, then perform follow-up Post Hoc pairwise comparison tests to determine which pairs of means are different.
(We’ll continue assuming our assumptions are met, for the sake of completing the exercise. We understand that we’ll want to use a nonparametric test, instead.)
There are several multiple comparison methods described in the notes. Let’s use Tukey’s Honest Significant Difference (HSD) here to test which pairs of populations differ.
# A tibble: 5 x 4
Ethnicity m s n
<fct> <dbl> <dbl> <int>
1 Hisp 12.7 7.64 335
2 Asia 14.5 5.98 58
3 AfAm 14.5 6.94 211
4 NaAm 15.9 7.64 42
5 Cauc 17.4 8.16 1060 Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = TotalCigsSmokedSqrt ~ Ethnicity, data = nesarc_sub)
$Ethnicity
diff lwr upr p adj
AfAm-Cauc -2.84376899 -4.461685 -1.22585334 0.0000171
NaAm-Cauc -1.44066312 -4.809443 1.92811662 0.7698046
Asia-Cauc -2.92956874 -5.817125 -0.04201297 0.0448071
Hisp-Cauc -4.66440504 -6.010290 -3.31851959 0.0000000
NaAm-AfAm 1.40310587 -2.215808 5.02201991 0.8275803
Asia-AfAm -0.08579975 -3.261621 3.09002144 0.9999933
Hisp-AfAm -1.82063605 -3.707239 0.06596660 0.0645798
Asia-NaAm -1.48890562 -5.826747 2.84893557 0.8823694
Hisp-NaAm -3.22374192 -6.729494 0.28200985 0.0886179
Hisp-Asia -1.73483630 -4.781079 1.31140678 0.5266047Summarize results by ordering the means and grouping pairs that do not differ (see notes for examples). I find it easiest to look for the large p-values (which pairs are not significantly different) and group those first, then look at the significant pairs and make sure none are connected.
Here’s two ways to do this (vertical is useful when there are many groups, horizontal is easier to read).
(numerical sumamries are a nice addition)
NotSigDiff Mean Ethnicity m s n
| 12.7 Hispanic 12.7 7.64 335
|| 14.5 Asian 14.5 6.94 211
|| 14.5 African American 14.5 5.98 58
|| 15.9 Native American 15.9 7.64 42
| 17.4 Caucasian 17.4 8.16 1060
Ethnicity: Hispanic Asian African American Native American Caucasian
Mean: 12.7 14.5 14.5 15.9 17.4
-------------------------------------
--------------------------------------------
----------------------------5.6 Week10: Nonparametric methods and Binomial and multinomial proportion tests
5.6.1 Tuesday ———
5.6.2 Thursday ———
5.6.3 Multinomial goodness-of-fit: Ethnicity at US census proportions
Question: Did NESARC sample at random from the US population so that ethnicities are proportionally represented?
To answer this question, I need to go back to the original NESARC dataset. If I use my subset, I may have removed observations that will bias my answer to this question.
5.6.3.1 Observed
NB. All the results in this section depend on nesarc_eth and US_eth the frequency tables of the full NESARC dataset and the 2010 US Census.
# create a temporary Ethnicity factor variable from the orginal dataset
nesarc_eth <-
NESARC %>%
select(ETHRACE2A) %>%
mutate(
Ethnicity = factor(ETHRACE2A
, labels = c("Caucasian"
, "African American"
, "Native American"
, "Asian"
, "Hispanic"
)
)
) %>%
select(Ethnicity) %>%
group_by(Ethnicity) %>%
summarise(obs = n()) %>%
mutate(obs_prop = obs / sum(obs)) %>%
ungroup()
nesarc_eth# A tibble: 5 x 3
Ethnicity obs obs_prop
<fct> <int> <dbl>
1 Caucasian 24507 0.569
2 African American 8245 0.191
3 Native American 701 0.0163
4 Asian 1332 0.0309
5 Hispanic 8308 0.193 5.6.3.2 Expected
The US Ethnicity proportions are listed below with numbered categories consistent with NESARC.
ETHRACE2A (Ethnicity)
IMPUTED RACE/ETHNICITY
Nominal
24507 1. White, Not Hispanic or Latino
8245 2. Black, Not Hispanic or Latino
701 3. American Indian/Alaska Native, Not Hispanic or Latino
1332 4. Asian/Native Hawaiian/Pacific Islander, Not Hispanic or Latino
8308 5. Hispanic or Latino
https://en.wikipedia.org/wiki/Demographics_of_the_United_States#Race_and_ethnicity
The U.S. population's distribution by race and ethnicity in 2010 was as follows:
1. Non-Hispanic White 196,817,552 63.7 %
2. Non-Hispanic Black or African American 37,685,848 12.2 %
3. Non-Hispanic American Indian or Alaska Native 2,247,098 0.7 %
4. Non-Hispanic Asian 14,465,124 4.7 %
4. Non-Hispanic Native Hawaiian or other Pacific Islander 481,576 0.2 %
4. Non-Hispanic some other race 604,265 0.2 %
4. Non-Hispanic two or more races 5,966,481 1.9 %
5. Hispanic or Latino 50,477,594 16.4 %Expected numbers and combining observed and expected into a table.
US_eth <- c(196817552
, 37685848
, 2247098
, 14465124 + 481576 + 604265 + 5966481
, 50477594
)
nesarc_eth <-
nesarc_eth %>%
mutate(
exp_prop = US_eth / sum(US_eth)
# Expected number of eth
, exp = sum(obs) * exp_prop
)
nesarc_eth# A tibble: 5 x 5
Ethnicity obs obs_prop exp_prop exp
<fct> <int> <dbl> <dbl> <dbl>
1 Caucasian 24507 0.569 0.637 27471.
2 African American 8245 0.191 0.122 5260.
3 Native American 701 0.0163 0.00728 314.
4 Asian 1332 0.0309 0.0697 3003.
5 Hispanic 8308 0.193 0.163 7045.5.6.4 Perform \(\chi^2\) Goodness-of-fit test
- Set up the null and alternative hypotheses in words and notation.
- In words: “The NESARC and US Census proportions are equal for people of each Ethnicity.”
- In notation: \(H_0: p_{1} = 0.637, p_{2} = 0.122, p_{3} = 0.007, p_{4} = 0.07, \textrm{and} p_{5} = 0.163\) versus \(H_A: \textrm{not } H_0\)
- Choose the significance level of the test, such as \(\alpha=0.05\).
# calculate chi-square goodness-of-fit test
x_summary <- chisq.test(nesarc_eth$obs, correct = FALSE, p = nesarc_eth$exp_prop)
# print result of test
x_summary
Chi-squared test for given probabilities
data: nesarc_eth$obs
X-squared = 3648.5, df = 4, p-value < 2.2e-16- Compute the test statistic.
The test statistic is \(X^2 = 3648.5\).
- Compute the \(p\)-value.
The p-value \(= 0\).
- (1 p) State the conclusion in terms of the problem.
- We reject \(H_0\) in favor of \(H_A\) if \(p\textrm{-value} < \alpha\).
- Fail to reject \(H_0\) if \(p\textrm{-value} \ge \alpha\). (Note: We DO NOT accept \(H_0\).)
Because the p-value \(= 0 < 0.05 = \alpha\), we reject \(H_0\) in favor of \(H_A\) concluding that the proportions of ethnicities sampled in NESARC differs from the US Census.
- (1 p) Check assumptions of the test (see notes).
Because there is an expected frequency of at least 5 for each category, the assumptions are met for this test.
5.6.4.1 Chi-sq statistic helps us understand the deviations
# use output in x_summary and create table
x_table <-
data.frame(
Ethnicity = nesarc_eth$Ethnicity
, obs = x_summary$observed
, exp = x_summary$expected
, res = x_summary$residuals
, chisq = x_summary$residuals^2
, stdres = x_summary$stdres
)
x_table Ethnicity obs exp res chisq stdres
1 Caucasian 24507 27470.7088 -17.88136 319.7431 -29.69830
2 African American 8245 5259.9829 41.15804 1693.9841 43.92609
3 Native American 701 313.6376 21.87276 478.4174 21.95279
4 Asian 1332 3003.2865 -30.49668 930.0473 -31.61836
5 Hispanic 8308 7045.3843 15.04246 226.2756 16.44690Plot observed vs expected values to help identify Ethnicity groups that deviate the most. Plot contribution to chi-square values to help identify Ethnicity groups that deviate the most. The term “Contribution to Chi-Square” (chisq) refers to the values of \(\frac{(O-E)^2}{E}\) for each category. \(\chi^2_{\textrm{s}}\) is the sum of those contributions.
These plots suggest that African Americans and Native Americans were oversampled and Asians were undersampled.
library(tidyr)
x_table_obsexp <-
x_table %>%
gather(
# list the source columns to combine into two columns
obs
, exp
# key is the name of the destination column identifying each
# original column that the measurement came from
, key = "stat"
# value is the column name for values in table
, value = "value"
)
# Observed vs Expected counts
library(ggplot2)
p1 <- ggplot(x_table_obsexp, aes(x = Ethnicity, fill = stat, weight=value))
p1 <- p1 + geom_bar(position="dodge")
p1 <- p1 + labs(title = "Observed and Expected frequencies")
p1 <- p1 + xlab("Ethnicity category")
print(p1)
# Contribution to chi-sq
# pull out only the Ethnicity and chisq columns
x_table_chisq <-
x_table %>%
select(
Ethnicity
, chisq
) %>%
# reorder the Ethnicity categories to be descending relative to the chisq statistic
arrange(
-chisq
)
x_table_chisq Ethnicity chisq
1 African American 1693.9841
2 Asian 930.0473
3 Native American 478.4174
4 Caucasian 319.7431
5 Hispanic 226.2756 # reorder Ethnicity descending by chisq statistic
p2 <- ggplot(x_table_chisq, aes(x = reorder(Ethnicity, -chisq), weight = chisq))
p2 <- p2 + geom_bar()
p2 <- p2 + labs(title = "Contribution to Chi-sq statistic")
p2 <- p2 + xlab("Sorted Ethnicity")
p2 <- p2 + ylab("Chi-sq")
print(p2)
5.6.4.2 Multiple Comparisons
The goodness-of-fit test suggests that at least one of the Ethnicity category proportions for the observed number og eth differs from the expected number. A reasonable next step in the analysis would be to separately test the five hypotheses to see which Ethnicity categories led to rejecting the goodness-of-fit null hypothesis:
- \(H_0: p_{1} = 0.637\)
- \(H_0: p_{2} = 0.122\)
- \(H_0: p_{3} = 0.007\)
- \(H_0: p_{4} = 0.07\)
- \(H_0: p_{5} = 0.163\)
A Bonferroni comparison with a Family Error Rate \(\alpha=0.05\) sets each test’s significance level to \(\alpha = 0.05 / 5 = 0.01\), which corresponds to 99% two-sided CIs for the individual proportions.
Below I perform exact binomial tests of proportion for each of the five Ethnicity categories at the Bonferroni-adjusted significance level. I save the p-values and confidence intervals in a table along with the observed and census proportions in order to display the table below.
b_sum <-
nesarc_eth %>%
group_by(Ethnicity) %>%
do(
b_out = binom.test(.$obs, sum(nesarc_eth$obs), p = .$exp_prop, alternative = "two.sided", conf.level = 1-0.05 / nrow(.))
) %>%
summarise(
p.value = b_out$p.value
, CI.lower = b_out$conf.int[1]
, CI.upper = b_out$conf.int[2]
) %>%
ungroup()
b_sum$Ethnicity <- nesarc_eth$Ethnicity
b_sum$Observed <- x_table$obs / sum(x_table$obs)
b_sum$exp_prop <- nesarc_eth$exp_prop
b_sum# A tibble: 5 x 6
p.value CI.lower CI.upper Ethnicity Observed exp_prop
<dbl> <dbl> <dbl> <fct> <dbl> <dbl>
1 2.25e-189 0.564 0.573 Caucasian 0.569 0.637
2 Inf.Nae-324 0.188 0.195 African American 0.191 0.122
3 1.51e- 79 0.0151 0.0175 Native American 0.0163 0.00728
4 1.48e-272 0.0293 0.0326 Asian 0.0309 0.0697
5 3.51e- 58 0.189 0.197 Hispanic 0.193 0.163 All five hypothesis tests have p-values less than \(\alpha = 0.01\), thus they all reject \(H_0\), concluding that each Ethnicity was not sampled in perfect proportion with the 2010 US Census. In a survey like this, typically it is desirable to oversample smaller groups, such as Asian and Native American, and this was not done for Asian.
# Observed vs Expected proportions
library(ggplot2)
p <- ggplot(b_sum, aes(x = Ethnicity, y = Observed))
# observed
p <- p + geom_point(size = 4)
p <- p + geom_line(aes(group = 1))
p <- p + geom_errorbar(aes(ymin = CI.lower, ymax = CI.upper), width = 0.1)
# expected
p <- p + geom_point(aes(y = exp_prop), colour = "red", shape = "x", size = 6)
p <- p + geom_line(aes(y = exp_prop, group = 1), colour = "red")
p <- p + labs(title = "NESARC (w/99% CI) and 2010 US Census (red x) proportions")
p <- p + expand_limits(y = 0)
print(p)
5.7 Week11: Two-way categorical tables and Simple linear regression, inference
5.7.1 Tuesday ———
5.7.2 Two-way categorical analysis.
Using two categorical variables with two to five levels each, specify a hypothesis test for homogeneity of proportions associated with your research questions.
5.7.2.1 (1 p) Specify the hypotheses in words and notation.
* In words: "There is an association between TobaccoDependence and Ethnicity."
* In notation: $H_0: p(i \textrm{ and } j) = p(i)p(j)$ versus $H_A: p(i \textrm{ and } j) \ne p(i)p(j)$, for all row categories $i$ and column categories $j$.5.7.2.2 (1 p) State the conclusion of the test in the context of the problem.
# Tabulate by two categorical variables:
tab_E_TD <- xtabs(~ Ethnicity + TobaccoDependence, data = nesarc_sub)
tab_E_TD TobaccoDependence
Ethnicity Nicotine Dependence No Nicotine Dependence
Cauc 591 469
AfAm 104 107
NaAm 27 15
Asia 35 23
Hisp 139 196 TobaccoDependence
Ethnicity Nicotine Dependence No Nicotine Dependence
Cauc 0.65959821 0.57901235
AfAm 0.11607143 0.13209877
NaAm 0.03013393 0.01851852
Asia 0.03906250 0.02839506
Hisp 0.15513393 0.24197531
Pearson's Chi-squared test
data: tab_E_TD
X-squared = 25.423, df = 4, p-value = 4.135e-05The test statistic is \(X^2 = 25.4\).
Because the p-value \(= 4.14\times 10^{-5} < 0.05\) we reject the null hypothesis concluding that there is an association between TobaccoDependence and Ethnicity.
TobaccoDependence
Ethnicity Nicotine Dependence No Nicotine Dependence
Cauc 556.71747 503.28253
AfAm 110.81829 100.18171
NaAm 22.05862 19.94138
Asia 30.46190 27.53810
Hisp 175.94373 159.05627[1] 19.94138The model assumptions are met since the expected count for each cell is at least 5.
5.7.2.3 (1 p) Plot a mosaic plot of the data and Pearson residuals.
TobaccoDependence
Ethnicity Nicotine Dependence No Nicotine Dependence
Cauc 1.4529661 -1.5281534
AfAm -0.6476937 0.6812102
NaAm 1.0521058 -1.1065497
Asia 0.8222344 -0.8647830
Hisp -2.7851836 2.9293099# The sum of the squared residuals is the chi-squared statistic:
#chisq_E_TD$residuals^2
#sum(chisq_E_TD$residuals^2)# mosaic plot
library(vcd)
# this layout gives us the interpretation we want:
mosaic(tab_E_TD, shade=TRUE, legend=TRUE)
5.7.2.4 (1 p) Interpret the mosaic plot with reference to the Pearson residuals.
The primary cause of rejecting \(H_0\) is that Hispanics have less nicotine dependence than expected.
5.7.3 Thursday ———
5.7.4 Simple linear regression.
Select two numerical variables.
I don’t have great variables in my dataset related to my research question for this. I select TotalCigsSmokedLog10 and DaysSmoke.
5.7.4.1 (1 p) Plot the data and, if required, transform the variables so a roughly linear relationship is observed. All interpretations will be done on this scale of the variables.
See below.
5.7.4.2 (0 p) Fit the simple linear regression model.
(Intercept) DaysSmoke
0.4663599 0.6131081 Note that prediction intervals should be bounded below at 0.
library(ggplot2)
p <- ggplot(nesarc_sub, aes(x = DaysSmoke, y = TotalCigsSmokedSqrt))
p <- p + geom_vline(xintercept = 0, alpha = 0.25)
p <- p + geom_hline(yintercept = 0, alpha = 0.25)
# prediction bands
p <- p + geom_ribbon(aes(ymin = predict(lm_fit, data.frame(DaysSmoke)
, interval = "prediction", level = 0.95)[, 2],
ymax = predict(lm_fit, data.frame(DaysSmoke)
, interval = "prediction", level = 0.95)[, 3],)
, alpha=0.1, fill="darkgreen")
# linear regression fit and confidence bands
p <- p + geom_smooth(method = lm, se = TRUE)
# jitter a little to uncover duplicate points
p <- p + geom_jitter(position = position_jitter(.1), alpha = 0.75)
p <- p + labs(title = "Regression with confidence and prediction bands")
print(p)
5.7.4.3 (1 p) Assess the residuals for lack of fit (interpret plots of residuals vs fitted and \(x\)-value).
# plot diagnistics
par(mfrow=c(2,3))
plot(lm_fit, which = c(1,4,6))
# residuals vs DaysSmoke
# can only plot observations that were did not include NA for the lm_fit
plot(nesarc_sub %>% drop_na(TotalCigsSmokedSqrt, DaysSmoke) %>% pull(DaysSmoke), lm_fit$residuals, main="Residuals vs DaysSmoke")
# horizontal line at zero
abline(h = 0, col = "gray75")
# Normality of Residuals
library(car)
qqPlot(lm_fit$residuals, las = 1, id = list(n = 3), main="QQ Plot") 15 338 91
15 336 90 # # residuals vs order of data
# plot(lm_fit$residuals, main="Residuals vs Order of data")
# # horizontal line at zero
# abline(h = 0, col = "gray75")
hist(lm_fit$residuals, breaks=10, main="Residuals")
There is some nonrandom pattern. It almost appears that there’d be a shallower slope if it wasn’t for the 30 days/month smokers.
5.7.4.4 (1 p) Assess the residuals for normality (interpret QQ-plot and histogram).
Normality is not perfect, with some slight right skewness, but probably about as good as it gets for this dataset. The bigger problem is the mean pattern (mentioned above).
5.7.4.5 (1 p) Assess the relative influence of points.
There are a few points that have large Cook’s D, which are those few people who smoke 30 days/month and smoke the most and a person who smokes few days but smokes many times more. It may be worth seeing the effect of removing those four observations.
Seeing affect of removing 4 influential observations…
1704 1353 91 540 338 15
0.002887831 0.004125510 0.006439934 0.006439934 0.010526496 0.014499656 As seen below, removing the 4 largest Cook’s D observations (which is a rash move – usually we only remove one at a time since which points are influential can change based on removing a single observation) does not greatly improve the model fit or normality.
# data with deleted points with large Cook's D
nesarc_sub_D <- nesarc_sub %>% filter(!(rownames(nesarc_sub) %in% names(tail(sort(cooks.distance(lm_fit)), 4))))
dim(nesarc_sub)[1] 1706 25[1] 1702 25lm_fit_D <- lm(TotalCigsSmokedSqrt ~ DaysSmoke, data = nesarc_sub_D)
library(ggplot2)
p <- ggplot(nesarc_sub_D, aes(x = DaysSmoke, y = TotalCigsSmokedSqrt))
p <- p + geom_vline(xintercept = 0, alpha = 0.25)
p <- p + geom_hline(yintercept = 0, alpha = 0.25)
# prediction bands
p <- p + geom_ribbon(aes(ymin = predict(lm_fit_D, data.frame(DaysSmoke)
, interval = "prediction", level = 0.95)[, 2],
ymax = predict(lm_fit_D, data.frame(DaysSmoke)
, interval = "prediction", level = 0.95)[, 3],)
, alpha=0.1, fill="darkgreen")
# linear regression fit and confidence bands
p <- p + geom_smooth(method = lm, se = TRUE)
# jitter a little to uncover duplicate points
p <- p + geom_jitter(position = position_jitter(.1), alpha = 0.75)
p <- p + labs(title = "Regression with confidence and prediction bands")
print(p)
# plot diagnistics
par(mfrow=c(2,3))
plot(lm_fit_D, which = c(1,4,6))
# residuals vs DaysSmoke
# can only plot observations that were did not include NA for the lm_fit
plot(nesarc_sub_D %>% drop_na(TotalCigsSmokedSqrt, DaysSmoke) %>% pull(DaysSmoke), lm_fit_D$residuals, main="Residuals vs DaysSmoke")
# horizontal line at zero
abline(h = 0, col = "gray75")
# Normality of Residuals
library(car)
qqPlot(lm_fit_D$residuals, las = 1, id = list(n = 3), main="QQ Plot") 55 206 208
54 204 206 # # residuals vs order of data
# plot(lm_fit_D$residuals, main="Residuals vs Order of data")
# # horizontal line at zero
# abline(h = 0, col = "gray75")
hist(lm_fit_D$residuals, breaks=10, main="Residuals")
Continuing with the original data set …
5.7.4.6 (1 p) Test whether the slope is different from zero, \(H_A: \beta_1 \ne 0\).
Call:
lm(formula = TotalCigsSmokedSqrt ~ DaysSmoke, data = nesarc_sub)
Residuals:
Min 1Q Median 3Q Max
-13.382 -3.368 -0.688 5.635 35.362
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.46636 0.39314 1.186 0.236
DaysSmoke 0.61311 0.01456 42.112 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.637 on 1695 degrees of freedom
(9 observations deleted due to missingness)
Multiple R-squared: 0.5113, Adjusted R-squared: 0.511
F-statistic: 1773 on 1 and 1695 DF, p-value: < 2.2e-16First note that we have insufficient evidence to conclude that \(\beta_0\) is different from 0, which makes sense — if you smoke 0 days, then 0 should be the number of cigarettes you smoke.
We have sufficient evidence to conclude the slope is different from zero. Furthermore, for every additional day increase per month a person smokes, there is a predicted 0.6 square-root of cigarettes smoked. When this is squared to be on the original scale, the increase is scary.
5.7.4.7 (1 p) Interpret the \(R^2\) value.
The regression model explains \(R^2 = 0.511\) proportion of the variation in the response above using the mean value alone.
5.8 Week12: Logistic regression and Experiments vs Observational studies
5.8.1 Tuesday ———
5.8.2 Logistic Regression
If I provide the solutions here, then there’s nothing for you to do for the in-class assignment …
5.8.3 Thursday ———
5.8.4 Experiments and observational studies
- (2 p) Write a few sentences describing these characteristics of the data you’re using for your project.
- What is the sampling design?
- What is the target population and sampling frame?
- How was nonresponse handled? Answer this both in terms of individual nonresponse (when a respondant is lost to follow up, if this applies for you (e.g., AddHealth)) and question nonresponse (when a question is skipped).
You started doing this earlier in Week07 Tuesday, when you described the sampling.
6 Poster presentation
6.1 Week13: Complete poster in HW document
6.1.1 Tuesday ———
6.1.2 Thursday ———
6.1.2.1 Title
Smoking behavior is (barely) associated with major depression in young adults
6.1.2.2 1. (1 p) Introduction
Smoking behavior is associated with major depression.
Longitudinal investigations have shown that depression increases risk of later smoking (Breslau et al. 1998; Dierker et al. 2001).
Many subsequent studies find a role of major depression in increasing the probability and amount of smoking (Dierker et al. 2004; Rohde et al. 2004, 2003).
Smoking behavior and Nicotine dependence.
While smoking is a necessary requirement for nicotine dependence, smoking behavior is a poor predictor of developing nicotine dependence (Kandel and Chen 2000; Stanton, Lowe, and Silva 1995).
Heavy smokers may not have nicotine dependence, conversely light smokers may have nicotine dependence (Kandel and Chen 2000).
6.1.2.3 2. (1 p) Research Questions
- (Week02 Research Questions) 2 bullets, one for each research question, stated as the alternative hypothesis.
The goals of the analysis include answering these two questions:
Is there an association between major depression and nicotine dependence?
Does the prevalence of nicotine dependence differ by ethnicity?
6.1.2.4 3. (1 p) Methods
- (Week01 Personal Codebook) Data source(s).
A sample from the first wave of the National Epidemiologic Survey on Alcohol and Related Conditions (NESARC) (US adults). We focus on young adult (18–25) smokers (1706 of 43093 respondants).
In the Introduction to the National Epidemiologic Survey on Alcohol and Related Conditions they describe the sampling design. Adult respondants were sampled from a subset of households included in the US Census 2000 Supplementary Survey. To improve information about smaller groups, Blacks and Hispanics were oversampled by about 50%, as were young adults ages 18-24 by 125%. Using the US Census as a sampling frame helps to provide a representative sampling of the country. Oversampling smaller subsets of the population means that comparisons involving those subsets will be more reliable.
(Week01 Personal Codebook) Variables used.
Lifetime major depression (
Depression, 0=No, 1=Yes) (Grant et al. 2003, 1995).Total cigarettes per month (
TotalCigsSmoked= days smoked times daily cigs smoked, range: 0–NA), a combination of these questions: “About how often did you usually smoke in the past year?” and “On the days that you smoked in the last year, about how many cigarettes did you usually smoke?”Nicotine dependence (
TobaccoDependence, 0=No, 1=Yes) in the last 12 months.
(Various) Statistical methods used to answer the research questions.
A two-sample \(t\)-test comparing the mean total cigarettes smoked by depression status.
A \(\chi^2\) analysis of a two-way contingency table of nicotine dependence by ethnicity.
6.1.2.5 4. (1 p) Discussion (while this would follow your results, let’s put it here so you have a full column to show the results of the analysis of both research questions)
- Put the results you found for each research question in the context you provided in your introduction.
These results support the previous literature that individuals with major depression are more sensitive to the development of nicotine dependence regardless of how much they smoke. Thus, people with depression are an important population subgroup for targeted smoking intervention programs.
6.1.2.6 5. (1 p) Further directions or Future work or Next steps or something else that indicates there more to do and you’ve thought about it.
- What do these results lead you to want to investigate?
We can extend focus from young adults to all adults, and (with multiple regression techniques we’ll learn in ADA2) assess the particular depressed subpopulations most susceptible to developing nicotine dependence.
6.1.2.7 6. (1 p) References
By citing sources in your introduction, this section will automatically have your bibliography.
[Bibliography will go here – it’s currently at the bottom of the document]
6.1.2.8 7. (2 p) Results for your first research question.
- “Is the population mean square-root total cigarettes smoked different for those with depression or not?”
- \(H_0: \mu_{D} = \mu_{ND}\) versus \(H_A: \mu_{D} \ne \mu_{ND}\)
Welch Two Sample t-testdata: TotalCigsSmokedSqrt by Depression t = -1.7647, df = 763.27, p-value = 0.07802 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -1.67605335 0.08920479 sample estimates: mean in group No Depression mean in group Yes Depression 15.77958 16.57301
Welch Two Sample t-test
data: TotalCigsSmoked by Depression t = -1.8191, df = 733.22, p-value = 0.06931 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -59.186888 2.255475 sample estimates: mean in group No Depression mean in group Yes Depression 313.0782 341.5439
Model assumptions are met, the the sampling distribution of the difference in means is normal.
Because \(p=0.078 > 0.05\) (with \(t_{s} = -1.765\)), we have insufficient evidence to reject \(H_0\) at an \(\alpha=0.05\) significance level, concluding that the total cigarettes smoked does not differ by depression status.
The difference is very small, on the square-root scale it is 0.793 and on the natural scale it is 248.9952087, 274.6645376, a difference of 26 cigarettes each year.
6.1.2.9 8. (2 p) Results for your second research question.
- Is there an association between TobaccoDependence and Ethnicity?
The model assumptions are met since the expected count for each cell is at least 5.
Because the p-value \(= 4.14\times 10^{-5} < 0.05\) (\(X^2 = 25.4\)) we reject the null hypothesis concluding that there is an association between TobaccoDependence and Ethnicity.
The primary cause of rejecting \(H_0\) is that Hispanics have less nicotine dependence than expected by chance.
6.2 Week14: Posters, finishing up
6.2.1 Tuesday ———
6.2.2 Thursday ———
7 Additional items
7.1 Statistics summary tables
We won’t get to this in our class, but anyone in biostatistics or epidemiology will appreciate how easy it is to produce such an informative “Demographics” table.
library(moonBook)
out <-
moonBook::mytable(
Depression
+ TobaccoDependence
~
Age
+ SmokingFreq
+ DailyCigsSmoked
+ Ethnicity
+ Sex
+ DaysSmoke
+ TotalCigsSmoked
+ CigsSmokedFac
, data = nesarc_sub
)
myhtml(out)| Depression | No Depression | Yes Depression | ||||
|---|---|---|---|---|---|---|
| TobaccoDependence |
Nicotine Dependence (N=574) |
No Nicotine Dependence (N=686) |
p |
Nicotine Dependence (N=322) |
No Nicotine Dependence (N=124) |
p |
| Age | 21.5 ± 2.3 | 21.7 ± 2.1 | 0.145 | 21.5 ± 2.2 | 21.8 ± 2.2 | 0.205 |
| SmokingFreq | ||||||
| Never | 0 ( 0.0%) | 0 ( 0.0%) | 0 ( 0.0%) | 0 ( 0.0%) | ||
| Once a month or less | 6 ( 1.0%) | 43 ( 6.3%) | 1 ( 0.3%) | 21 (16.9%) | ||
| 2 to 3 Days/month | 8 ( 1.4%) | 41 ( 6.0%) | 4 ( 1.2%) | 12 ( 9.7%) | ||
| 1 to 2 Days/week | 12 ( 2.1%) | 60 ( 8.8%) | 7 ( 2.2%) | 9 ( 7.3%) | ||
| 3 to 4 Days/week | 22 ( 3.8%) | 50 ( 7.3%) | 10 ( 3.1%) | 9 ( 7.3%) | ||
| 5 to 6 Days/week | 16 ( 2.8%) | 34 ( 5.0%) | 11 ( 3.4%) | 7 ( 5.6%) | ||
| Every Day | 510 (88.9%) | 455 (66.6%) | 289 (89.8%) | 66 (53.2%) | ||
| DailyCigsSmoked | 13.6 ± 9.0 | 9.0 ± 7.5 | 0.000 | 13.8 ± 9.5 | 7.3 ± 6.1 | 0.000 |
| Ethnicity | 0.000 | 0.053 | ||||
| Cauc | 377 (65.7%) | 402 (58.6%) | 214 (66.5%) | 67 (54.0%) | ||
| AfAm | 77 (13.4%) | 86 (12.5%) | 27 ( 8.4%) | 21 (16.9%) | ||
| NaAm | 12 ( 2.1%) | 9 ( 1.3%) | 15 ( 4.7%) | 6 ( 4.8%) | ||
| Asia | 25 ( 4.4%) | 20 ( 2.9%) | 10 ( 3.1%) | 3 ( 2.4%) | ||
| Hisp | 83 (14.5%) | 169 (24.6%) | 56 (17.4%) | 27 (21.8%) | ||
| Sex | 0.633 | 0.031 | ||||
| Female | 263 (45.8%) | 304 (44.3%) | 194 (60.2%) | 89 (71.8%) | ||
| Male | 311 (54.2%) | 382 (55.7%) | 128 (39.8%) | 35 (28.2%) | ||
| DaysSmoke | 28.0 ± 6.2 | 22.8 ± 10.9 | 0.000 | 28.3 ± 5.6 | 19.1 ± 12.7 | 0.000 |
| TotalCigsSmoked | 395.5 ± 279.7 | 243.4 ± 238.0 | 0.000 | 404.4 ± 293.6 | 177.6 ± 195.5 | 0.000 |
| CigsSmokedFac | ||||||
| Q1 | 0 ( 0.0%) | 0 ( 0.0%) | 0 ( 0.0%) | 0 ( 0.0%) | ||
| Q2 | 0 ( 0.0%) | 0 ( 0.0%) | 0 ( 0.0%) | 0 ( 0.0%) | ||
| Q3 | 310 (54.0%) | 520 (76.6%) | 172 (53.6%) | 105 (85.4%) | ||
| Q4 | 264 (46.0%) | 159 (23.4%) | 149 (46.4%) | 18 (14.6%) | ||
References (from Week02)
Beckschäfer, P, L Fehrmann, Rd Harrison, J Xu, and C Kleinn. 2014. “Mapping Leaf Area Index in Subtropical Upland Ecosystems Using RapidEye Imagery and the randomForest Algorithm.” iForest - Biogeosciences and Forestry 7 (1): 1–11. https://doi.org/10.3832/ifor0968-006.
Breslau, Naomi. 1995. “Psychiatric Comorbidity of Smoking and Nicotine Dependence.” Behavior Genetics 25 (2): 95–101.
Breslau, N., E. L. Peterson, L. R. Schultz, H. D. Chilcoat, and P. Andreski. 1998. “Major Depression and Stages of Smoking. A Longitudinal Investigation.” Archives of General Psychiatry 55 (2): 161–66.
Caraballo, Ralph S., Scott P. Novak, and Katherine Asman. 2009. “Linking Quantity and Frequency Profiles of Cigarette Smoking to the Presence of Nicotine Dependence Symptoms Among Adolescent Smokers: Findings from the 2004 National Youth Tobacco Survey.” Nicotine & Tobacco Research, January, ntn008. https://doi.org/10.1093/ntr/ntn008.
Dean, Jared, and ebrary, Inc. 2014. Big Data, Data Mining, and Machine Learning Value Creation for Business Leaders and Practitioners. Wiley & SAS Business Series. Hoboken, NJ: Wiley.
Dierker, Lisa C., Shelli Avenevoli, Abbie Goldberg, and Meyer Glantz. 2004. “Defining Subgroups of Adolescents at Risk for Experimental and Regular Smoking.” Prevention Science: The Official Journal of the Society for Prevention Research 5 (3): 169–83.
Dierker, Lisa C., Shelli Avenevoli, Kathleen R. Merikangas, Brian P. Flaherty, and Marilyn Stolar. 2001. “Association Between Psychiatric Disorders and the Progression of Tobacco Use Behaviors.” Journal of the American Academy of Child & Adolescent Psychiatry 40 (10): 1159–67. https://doi.org/10.1097/00004583-200110000-00009.
Dierker, Lisa C., Eric Donny, Stephen Tiffany, Suzanne M. Colby, Nicholas Perrine, and Richard R. Clayton. 2007. “The Association Between Cigarette Smoking and DSM-IV Nicotine Dependence Among First Year College Students.” Drug and Alcohol Dependence 86 (2–3): 106–14. https://doi.org/10.1016/j.drugalcdep.2006.05.025.
Grant, B. F., T. C. Harford, D. A. Dawson, P. S. Chou, and R. P. Pickering. 1995. “The Alcohol Use Disorder and Associated Disabilities Interview Schedule (AUDADIS): Reliability of Alcohol and Drug Modules in a General Population Sample.” Drug and Alcohol Dependence 39 (1): 37–44.
Grant, Bridget F., Deborah A. Dawson, Frederick S. Stinson, Patricia S. Chou, Ward Kay, and Roger Pickering. 2003. “The Alcohol Use Disorder and Associated Disabilities Interview Schedule-IV (AUDADIS-IV): Reliability of Alcohol Consumption, Tobacco Use, Family History of Depression and Psychiatric Diagnostic Modules in a General Population Sample.” Drug and Alcohol Dependence 71 (1): 7–16.
Kandel, D. B., and K. Chen. 2000. “Extent of Smoking and Nicotine Dependence in the United States: 1991-1993.” Nicotine & Tobacco Research: Official Journal of the Society for Research on Nicotine and Tobacco 2 (3): 263–74.
Murphy, Kevin P. 2012. Machine Learning a Probabilistic Perspective. Adaptive Computation and Machine Learning Series. Cambridge, Mass: MIT Press.
Richert, Willi. 2013. Building Machine Learning Systems with Python. Birmingham, UK: Packt Publishing.
Rohde, Paul, Christopher W. Kahler, Peter M. Lewinsohn, and Richard A. Brown. 2004. “Psychiatric Disorders, Familial Factors, and Cigarette Smoking: II. Associations with Progression to Daily Smoking.” Nicotine & Tobacco Research 6 (1): 119–32. https://doi.org/10.1080/14622200310001656948.
Rohde, Paul, Peter M. Lewinsohn, Richard A. Brown, Jeffrey M. Gau, and Christopher W. Kahler. 2003. “Psychiatric Disorders, Familial Factors and Cigarette Smoking: I. Associations with Smoking Initiation.” Nicotine & Tobacco Research 5 (1): 85–98. https://doi.org/10.1080/1462220031000070507.
Stanton, Warren R., John B. Lowe, and Phil A. Silva. 1995. “Antecedents of Vulnerability and Resilience to Smoking Among Adolescents.” Journal of Adolescent Health 16 (1): 71–77. https://doi.org/10.1016/1054-139X(94)00051-F.