[ADMB Users] quick, stupid question: retrieve standard errors for some parameters when bad hessian?

[Ben Bolker]

In the continuing Zhang et al 2011 saga: one of the conditions under
which Zhang et al did the simulations was tau=0.001 (i.e. very low
among-individual variation) -- i.e. the problem really reverts to a
non-mixed-model problem.

When I try to do this one, ADMB perfectly reasonably complains that
the Hessian is not positive definite. Therefore, the std file doesn't
get created. I'd like to try to retrieve the Hessian and get the
inverse of specified columns myself (or do a generalized inverse). I see
that there are files produced:

binary: admodel.hes, analyzer.rhes, hesscheck, hessian.bin

are these documented somewhere, or has someone figured out what they
are, or do I need to work through an example myself to figure it out ... ?

cheers
Ben Bolker

[dave fournier]

On 11-11-02 09:25 AM, Ben Bolker wrote:

That is the least interesting case. However what one can do is to
ocnvince oneself that
when the model wants to set the variance equal to 0 you get the same
estimates by
fixinf the variance at a small number (and fixing th covariance too).

init_vector log_tau(1,2,-1)
init_bounded_number rho(-3.,3.,-1)

then make a pin file with the same values as in the ifnal estimates
and run the model

# log_tau:
-6.29298 -5.67007
# rho:
0.00000000000
# beta:
0 0
# bb:
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0

and now you get the std devs for beta

index name value std dev
1 beta 9.9124e-01 2.1486e-01
2 beta 1.1868e+00 2.4379e-01

So how to convince yourself that these are the same
std devs that you would get for the original model.

You cad add penalties on the parameters

f+=100*square(log_tau(1)+6.2929)

etc in the original model.

this adds a term to the Hessian which will make it positive definite.
You will see that it does not change the estimate for the std dev of the
betas.

The point is that the estimates for beta and log_tau are uncorrelated when
tau is (almost) 0.

_______________________________________________
Users mailing list
Us...@admb-project.org
http://lists.admb-project.org/mailman/listinfo/users

[dave fournier]

Anoither way to do this is to add a tiny quadratic penalty to the
original model

SEPARABLE_FUNCTION void f1(int& i,const dvar_vector& log_tau,const
prevariable & rho,const dvar_vector & bbi,const dvar_vector & beta)
dvar_vector tau=exp(log_tau);
if (i==1)
{
ll+=1.e-6*norm2(log_tau);
ll+=1.e-6*square(rho);
}

You may need to use something like

-crit 1.e-8

I have done enough simulations to see a pattern.

Laplace approx

ADMB-RE 0.12426 169 0.973034 0.159918 0.141409
lmer 0.252941 170 0.985313 0.200662 0.118773

so the type 1 error rates are approx 0.12 and 0.25

and for 5 point AGHQ

ADMB-RE 0.0666667 105 0.986444 0.129637 0.140441
lmer 0.333333 105 0.872585 0.167648 0.107097

Of course these samples are a bit small to say anything definitive.

It would be good to modify the R code so that the random numbers produced
can be controlled so that we all get the same results for the same seed.

[dave fournier]

Just to round this out I modified the ADMB-RE fits from AGHQ to -is 100.
That is importance sampling with 100 points. This is useful because it
can be used for
crossed random effects where AGHQ is impractical.
The type 1 error for ADMB-RE is 5.6%, near the theoretical value.
lmer is 32.4%.

ADMB_RE -is 100 0.056 250 0.994368 0.149674 0.143516
lmer AGHQ=5 0.324 250 0.898267 0.160694 0.10748

[dave fournier]

This turns out to be remarkably difficult and boring to deal with.

The idea is that first you fit the model and then check for very small
variances

tau1 and tau2. If one is small you make it inactive and then
calculate the
Hessian for the rest. Easiest is to run the model twice. Before running
it the first time
get rid of analyzer.par

rm analyzer.par
./analyzer -ind analyzer.dat -mno 100000 -gbs 200000000 -cbs 100000000
-nohess
./analyzer -ind analyzer.dat -mno 100000 -gbs 200000000 -cbs 100000000
-phase 10 -ainp analyzer.par
grep "beta2" *.std >> b1

This extra code is added to the model to determine what variables to
turn off

DATA_SECTION
init_int n
int lflag1
int lflag2
init_matrix x(1,n,1,3)
init_matrix Y(1,n,1,3)
LOC_CALCS
lflag1=2;
lflag2=2;
double lt1,lt2;
lt1=0;
lt2=0;
cifstream cif("analyzer.par");
if (!(!cif))
{
cif >> lt1 >> lt2;
if (lt1<-3.9)
lflag1=-1;
if (lt2<-3.9)
lflag2=-1;
}

and

init_bounded_number log_tau1(-4.0,4.0,lflag1)
init_bounded_number log_tau2(-4.0,4.0,lflag2)

I attached the tpl file for the modified model.