GLMMs and GEEs with Binary Data

Data on attitudes about abortion from Agresti (2013)

Subjects were asked whether they supported legalizing abortion in each of three situations:

1. If the family has a very low income and cannot afford any more children.
2. When the woman is not married and does not want to marry the man.
3. When the woman wants it for any reason.

Gender was also measured.

Note: Since the model implies nonnegative association among responses, situations and scales should correspond accordingly. For example, with scale (yes, no), it would not make sense to ask “Should abortion be legal…?” and ” Should abortion be illegal…?” as two different situations.

Read in data

Variables:

* gender = 1 if female; 0 if male
* response = 1 if support; 0 if do not support
* question = situation (1, 2, or 3)
* case = subject id

attitudes <- read.table("https://math.montana.edu/shancock/data/abortion.txt", header=TRUE)
# Code situation as factor:
attitudes$question = factor(attitudes$question)
# Make situation 3 the reference situation:
attitudes$question = relevel(attitudes$question, ref="3")
# Note that we can leave gender and response as numeric
# since they are already coded as 0 or 1.

Exploratory data analysis

Plot of proportion of support in each covariate pattern:

means <- tapply(attitudes$response, list(attitudes$question,attitudes$gender), mean)
plot(c(3,1,2),means[1:3,1],type="p",col="red",pch=16,
    xlab="Situation",ylab="Proportion Support",
    main="Proportion Support by Situation and Gender",
    cex=2,ylim=c(0.45,.52))
points(c(3,1,2),means[1:3,2],col="darkblue",cex=2)  
legend("bottomleft", c("Male","Female"), col=c("red","darkblue"), pch=c(16,1))

Barplot of same information:

barplot(t(means), beside=TRUE, col=c("Red","Blue"),
    ylim=c(0,0.6), legend = c("Male","Female"),
    xlab="Situation",ylab="Proportion Support",
    main="Proportion Support by Situation and Gender")

# Convert data to long form - currently first col is male, second is female
means_long <- data.frame(Proportion = c(means[,1], means[,2]))
# Add gender and situation variables
means_long$Gender <- rep(c("Male", "Female"), each = 3)
means_long$Situation <- rep(c("Any reason", "Income", "Unmarried"), 2)

means_long %>% ggplot(aes(x = Situation, fill = Gender, y = Proportion)) +
  geom_col(position = "dodge")

barplot(t(means), beside=TRUE, col=c("Red","Blue"),
    ylim=c(0,0.6), legend = c("Male","Female"),
    xlab="Situation",ylab="Proportion Support",
    main="Proportion Support by Situation and Gender")

Contingency table summary of data by situation:

my.tab <- xtabs(~response+gender+question, data=attitudes)
my.tab

## , , question = 3
## 
##         gender
## response   0   1
##        0 435 547
##        1 378 490
## 
## , , question = 1
## 
##         gender
## response   0   1
##        0 402 511
##        1 411 526
## 
## , , question = 2
## 
##         gender
## response   0   1
##        0 418 540
##        1 395 497

Visualize three-way contingency table:

pairs(my.tab, diag_panel = pairs_barplot,  labeling=TRUE)

Since we have binary data on three occasions, there are only 8 possible values for the response vector.

# Convert data to wide form
att_wide <- attitudes %>% pivot_wider(
  id_cols = "case",
  names_from = "question",
  names_prefix = "Sit",
  values_from = "response",
  unused_fn = list("gender" = max)
)
n <- dim(att_wide)[1]

The number of rows in wide form tell us the number of subjects questioned: 1850

How many supported legalizing abortion in all three situations, i.e., responded (1,1,1)?

# Create sequence or responses variable:
# Hard way -
att_wide <- att_wide %>%
  mutate(response_seq = case_when(
    Sit1 == 1 & Sit2 == 1 & Sit3 == 1 ~ "111",
    Sit1 == 1 & Sit2 == 1 & Sit3 == 0 ~ "110",
    Sit1 == 1 & Sit2 == 0 & Sit3 == 1 ~ "101",
    Sit1 == 0 & Sit2 == 1 & Sit3 == 1 ~ "011",
    Sit1 == 1 & Sit2 == 0 & Sit3 == 0 ~ "100",
    Sit1 == 0 & Sit2 == 1 & Sit3 == 0 ~ "010",
    Sit1 == 0 & Sit2 == 0 & Sit3 == 1 ~ "001",
    Sit1 == 0 & Sit2 == 0 & Sit3 == 0 ~ "000",
    )
  )
