[ADMB Users] crappy software
Mark Maunder
mmaunder at iattc.org
Tue Nov 1 09:05:32 PDT 2011
Just a useful comparison of methods partly to show how to convert from BUGS to ADMB or vice-versa.
-----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 9:00 AM
To: users at admb-project.org
Subject: Re: [ADMB Users] crappy software
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 > > > > > >
_______________________________________________
Users mailing list
Users at admb-project.org
http://lists.admb-project.org/mailman/listinfo/users
More information about the Users
mailing list