[R-sig-ME] Fwd: Request for help using a generalized linear mixed model in correlated data

Maria Eugenia Castellanos mec@@te||@no@8 @end|ng |rom gm@||@com
Tue Jan 7 18:23:12 CET 2020

Dear R-sig-mixed-models group,

Happy new year and hope you are fine. I will appreciate a lot your help in
this analysis in which I have spent several weeks.

I am an epidemiologist, trying to estimate risk differences between the
risk of tuberculosis (TB) in  2 groups, contacts of TB cases and contacts
of community controls.

My sample is 1043 contacts of controls and 1002 contacts of cases.

The variable index_contact1  (exposure) indicates if an observation comes
from a contact of a TB case or it is a contact of a community control
(coded 1 and 0 respectively). I want to adjust by age  (enr_age, continuous
variable), sex (enr_sex coded 2=women, 1=men) and HIV status (hiv,
1=infected, 0=non infected).

My binary outcome is called ‘tst’ and it is coded as 1 (infected) or 0 (not

One of the reviewers asked me to account for the clustering of these
contacts. So I have a variable called enr_lnkid that indicates the ID of a
TB case (n=122) or a community control (n=124).

All my variables are coded as factors, except for tst and for enr_age,
which they are coded as numeric. I had to write tst as a numeric variable
in order to work with the regression models.

First, I used a GEE model and it worked:

allnoc <- gee(tst ~ index_contact1 + enr_age + enr_sex + hiv, id=netid,
family=binomial ('identity'), corstr = "exchangeable", data = c3)



gee S-function, version 4.13 modified 98/01/27 (1998)


Link:                      Identity

 Variance to Mean Relation: Binomial

 Correlation Structure:     Exchangeable


gee(formula = tst ~ index_contact1 + enr_age + enr_sex + hiv,

    id = netid, data = c3, family = binomial("identity"), corstr =

Summary of Residuals:

       Min         1Q     Median         3Q        Max

-0.9229960 -0.4073128 -0.2338467  0.5127307  0.8636212


                  Estimate   Naive S.E.    Naive z  Robust S.E.   Robust z

(Intercept)      0.1991535 0.0254609144  7.8219323 0.0254041176  7.8394200

index_contact11  0.1600779 0.0208047552  7.6942957 0.0210721716  7.5966513

enr_age          0.0086733 0.0007726673 11.2251418 0.0009238111  9.3886077

enr_sex2        -0.1061413 0.0208703279 -5.0857503 0.0206749937 -5.1337997

hiv1            -0.0393979 0.0414949448 -0.9494627 0.0418152072 -0.9421908

Estimated Scale Parameter:  1.002437

Number of Iterations:  3

Working Correlation


[1,]    1


The result was 0.16 (16%), which is very close to the crude risk difference
(approx.. 14%).

I am trying to replicate my findings, using a generalized linear mixed
model, so I use this code:

  Model1 <- glmer(tst ~  index_contact1 + enr_age + enr_sex + hiv +  (1 |
enr_lnkid), family=binomial, data = c3, control = glmerControl(optimizer =
"bobyqa"), nAGQ=0)

And it works

Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 0) ['glmerMod']

 Family: binomial  ( logit )

Formula: tst ~ index_contact1 + enr_age + enr_sex + hiv + (1 | enr_lnkid)

   Data: c3

      AIC       BIC    logLik  deviance  df.resid

 2584.565  2618.304 -1286.283  2572.565      2039

Random effects:

 Groups    Name        Std.Dev.

 enr_lnkid (Intercept) 0.7106

Number of obs: 2045, groups:  enr_lnkid, 246

Fixed Effects:

    (Intercept)  index_contact11          enr_age         enr_sex2

       -1.39891          0.75106          0.04386         -0.48146

But my results come as link=logit whereas I need to have link=identity to
get risk differences.

I tried then this, using Poisson as the family binomial did not accept


  Model2 <- glmer(tst ~  index_contact1 + enr_age + enr_sex + hiv + (1 |
enr_lnkid), family=poisson(link=identity), data = c3, control =
glmerControl(optimizer = "bobyqa"), nAGQ=0)

But I get this error:

Error in (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L,
maxit = 100L,  :

  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

I want to ask two things:

   1. Is there a way that I can convert the results for model 1 from
   link=logit to link=identity?
   2. Or how I can solve these error with the Poisson family?

Thanks for the help you can provide me, thank you!

Maria Eugenia Castellanos

Global Health Institute

College of Public Health

University of Georgia

Athens, GA, USA

	[[alternative HTML version deleted]]

More information about the R-sig-mixed-models mailing list