[R] simulation study using AD model builder to fit a GLMM under the binomial probit link
Roel de Jong
dejongroel at gmail.com
Fri Nov 11 20:39:32 CET 2005
I recently took Dave Fournier up on his offer to evaluate his AD Model
builder package (http://otter-rsch.com/admodel.html) when fitting a GLMM
under the binomial probit link.
I conducted a simulation study in which I drawed 500 samples each
containing 1500 observations from the following model specification:
y = (intercept*f1+pred2*f2+pred3*f3)+(intercept*ri+pred2*rs)
where pred2 and pred3 are predictors distributed N(0,1)
f1..f3 are fixed effects, f1=-1, f2=1.5, f3=0.5
ri is random intercept with associated variance var_ri=0.2
rs is random slope with associated variance var_rs=0.4
the covariance between ri and rs=0
we have 50 level 2 units, so 30 observations/level 2 unit
I then proceeded with the analysis of the 500 samples with the AD Model
builder package. To check for bias, I calculated the average of the
parameter estimates of the 500 samples and compared them to the true
population parameters. There was virtually no bias:
parameter average parameter estimate true value
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
f1 -1.001 -1.000
f2 1.510 1.500
f3 0.499 0.500
var_ri 0.197 0.200
var_rs 0.396 0.400
Then I checked the coverage with alpha=0.95, where asymmetrical
confidence intervals were calculated for the variance components:
parameter coverage (alpha=0.95)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
f1 .928
f2 .948
f3 .956
var_ri .960
var_rs .984
The coverages are quite good, only the variance of the random slope is
high, which suggests that the associated standard error is too large.
Where AD model builder really shines is the fact that convergence was
reached without problems in all 500 samples, where R alternatives like
lmer and glmmPQL, which use Penalized Quasi Likelihood, tend to run in
computational problems. I therefore highly recommend the software for
analyzing binomial mixed models, and I encourage Dave to add it to his
existing negative binomial package for R.
Regards,
Roel de Jong
More information about the R-help
mailing list