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
```

- case: A respondent individual
- situation: one of the three situations
- response: yes (1) or no (0)
- gender: the gender of the respondent

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
```

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")
```