[ADMB Users] crappy software

[dave fournier]

I should set up my own server at home and talk to myself.  Oh wait ...

Anyway I did 250 simulation with 5 point AHGQ.  I got the same
results as Zhang et al in that lmer did not get any better. In fact
is is arguably worse so maybe the AGHQ is a bit wonky as well.

The type 1 error for rate lmer is .304. there also appears to be some
bias as the mean parameter estimate is .8977

lmer      0.304 250 0.897707 0.165471 0.108108


In contrast admb-re   has a type 1 error rate of 0.048
very close to the theoretical value of 0.05
  The calculated mean std dev is also almost identical to the std dev of
the parameter estimates  i.e around 0.16


admb-re      0.048 250 0.985706 0.165821 0.162532

It would be nice if another person with obviously nothing useful to
occupy their time could verify the R stuff.  I'm not competent to
do serious R work.

I won't bother you again for a while.

Linux bash script to run the simulations

for i in {1..250}
do
#  ./simulator
  R CMD BATCH ./bolker-par.r
   ./analyzer -mno 10000 -ams 10000000 -noinit -gh 5 -gbs 200000000 -cbs 
100000000
   grep "^ *4 *b" *.std >> b1
done

R script bolker-par

ibrary("lme4")
simfun <- function(n,beta=c(1,1),tau=2.0,
                    Cor=matrix(c(1,0.25,0.25,1),nrow=2)) {
     D <- expand.grid(id=factor(1:n),t=1:3)
     D$x <- rnorm(3*n)  # x_{it} ~ N(0,1)
     X <- model.matrix(~x,data=D)
     Z <- model.matrix(~id-1+id:x,data=D)
     b <- MASS::mvrnorm(n,mu=c(0,0),Sigma=tau^2*Cor)
     D$Y <- rbinom(nrow(D),
                   prob=plogis(X %*% beta + Z %*% c(b)),
                   size=1)
     D
}
   dat=simfun(500)
   fm <- glmer(Y~ 1 + x + (1 + x|id), data=dat,family=binomial,nAGQ=5)
    sink("r.out")
    fm
    sink()
    system("grep \"^x   \"    r.out >> xvalues ",intern=TRUE)
   sdat=dat[sort.list(dat[,1]), ]
   n=500
   write("#n", file = "analyzer.dat", sep = " ", append = FALSE)
   write(n, file = "analyzer.dat", sep = " ", append = TRUE)
   write("#x", file = "analyzer.dat", sep = " ", append = TRUE)
   write(sdat$x, file = "analyzer.dat", sep = " ", append = TRUE)
   write("#Y", file = "analyzer.dat", sep = " ", append = TRUE)
   write(sdat$Y, file = "analyzer.dat", sep = " ", append = TRUE)