Calling ADMB-RE from R |
AD model Builder has existed for years as a program to produce stand alone executables on Winows and Linux. Modifying it to seamlessly produce shared libraries for R can be expected to produce a few wrinkles so
please be patient and give us feedback.
To run the examples in R you need to download (and source into R) the file
glmmadmb.s which defines the driver
function glmm.admb() and a modified versions "epil2" of the epil-dataset from MASS.
You also need to download the library file (nbmm.dll if you
run Windows, or nbmm.so
if you run Linux) and save it in
the directory where you run R.
You should note that the ADMB-RE executables create temporary files (sometimes large), so you should probably start R in a
specially dedicated directory.
Standard deviations of parameter estimates can be found in the file "nbmm.std"
| |
|
Model description
The negative binomial distribution can be used instead of the Poisson distribution to investigate whether
there is overdispersion in the data, that is whether the variance of the
observations is greater than that which would be expected for a Poisson
distribution. Parameter estimation for such models is
generally claimed to be difficult. See for example R-help
the mailing list archives of the statistical modeling language R. The data
used in this example are the epilepsy data considered in Venables and
Ripley Modern applied
statics with S 4th edition. and by Booth et al. Negative Binomial
Loglinear Mixed Models.
Implementation in ADMB-RE callable from R
We coded up the model in ADMB-RE (nbmm.tpl)
with flexible linear predictors for both fixed and random effects. The program was then compiled into a DLL that can be called from R
via the R-function glmm.admb(). Examples of how to use this function are given below. Currently glmm.admb()
only allows negative binomial, but implementing other distributions like Bernoulli and Poisson is just a question of adding
a few lines of code to nbmm.tpl.
Comparison with SAS NLMIXEDBooth et al. attempt
to fit two negative binomial loglinear mixed models to the data. They refer to these models as the full model and a simpler model.
For the full model they report:
The fit of the full negative binomial model using NLMIXED was
very unstable. Different starting values led to different estimates and
very different standard errors.
Booth et al also apply a Monte Carlo EM algorithm (MCEM) to the full model and report:
Application of the MCEM algorithm in this problem suggest that
the random slope is 0. The MCEM algorithm was run for a large number of
iterations with all of the estimates except for slope variance and the
covariance converging quickly. These latter two estimates appear to be
slowly converging toward 0.
The full model of Booth et al. is specified as:
glmm.admb(y~Base*trt+Age+Visit,random=~Visit,group="subject",data=epil2)
This model converges quickly (30 seconds), with the ML estimate of the variance of the random slope being equal to zero (or extremely small),
and as a consequence of this there is very little information about the correlation parameter (between the random intercept and the slope).
The standard eviations of the parameter estimates including those of the random effects are found in the file nbmm.std.
We also fitted the simpler model of Booth et al.:
glmm.admb(y~Base*trt+Age+Visit,random=~1,group="subject",data=epil2)
This model converged quickly (15 seconds) to the ML estimates. We used different starting
values to investigate the stability of the model and found that it converged to the same values each time (provided that the chosen
initial values exceeded a minimum level of overdispersion).
Thus it appears that the performance of ADMB-RE is superior to SAS NLMIXED for this problem.
|