Social Survey Example Revisit

[Data] In the 2000 General Social Survey, 1144 subjects were asked whether, to help the environment, they would be willing to (1) pay higher taxes or (2) accept a cut in living standards.

Social Survey Example
		Cut Living Standards
Pay Higher Taxes	Yes	No	Total
Yes	227	132	359
No	107	678	785
Total	334	810	1144

[Goal] How can we compare the probabilities of a ``yes" outcome for the two environmental questions?

[Code]

Note that specifies the number of quadrature points to be 10

library(lme4)   # install packages by install.packages("lme4")

## Loading required package: Matrix

survey = NULL
for (subject in 1:227){
    survey = rbind(survey, c(1, 1, subject))
    survey = rbind(survey, c(1, 2, subject))
}
for (subject in (227+1):(227+132)){
    survey = rbind(survey, c(1, 1, subject))
    survey = rbind(survey, c(0, 2, subject))
}
for (subject in (227+132+1):(227+132+107)){
    survey = rbind(survey, c(0, 1, subject))
    survey = rbind(survey, c(1, 2, subject))
}
for (subject in (227+132+107+1):(227+132+107+678)){
    survey = rbind(survey, c(0, 1, subject))
    survey = rbind(survey, c(0, 2, subject))
}
survey = as.data.frame(survey)
names(survey) = c("response", "question", "subject")

R.fit = glmer(response ~ question+(1|subject), 
              data=survey, family=binomial, 
              control = glmerControl(optimizer = "bobyqa"), nAGQ = 10)
summary(R.fit)

## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
##  Family: binomial  ( logit )
## Formula: response ~ question + (1 | subject)
##    Data: survey
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   2527.4   2544.6  -1260.7   2521.4     2285 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.8866 -0.2708 -0.2438  0.4666  1.2528 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  subject (Intercept) 8.035    2.835   
## Number of obs: 2288, groups:  subject, 1144
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.4027     0.2350  -5.969 2.39e-09 ***
## question     -0.2099     0.1300  -1.614    0.107    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## question -0.794

Basketball player example

library(lme4)

nn = c(13, 10, 15, 14, 6, 10, 10, 4, 11, 10, 8, 9, 9, 8, 6)
pp = c(0.769, 0.9, 0.667, 0.643, 0.667, 0.9, 0.6, 1, 0.545, 0.9, 0.5, 0.889, 0.778, 0.625, 0.167)
success = round(nn*pp, 0)
failure = nn - success
mydata = data.frame(subject = c(rep(1:15, success), rep(1:15, failure)), 
                    score = c(rep(1, sum(success)), rep(0, sum(failure)))
                   )

mydata[1:15,]

##    subject score
## 1        1     1
## 2        1     1
## 3        1     1
## 4        1     1
## 5        1     1
## 6        1     1
## 7        1     1
## 8        1     1
## 9        1     1
## 10       1     1
## 11       2     1
## 12       2     1
## 13       2     1
## 14       2     1
## 15       2     1

R.fit = glmer(score ~ 1 + (1|subject), 
              data=mydata, family=binomial, 
              control = glmerControl(optimizer = "bobyqa"), nAGQ = 10)
summary(R.fit)

## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
##  Family: binomial  ( logit )
## Formula: score ~ 1 + (1 | subject)
##    Data: mydata
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##    176.7    182.7    -86.4    172.7      141 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7821 -1.4021  0.5704  0.6607  0.8037 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  subject (Intercept) 0.1779   0.4217  
## Number of obs: 143, groups:  subject, 15
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   0.9076     0.2244   4.044 5.24e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R.fit0 = glm(score ~ 1, data=mydata, family=binomial)

Reading the output,

\(\hat\alpha\) = 0.9076
\(\hat\sigma^2\) = 0.1779
Randome effects \(u_i\) (as in \(\alpha+u_i\)) can be viewed by ranef(R.fit).

ranef(R.fit)

## $subject
##    (Intercept)
## 1   0.08956086
## 2   0.24804972
## 3  -0.07862137
## 4  -0.11393761
## 5  -0.04008209
## 6   0.24804972
## 7  -0.14550175
## 8   0.17941600
## 9  -0.23025489
## 10  0.24804972
## 11 -0.23168174
## 12  0.21513597
## 13  0.07901647
## 14 -0.09597728
## 15 -0.47051178

The random effect can be tested by

2*(logLik(R.fit) - logLik(R.fit0))

## 'log Lik.' 0.4234963 (df=2)

0.5*(1-pchisq(0.42,df=1))

## [1] 0.2584685

Treatment for Depression Example Revisit

[Data] Longitudinal Study of Treatment for Depression

Subjects were classified into two groups based on their diagnosis severity (mild/severe)
In each group, subjects were randomly assigned to one of two drugs
Dichotomous assessment of each subject’s extent of suffering from mental depression was made at weeks 1, 2, and 4

depression = read.csv("../data/depression.csv",header=TRUE)
head(depression)

##       drug severity response time subject
## 1 Standard     Mild        N    0       1
## 2 Standard     Mild        N    1       1
## 3 Standard     Mild        N    2       1
## 4 Standard     Mild        N    0       2
## 5 Standard     Mild        N    1       2
## 6 Standard     Mild        N    2       2

depression$drug <- relevel(depression$drug , ref = "Standard")

R.fit = glmer((response=="N") ~ 1 + severity + drug*time + (1|subject),
              data=depression, family=binomial,
              control = glmerControl(optimizer = "bobyqa"), nAGQ = 10)

summary(R.fit)

## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod]
##  Family: binomial  ( logit )
## Formula: (response == "N") ~ 1 + severity + drug * time + (1 | subject)
##    Data: depression
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   1173.9   1203.5   -581.0   1161.9     1014 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2835 -0.8264  0.2324  0.7963  2.0191 
## 
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  subject (Intercept) 0.004331 0.06581 
## Number of obs: 1020, groups:  subject, 340
## 
## Fixed effects:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -0.02795    0.16412  -0.170    0.865    
## severitySevere    -1.31521    0.15465  -8.505  < 2e-16 ***
## drugNew Drug      -0.05970    0.22246  -0.268    0.788    
## time               0.48284    0.11596   4.164 3.13e-05 ***
## drugNew Drug:time  1.01842    0.19241   5.293 1.20e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) svrtyS drgNwD time  
## severitySvr -0.385                     
## drugNewDrug -0.614 -0.004              
## time        -0.671 -0.133  0.522       
## drgNwDrg:tm  0.461 -0.134 -0.740 -0.552

Parameter estimates are rather close to their GEE counterparts.

Random effects: Generalized Linear Mixed Models

Social Survey Example Revisit

Basketball player example

Treatment for Depression Example Revisit