[ADMB Users] crappy software
dave fournier
davef at otter-rsch.com
Tue Nov 1 07:57:39 PDT 2011
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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: analyzer.tpl
URL: <http://lists.admb-project.org/pipermail/users/attachments/20111101/edeba866/attachment.ksh>
More information about the Users
mailing list