comparing results across software packages

Price, Emily <ep311508 at ...> writes:

Dear R-sig-mixed-models,

My colleague and I are trying to verify results of running
hierarchical generalized linear models in STATA and R.  In STATA to
designate the number of trials we used binomial(10). To do this in R
we created a variable that was 10 for every observation call trial.
We ran the model two ways in R with and without the offset
command. The following are our commands in STATA and R respectively:

Xtmelogit dv classb trt_lam, || id_lam: trt_lam, 

covariance(unstructured) binomial(10)

glm22 <- glmmadmb(formula=dv~trt.f + classb.f + offset(trial) +
 (trt.f|id.f), data=lam_hlm, family="binomial",
 link="logit",corStruct="full")

glm_off <- glmmadmb(formula=dv~trt.f + classb.f + (trt.f|id.f), 

data=lam_hlm, family="binomial", link="logit",corStruct="full")

Our results for the fixed effects, random effects and log likelihood
are different between R and STATA for both R models. We are not sure
why they are so different and were wondering if this issue has been
encountered before and if additional assistance could help explain
this.

Emily Price and Christine Crumbacher

Three comments here:

(1) I don't think a logit-offset is the way to add information about
the number of trials.  The standard R approach is to use

glmmadmb(formula=cbind(dv,trial-dv)~trt.f + classb.f + 
  (trt.f|id.f), 
 data=lam_hlm, family="binomial", link="logit",corStruct="full")

(2) There is a bug in the current (and previous!!) versions of glmmADMB,
which I haven't dealt with yet -- it doesn't really fit unstructured
variance-covariance matrices, it always defaults to diagonal
v-cov matrices.  You can wait a few days (I will bump this up the
priority list); in the meanwhile it might be worth comparing with
the diagonal covariance structure to see if you get the right
STATA/R comparison.

(3) lme4 should also be able to handle this model, either with


glmmadmb(cbind(dv,trial-dv)~trt.f + classb.f + 
  (trt.f|id.f), 
 data=lam_hlm, family=binomial(link="logit"))

(unstructured variance-covariance matrix is the default) or

glmmadmb(dv/trial~trt.f + classb.f + 
  (trt.f|id.f), weights=trial,
 data=lam_hlm, family=binomial(link="logit"))

(I think glmmML can also do this one).

  Ben Bolker