[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