[ADMB Users] crappy software

[Ben Bolker]

On 11-11-02 11:14 AM, mintoc at mathstat.dal.ca wrote:
> 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

  My two cents would be that I think you're right.

   Can you give me a version/page number reference for the manual?  I'm
not finding anything but I'm probably just looking in the wrong place ...

   Ben

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