Section 9 Box 7: Estimate fluctuation parameter \(\epsilon\)

eps_mod <- glm(Y.bounded ~ -1 + H.A1L + H.A0L +  
                 offset(qlogis(Pred.initialQ0)), 
               family = "binomial",
               data = ObsData)
epsilon <- coef(eps_mod)  
epsilon["H.A1L"]

##       H.A1L 
## 0.005867462

epsilon["H.A0L"] 

##       H.A0L 
## 0.008207235