[ADMB Users] crappy software

dave fournier davef at otter-rsch.com
Tue Nov 1 08:59:33 PDT 2011


What's bugs got to do with it? The point is to compare point estimates and
estimated std devs from two models using the same procedure (presumably
I don't see why the point estimates are so different.  Way back when 
Hans and I compared
point estimates for ADMB and SAS NLMIXED they were essentially
identical.  I'm still slightly suspicious of my conversion of the data to a
format suitable for ADMB. R sucks because you can use the transpose of
a non square matrix instead of the matrix by mistake and it just
quietly screws you.)

 > 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=tau2*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