# ~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~ CRP 241 Module 4 Day 8 ~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~
# Example 0
# Download and load the linearity dataset:
download.file("http://www.duke.edu/~sgrambow/crp241data/linearity.RData",
destfile = "linearity.RData",quiet=TRUE,mode="wb",cacheOK=FALSE)
load("linearity.RData")
# Do the following:
# (1) plot the data
# (2) run a simple linear regression model
# (3) check the diagnostics
# Example 1
# FEV1 and Genetic Variation:
# A study was conducted to determine if there was an association between a genetic
# variant and forced expiratory volume in one second (FEV1) in patients with COPD.
# FEV1 was measured in liters. Genotypes of interest were wild type vs. mutant
# (i.e. heterozygous or homozygous for risk allele). Sex was also recorded.
# 100 patients were randomly selected from the researchers clinical practice.
# Data Dictionary:
# 1. FEV1 forced expiratory volume in one second (liters)
# 2. GENO patient genotyps (1 = Mutant; 0 = Wild Type)
# 3. SEX patient sex (1 = Female; 0 = Male)
# Download and load the lead dataset used in lecture:
download.file("http://www.duke.edu/~sgrambow/crp241data/fev1_geno.RData",
destfile = "fev1_geno.RData",quiet=TRUE,mode="wb",cacheOK=FALSE)
load("fev1_geno.RData")
# (1) Two Sample T-test Analysis
# - Assuming popluation variances are equal
t.test(fgdata$FEV1~fgdata$GENO,var.equal=T)
##
## Two Sample t-test
##
## data: fgdata$FEV1 by fgdata$GENO
## t = 2.7962, df = 198, p-value = 0.005681
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.1228584 0.7108077
## sample estimates:
## mean in group 0 mean in group 1
## 3.641340 3.224507
3.641340 -3.224507
## [1] 0.416833
# Is sex a confounder of the relationship between FEV1 and genotype?
# - Add labels to sex and genotype variables to make ouptut easier to interpret
fgdata$fSEX <- factor(fgdata$SEX,labels=c('Male','Female'))
fgdata$fGENO <- factor(fgdata$GENO,labels=c('Wild Type','Mutant'))
# - Check distribution of sex by each genotype
cross.tab <- table(fgdata$fSEX,fgdata$fGENO) # Counts
cross.tab
##
## Wild Type Mutant
## Male 70 30
## Female 42 58
prop.table(cross.tab) # Proportions among genotype
##
## Wild Type Mutant
## Male 0.35 0.15
## Female 0.21 0.29
# - Check distribution of FEV1 by each genotype
boxplot(fgdata$FEV1~fgdata$fSEX,
ylab='FEV1 Level (liters)')
# (2) Unadjusted Linear Regression Analysis
ufit <- lm(FEV1 ~ GENO, data=fgdata)
summary(ufit)
##
## Call:
## lm(formula = FEV1 ~ GENO, data = fgdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.35865 -0.68975 -0.02921 0.59003 2.63869
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.64134 0.09888 36.824 < 2e-16 ***
## GENO -0.41683 0.14907 -2.796 0.00568 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.046 on 198 degrees of freedom
## Multiple R-squared: 0.03799, Adjusted R-squared: 0.03313
## F-statistic: 7.819 on 1 and 198 DF, p-value: 0.005681
confint(ufit)
## 2.5 % 97.5 %
## (Intercept) 3.4463389 3.8363403
## GENO -0.7108077 -0.1228584
# (3) Adjusted Linear Regression Analysis
afit <- lm(FEV1 ~ GENO + SEX, data=fgdata)
summary(afit)
##
## Call:
## lm(formula = FEV1 ~ GENO + SEX, data = fgdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.23691 -0.69417 -0.01816 0.77985 2.36431
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9157 0.1081 36.213 < 2e-16 ***
## GENO -0.2090 0.1467 -1.425 0.156
## SEX -0.7317 0.1456 -5.025 1.12e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9877 on 197 degrees of freedom
## Multiple R-squared: 0.1473, Adjusted R-squared: 0.1386
## F-statistic: 17.02 on 2 and 197 DF, p-value: 1.526e-07
confint(afit)
## 2.5 % 97.5 %
## (Intercept) 3.7024806 4.12896143
## GENO -0.4981893 0.08025266
## SEX -1.0188149 -0.44455282
# (4) Stratified Analysis
females <- subset(fgdata,SEX==1)
males <- subset(fgdata,SEX==0)
# - How does "adjusting" work
by(females$FEV1,females$GENO,summary) # Mean difference among females
## females$GENO: 0
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.283 2.420 3.152 3.141 3.982 4.947
## --------------------------------------------------------
## females$GENO: 1
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.335 2.648 3.030 3.006 3.644 4.716
3.006 - 3.141
## [1] -0.135
by(males$FEV1,males$GENO,summary) # Mean difference among males
## males$GENO: 0
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.679 3.170 3.898 3.942 4.838 6.280
## --------------------------------------------------------
## males$GENO: 1
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.589 3.011 3.431 3.646 4.545 5.673
3.646 - 3.942
## [1] -0.296
(3.006 - 3.141 + 3.646 - 3.942)/2 # Average of strata-differences
## [1] -0.2155
# - (4a) Unadjusted Linear Regression Analysis Among Females
ffit <- lm(FEV1 ~ GENO, data=females)
summary(ffit)
##
## Call:
## lm(formula = FEV1 ~ GENO, data = females)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85813 -0.63437 0.02404 0.72909 1.80653
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1408 0.1354 23.202 <2e-16 ***
## GENO -0.1345 0.1777 -0.756 0.451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8773 on 98 degrees of freedom
## Multiple R-squared: 0.005805, Adjusted R-squared: -0.004339
## F-statistic: 0.5723 on 1 and 98 DF, p-value: 0.4512
confint(ffit)
## 2.5 % 97.5 %
## (Intercept) 2.8721889 3.4094579
## GENO -0.4871961 0.2182726
# - (4a) Unadjusted Linear Regression Analysis Among Males
mfit <- lm(FEV1 ~ GENO, data=males)
summary(mfit)
##
## Call:
## lm(formula = FEV1 ~ GENO, data = males)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.26284 -0.70681 -0.09617 0.91070 2.33838
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.9416 0.1303 30.249 <2e-16 ***
## GENO -0.2954 0.2379 -1.242 0.217
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.09 on 98 degrees of freedom
## Multiple R-squared: 0.01549, Adjusted R-squared: 0.005442
## F-statistic: 1.542 on 1 and 98 DF, p-value: 0.2173
confint(mfit)
## 2.5 % 97.5 %
## (Intercept) 3.6830593 4.2002394
## GENO -0.7675146 0.1767227
# Example 2
# Lead Study Data:
# A study was performed of the effects of exposure to lead on the psychological and
# neurological well-being of children. The data for this study are provided in the
# lead Dataset, posted on the course web site. In summary, a group of children living
# near a lead smelter in El Paso, Texas, were identified and their blood levels of
# lead were measured.
# An exposed group of 36 children were identified who had blood-lead levels >= 40 ug/ml.
# This group is defined by the variable Group = 2. A control group of 66 children
# were also identified who had blood-lead levels <40 ug/ml. This group is identified
# by the variable GROUP = 1. All children lived in close proximity to the lead smelter.
# One of the key outcome variables studied was the number of finger-wrist taps in the
# dominant hand (#taps in one ten second trial), a measure of neurological function (MAXFWT).
# Data Dictionary:
# There are numerous variables in the lead dataset but those of interest for the
# current exercise include:
# 1. ageyrs age of child in years xx.xx
# 2. Group exposure group (1= Control, 2= Exposed)
# 3. maxfwt finger-wrist tapping test in dominant hand (max of right and left hands)
# Download and load the lead dataset used in lecture:
download.file("http://www.duke.edu/~sgrambow/crp241data/lead.RData",
destfile = "lead.RData",quiet=TRUE,mode="wb",cacheOK=FALSE)
load("lead.RData")
# Note: Missing values in maxfwt are coded as a 99.
# - Recode as an NA so that R treats them correctly as missing values.
lead$maxfwt[lead$maxfwt==99] <- NA
# Question 1:
# Is age a confounder of the relationship between lead exposure and the
# score of the finger-wrist tapping test?
# Question 2:
# Estimate the unadjusted association between lead exposure and the
# score of the finger-wrist tapping test. Is it statistically significant?
# Question 3:
# Estimate the association between lead exposure and the score of the
# finger-wrist tapping test adjusted for age. Is it statistically significant?
# Question 4:
# What is the impact of adjusting for age on the conclusions of the anlysis?
# Question 5:
# How were the missing values handled in the regression analysis?
# End of Program