Example for GLMM: modeling correlated survey response

The data is a survey result for whether the respondents support legalizing abortion under three situations: 1) if the family has a very low income and cannot afford anymore children, 2) when the women is not married and does not want to marry the man, 3) when the women wants it for any reason. We aim to understand how one person’s opinion change under different situations and whether gender has an effect.

Abortion <- read.table("Abortion.dat",header=T)
Abortion

1.1 Applying the GLMM model

We want to use a GLMM model to take into account the correlation of the choices from the same individual.

library(glmmML) ## do not install package from source if you face an installation error
## one can also use the glmer function in the lme4 package
fit.glmm <- glmmML(response ~ gender + factor(situation),
                   cluster = case, 
                   data = Abortion,
                   family = binomial, 
                   method = "ghq", n.points = 70)
summary(fit.glmm)
## 
## Call:  glmmML(formula = response ~ gender + factor(situation), family = binomial,      data = Abortion, cluster = case, method = "ghq", n.points = 70) 
## 
## 
##                        coef se(coef)        z Pr(>|z|)
## (Intercept)         0.21596   0.3766  0.57347 5.66e-01
## gender              0.01257   0.4888  0.02572 9.79e-01
## factor(situation)2 -0.54230   0.1572 -3.45022 5.60e-04
## factor(situation)3 -0.83470   0.1601 -5.21346 1.85e-07
## 
## Scale parameter in mixing distribution:  8.736 gaussian 
## Std. Error:                              0.5421 
## 
##         LR p-value for H_0: sigma = 0:  0 
## 
## Residual deviance: 4579 on 5545 degrees of freedom   AIC: 4589

Intepretations of the results: - Intepretation of coefficients - individual heteogeneity: our estimate random effect is \(\hat \sigma_u = 8.736\), indicating strong correlation among the choices within one individual. - results change with n.points

Using the glmer function from lme4 will give you the same solution

library(lme4)
## Loading required package: Matrix
fit.glmm2 <- glmer(response ~ gender + factor(situation) + (1|case), data = Abortion,
                   family = binomial, nAGQ = 70)
summary(fit.glmm2)
## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 70) [glmerMod]
##  Family: binomial  ( logit )
## Formula: response ~ gender + factor(situation) + (1 | case)
##    Data: Abortion
## 
##      AIC      BIC   logLik deviance df.resid 
##   4588.5   4621.7  -2289.3   4578.5     5545 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7809 -0.1224 -0.1056  0.1397  1.7149 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  case   (Intercept) 76.32    8.736   
## Number of obs: 5550, groups:  case, 1850
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.21596    0.37684   0.573 0.566586    
## gender              0.01258    0.48922   0.026 0.979480    
## factor(situation)2 -0.54230    0.15720  -3.450 0.000561 ***
## factor(situation)3 -0.83470    0.16009  -5.214 1.85e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) gender fct()2
## gender      -0.727              
## fctr(sttn)2 -0.206  0.000       
## fctr(sttn)3 -0.206  0.000  0.512

1.2 Compare with the results from the standard logistic regression

Let’s look at the marginal model results. Here we just use standard GLM ignoring correlation structure. The book uses a more sophisticated marginal model which provides valid inference. Our standard glm here would not provide valid inference for the marginal model, but the estimates should still be good.

(Why? Brief explanation: though the correlation structure of the samples is not considered, the score equations from the standard GLM is typically still consistent. Specifically, the expectation of score equation at the true parameter value is 0, and when correlation across samples are weak, the consistency of score equation can be obtained from the LLN)

fit.glm <- glm(response ~ gender + factor(situation), data = Abortion, family = binomial)
summary(fit.glm)
## 
## Call:
## glm(formula = response ~ gender + factor(situation), family = binomial, 
##     data = Abortion)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.189  -1.148  -1.125   1.207   1.231  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)         0.023940   0.055528   0.431   0.6664  
## gender              0.003582   0.054138   0.066   0.9472  
## factor(situation)2 -0.097329   0.065783  -1.480   0.1390  
## factor(situation)3 -0.149347   0.065825  -2.269   0.0233 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7689.5  on 5549  degrees of freedom
## Residual deviance: 7684.2  on 5546  degrees of freedom
## AIC: 7692.2
## 
## Number of Fisher Scoring iterations: 3

The marginal model coefficients are very scaled down towards \(0\) compared with GLMM. We can check that the ratio is roughly \(1/\sqrt{1 + \sigma_u^2/1.7^2}\).

plot(fit.glm$coefficients, fit.glmm$coefficients, xlim = c(-0.5, 0.25), ylim = c(-1, 0.25))
abline(v = 0, lty = 2)
abline(h = 0, lty = 2)
abline(a = 0, b = sqrt(1 + (fit.glmm$sigma/1.7)^2), col = "red")