# Easy way -
att_wide <- att_wide %>%
  mutate(response_seq = paste0(Sit1, Sit2, Sit3))

xtabs(~gender + response_seq, data = att_wide)

##       response_seq
## gender 000 001 010 011 100 101 110 111
##      0 356  19  21   6  32  11  26 342
##      1 457  22  18  14  47  14  25 440

# Note that the majority of the sample was either 000 or 111.
# --> high association between repeated measurements.
# Sample correlations among repeated measurements:
cor(att_wide[,2:4])

##           Sit1      Sit2      Sit3
## Sit1 1.0000000 0.8248410 0.7958923
## Sit2 0.8248410 1.0000000 0.8312609
## Sit3 0.7958923 0.8312609 1.0000000

Modeling: GLMM

### Fit GLMM by ML ###
# No interaction between gender and question -->
# gender effect assumed to be identical for each situation.
mod1 <- glmer(response ~ gender+question+(1|case), 
family=binomial,
    data=attitudes, nAGQ = 25)
# Larger value of "nAGQ" =  greater accuracy in approximation
# --> longer computational time 
# Default is 1 (Laplace approximation) - gives very different results!
summary(mod1)

## Generalized linear mixed model fit by maximum likelihood (Adaptive
##   Gauss-Hermite Quadrature, nAGQ = 25) [glmerMod]
##  Family: binomial  ( logit )
## Formula: response ~ gender + question + (1 | case)
##    Data: attitudes
## 
##      AIC      BIC   logLik deviance df.resid 
##   4582.4   4615.5  -2286.2   4572.4     5545 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78590 -0.11563 -0.09964  0.13252  1.71897 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  case   (Intercept) 87.2     9.338   
## Number of obs: 5550, groups:  case, 1850
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.64517    0.40930  -1.576   0.1150    
## gender       0.01363    0.53197   0.026   0.9796    
## question1    0.84107    0.16063   5.236 1.64e-07 ***
## question2    0.29487    0.15733   1.874   0.0609 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) gender qustn1
## gender    -0.728              
## question1 -0.202  0.000       
## question2 -0.196  0.000  0.508

summary(update(mod1, nAGQ = 1))

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: response ~ gender + question + (1 | case)
##    Data: attitudes
## 
##      AIC      BIC   logLik deviance df.resid 
##   4942.7   4975.8  -2466.3   4932.7     5545 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6419 -0.1794 -0.1614  0.1996  1.6073 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  case   (Intercept) 31.24    5.59    
## Number of obs: 5550, groups:  case, 1850
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.46922    0.21722  -2.160   0.0308 *  
## gender       0.01073    0.28176   0.038   0.9696    
## question1    0.59851    0.08958   6.681 2.37e-11 ***
## question2    0.20949    0.08872   2.361   0.0182 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) gender qustn1
## gender    -0.727              
## question1 -0.207  0.000       
## question2 -0.207  0.000  0.503

mod2 <- glmer(response ~ gender*question+(1|case), family=binomial,
    data=attitudes, nAGQ = 25)
# Test for interaction effect between gender and question:  
anova(mod2)

## Analysis of Variance Table
##                 npar  Sum Sq Mean Sq F value
## gender             1  0.0019  0.0019  0.0019
## question           2 28.3869 14.1934 14.1934
## gender:question    2  0.9976  0.4988  0.4988

anova(mod1,mod2)

## Data: attitudes
## Models:
## mod1: response ~ gender + question + (1 | case)
## mod2: response ~ gender * question + (1 | case)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## mod1    5 4582.4 4615.5 -2286.2   4572.4                     
## mod2    7 4585.4 4631.7 -2285.7   4571.4 1.0227  2     0.5997

# EBLUPS same for all individuals with same response sequence
# and same gender, e.g.,
ranef(mod1)$case[att_wide$response_seq=="000", 1]

