[R-sig-ME] comparing results across software packages

[Ben Bolker]

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