[ADMB Users] crappy software

[mintoc at mathstat.dal.ca]

Hi Dave,

Thanks for the code and discussion.

The code:

    ...
    tau=exp(log_tau);
    dvar_matrix ch(1,2,1,2);
    ch(1,1)=1.0;
    ch(1,2)=0.0;
    ch(2,1)=rho;
    ch(2,2)=1.0;
    ch(2)/=norm(ch(2));
    ch*=tau;
    sigma=ch*trans(ch);
    ...

in 'analyzer.tpl' produces a 2X2 matrix in which the variances of both
random effects are constrained to equal tau^2. This is how the data are
generated (in 'simfun') but shouldn't we pretend we don't know this during
estimation?

Might something like (from the ADMB manual)

    ...
    dvar_matrix chQ(1,2,1,2);
    chQ.initialize();
    int ii=1;
    for (int i=1;i<=2;i++)
      for (int j=1;j<=i;j++)
        chQ(i,j)=Qcoff(ii++);
    sigma=chQ*trans(chQ);
    ...

be more comparable to 'glmer' in terms of number of parameters (3) used in
parameterizing the covariance matrix? Perhaps I am missing something and
I'm also not sure if this will make much of a difference but I wanted to
check so that the comparisons are comparable.

Thanks,
Coilin


>
> 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)
>
>
>
>
>
>
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
>