[ADMB Users] crappy software

Mark Maunder mmaunder at iattc.org
Tue Nov 1 08:47:46 PDT 2011


Anyone want to try this using BUGS?

-----Original Message-----
From: users-bounces at admb-project.org [mailto:users-bounces at admb-project.org] On Behalf Of dave fournier
Sent: Tuesday, November 01, 2011 7:58 AM
To: users at admb-project.org
Subject: Re: [ADMB Users] crappy software

On 11-10-31 12:26 PM, Ben Bolker wrote:

I modified Ben's R script to run both lmer and admb together as Hans suggested. There is not much mystery here. The R routine grossly underestimates the std dev of the parameter estimate for b1.
So the test

     ((b1-1)/std_dev)^2 > 3.84  (approx)

is accepted too often.  Based on MC results the std dev fo the estimate
b1 is
about 0.2.  ADMB estimates it at about 0.17

Here are a few lmer estimates

x             1.0459     0.1161   9.008 <2e-16 ***
x             0.7325     0.1113   6.578 4.76e-11 ***
x             1.1828     0.1356   8.722 <2e-16 ***
x             0.9313     0.1231   7.563 3.94e-14 ***
x             0.9223     0.1055   8.742 <2e-16 ***
x             0.7663     0.1025   7.477 7.63e-14 ***
x             1.0389     0.1333   7.792 6.58e-15 ***
x            0.88643    0.10603    8.36 <2e-16 ***
x             1.0123     0.1144   8.850 <2e-16 ***
x             0.8068     0.1243    6.49 8.57e-11 ***
x             1.0073     0.1201   8.385 <2e-16 ***
x            0.68238    0.10995   6.206 5.43e-10 ***
x             1.0341     0.1106   9.351 < 2e-16 ***
x             1.0262     0.1471   6.977 3.01e-12 ***
x             1.1407     0.1395   8.179 2.85e-16 ***
x             1.0953     0.1151    9.52 <2e-16 ***
x             1.1903     0.1218   9.774 <2e-16 ***
x             1.1091     0.1242   8.927 <2e-16 ***
x             1.3172     0.1657   7.951 1.85e-15 ***
x             1.3353     0.1431   9.329 <2e-16 ***
x             1.0300     0.1365   7.548 4.42e-14 ***
x             1.0520     0.1317    7.99 1.35e-15 ***
x             1.3979     0.1553   9.004 <2e-16 ***
x             1.2640     0.1346   9.391 <2e-16 ***
x             0.9209     0.1227   7.502 6.27e-14 ***
x             0.9244     0.1222   7.568  3.8e-14 ***
x            1.15014    0.11431  10.062 <2e-16 ***

and here are the corresponding ADMB values
    4   beta     1.0968e+00 1.7214e-01
      4   beta     8.3906e-01 1.7347e-01
      4   beta     1.1503e+00 1.7187e-01
      4   beta     8.9949e-01 1.5700e-01
      4   beta     9.7361e-01 1.4663e-01
      4   beta     8.5326e-01 1.5109e-01
      4   beta     1.0206e+00 1.7697e-01
      4   beta     8.7089e-01 1.3210e-01
      4   beta     1.0151e+00 1.5242e-01
      4   beta     8.4811e-01 1.6853e-01
      4   beta     1.0437e+00 1.7211e-01
      4   beta     6.9127e-01 1.4248e-01
      4   beta     1.0806e+00 1.5592e-01
      4   beta     9.7461e-01 1.7540e-01
      4   beta     1.0925e+00 1.7724e-01
      4   beta     1.1178e+00 1.6025e-01
      4   beta     1.1760e+00 1.6975e-01
      4   beta     1.1248e+00 1.7232e-01
      4   beta     1.2356e+00 2.2140e-01
      4   beta     1.3575e+00 2.1031e-01
      4   beta     9.7267e-01 1.6777e-01
      4   beta     1.0119e+00 1.6961e-01
      4   beta     1.3383e+00 2.1261e-01
      4   beta     1.2759e+00 1.8921e-01
      4   beta     9.0485e-01 1.6004e-01
      4   beta     9.4662e-01 1.6704e-01
      4   beta     1.1380e+00 1.5874e-01

This is the batch script needed to run the simulation on Linux.
for i in {1..25}
do
#  ./simulator
  R CMD BATCH ./bolker-par.r
   ./analyzer -mno 10000 -ams 10000000 -noinit
   grep "^ *4 *b" *.std >> b1
done

This is Ben's modified R script.

library("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)
    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)

I attached the  the tpl file for the analyzer modified to get the sign right.

I think anyone can run this simulation using these files.
change for i in {1..25}
to say
for i in {1..1000}

and go do something else. These are for the Laplace approx (I think). Then the scripts can be modified to do AGH








More information about the Users mailing list