Week 8 Case-Control RQ
R version 3.4.4 (2018-03-15) -- "Someone to Lean On"
> install.packages("epiDisplay") ## successor to Epicalc package (see epicalc book)
> data(VC1to1) # matched data in long-form, as would be produced from pair-matching with Matchit
# in case control studies the pair matching (on covariates) uses the outcome (cancer or not) not the treatment as in our prior matching analyses
# then associations with possible causes (exposures) are sought)-- here smoking alcohol rubber
> head(VC1to1)
matset case smoking rubber alcohol
1 1 1 1 0 0
2 1 0 1 0 0
3 2 1 1 0 1
4 2 0 1 1 0
5 3 1 1 1 0
6 3 0 1 1 0
# the epicalc book uses a builtin function matchTab to get odds ratios-- we can also use glm or glmer
> with(VC1to1, matchTab(case, alcohol, matset))
Exposure status: alcohol = 1
Total number of match sets in the tabulation = 26
Number of controls = 1
No. of controls exposed
No. of cases exposed 0 1
0 7 2
1 9 8
Odds ratio by Mantel-Haenszel method = 4.5
Odds ratio by maximum likelihood estimate (MLE) method = 4.5
95%CI= 0.972 , 20.827
> with(VC1to1, matchTab(case, smoking, matset))
Exposure status: smoking = 1
Total number of match sets in the tabulation = 26
Number of controls = 1
No. of controls exposed
No. of cases exposed 0 1
0 0 5
1 5 16
Odds ratio by Mantel-Haenszel method = 1
Odds ratio by maximum likelihood estimate (MLE) method = 1
95%CI= 0.29 , 3.454
> with(VC1to1, matchTab(case, rubber, matset))
Exposure status: rubber = 1
Total number of match sets in the tabulation = 26
Number of controls = 1
No. of controls exposed
No. of cases exposed 0 1
0 13 5
1 4 4
Odds ratio by Mantel-Haenszel method = 0.8
Odds ratio by maximum likelihood estimate (MLE) method = 0.8
95%CI= 0.215 , 2.979
# so alcohol appears to be the most potent cause, not quite significant
## try some glm or glmer, look for association (odds ratio away from 1)
> or1a = glm(case~ smoking , data = VC1to1, family = binomial) # smoking not important
> summary(or1a)
Call:
glm(formula = case ~ smoking, family = binomial, data = VC1to1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.177 -1.177 0.000 1.177 1.177
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.262e-15 6.325e-01 0 1
smoking 1.563e-15 7.037e-01 0 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 72.087 on 51 degrees of freedom
Residual deviance: 72.087 on 50 degrees of freedom
AIC: 76.087
Number of Fisher Scoring iterations: 2
> or1b = glm(case~ rubber , data = VC1to1, family = binomial)
> summary(or1b)
Call:
glm(formula = case ~ rubber, family = binomial, data = VC1to1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.20178 -1.20178 0.01271 1.15324 1.22782
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.05716 0.33820 0.169 0.866
rubber -0.17494 0.59202 -0.295 0.768
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 72.087 on 51 degrees of freedom
Residual deviance: 72.000 on 50 degrees of freedom
AIC: 76
Number of Fisher Scoring iterations: 3
> or1c = glm(case~ alcohol , data = VC1to1, family = binomial) # alcohol nearly significant
> summary(or1c)
Call:
glm(formula = case ~ alcohol, family = binomial, data = VC1to1)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.40943 -0.94476 0.00857 0.96190 1.42944
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5754 0.4167 -1.381 0.1673
alcohol 1.1060 0.5766 1.918 0.0551 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 72.087 on 51 degrees of freedom
Residual deviance: 68.265 on 50 degrees of freedom
AIC: 72.265
Number of Fisher Scoring iterations: 4
> exp(coef(or1c))
(Intercept) alcohol
0.562500 3.022222
> exp(confint(or1c))
Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 0.2380299 1.247758
alcohol 0.9972608 9.708124
> cbind(exp(coef(or1c)), exp(confint(or1c))) #here's a nice thing students did on exams, presume some class did this, good for them
Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 0.562500 0.2380299 1.247758
alcohol 3.022222 0.9972608 9.708124
> library(lme4)
Loading required package: Matrix
> or1 = glmer(case~ smoking + (smoking|matset), data = VC1to1, family = binomial) #smoking bombs, nothing there
Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GHrule(0L), compDev = compDev, :
Downdated VtV is not positive definite
> or3 = glmer(case~ alcohol + (alcohol|matset), data = VC1to1, family = binomial) # get the same as glm here, but a little diff than Mantel-Haensel
> summary(or3)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial ( logit )
Formula: case ~ alcohol + (alcohol | matset)
Data: VC1to1
AIC BIC logLik deviance df.resid
78.3 88.0 -34.1 68.3 47
Scaled residuals:
Min 1Q Median 3Q Max
-1.30384 -0.75000 0.00848 0.76696 1.33333
Random effects:
Groups Name Variance Std.Dev. Corr
matset (Intercept) 5.413e-15 7.357e-08
alcohol 1.648e-14 1.284e-07 -1.00
Number of obs: 52, groups: matset, 26
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.5754 0.4167 -1.381 0.1673
alcohol 1.1060 0.5766 1.918 0.0551 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr)
alcohol -0.723
> exp(fixef(or3))
(Intercept) alcohol
0.562500 3.022222
>