##   [1] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##   [8] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [15] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [22] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [29] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [36] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [43] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [50] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [57] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [64] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [71] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [78] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [85] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [92] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
##  [99] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [106] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [113] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [120] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [127] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [134] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [141] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [148] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [155] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [162] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [169] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [176] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [183] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [190] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [197] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [204] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [211] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [218] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [225] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [232] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [239] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [246] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [253] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [260] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [267] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [274] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [281] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [288] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [295] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [302] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [309] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [316] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [323] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [330] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [337] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [344] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [351] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [358] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [365] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [372] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [379] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [386] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [393] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [400] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [407] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [414] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [421] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [428] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [435] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [442] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [449] -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040 -3.978040
## [456] -3.978040 -3.978040 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [463] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [470] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [477] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [484] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [491] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [498] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [505] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [512] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [519] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [526] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [533] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [540] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [547] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [554] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [561] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [568] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [575] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [582] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [589] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [596] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [603] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [610] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [617] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [624] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [631] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [638] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [645] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [652] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [659] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [666] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [673] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [680] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [687] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [694] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [701] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [708] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [715] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [722] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [729] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [736] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [743] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [750] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [757] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [764] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [771] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [778] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [785] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [792] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [799] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [806] -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191 -3.967191
## [813] -3.967191

att_wide$gender[att_wide$response_seq=="000"]

##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [445] 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [482] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [519] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [556] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [593] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [630] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [667] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [704] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [741] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [778] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

# Note different answers if we use a different R library:
library(glmmML) # May need to install package
# glmmML function only allows for a random intercept
mod1b = glmmML(response ~ gender+question, cluster=case,
    family=binomial, method="ghq", 
    n.points=70, start.sigma=9, data=attitudes)
summary(mod1b)  

## 
## Call:  glmmML(formula = response ~ gender + question, family = binomial,      data = attitudes, cluster = case, start.sigma = 9, method = "ghq",      n.points = 70) 
## 
## 
##                 coef se(coef)        z Pr(>|z|)
## (Intercept) -0.61874   0.3777 -1.63839 1.01e-01
## gender       0.01259   0.4888  0.02575 9.79e-01
## question1    0.83470   0.1601  5.21347 1.85e-07
## question2    0.29240   0.1567  1.86622 6.20e-02
## 
## 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

Modeling: GEE

mod2 <- gee(response ~ gender+question, family=binomial,
            corstr="exchangeable", id=case, data=attitudes)

## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27

## running glm to get initial regression estimate

##  (Intercept)       gender    question1    question2 
## -0.125407576  0.003582051  0.149347113  0.052017989

summary(mod2)   

## 
##  GEE:  GENERALIZED LINEAR MODELS FOR DEPENDENT DATA
##  gee S-function, version 4.13 modified 98/01/27 (1998) 
## 
## Model:
##  Link:                      Logit 
##  Variance to Mean Relation: Binomial 
##  Correlation Structure:     Exchangeable 
## 
## Call:
## gee(formula = response ~ gender + question, id = case, data = attitudes, 
##     family = binomial, corstr = "exchangeable")
## 
## Summary of Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.5068644 -0.4825396 -0.4687095  0.5174604  0.5312905 
## 
## 
## Coefficients:
##                 Estimate Naive S.E.     Naive z Robust S.E.    Robust z
## (Intercept) -0.125325730 0.06782579 -1.84775925  0.06758212 -1.85442135
## gender       0.003437873 0.08790630  0.03910838  0.08784072  0.03913758
## question1    0.149347107 0.02814374  5.30658404  0.02973865  5.02198729
## question2    0.052017986 0.02815145  1.84779075  0.02704703  1.92324179
## 
## Estimated Scale Parameter:  1.000721
## Number of Iterations:  2
## 
## Working Correlation
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.8173308 0.8173308
## [2,] 0.8173308 1.0000000 0.8173308
## [3,] 0.8173308 0.8173308 1.0000000

Note since GEE is really an estimating procedure (similar to least squares) and there is no likelihood, we cannot carry out likelihood ratio tests to compare models.