H. Skaug hskaug at gmail.com
Wed Mar 16 02:33:20 PDT 2011

```Chris,

I think everything is OK with the variance formula in ADMB.

Informal proof:
I prefer to work with the much simpler example:

datum = b + u, where b=parameter, u=random effect and datum=0 is
observed (see attached files).

When b is fixed at its MLE b=0 I get the std-file (ignore dummy)

index   name   value      std dev
1   dummy  2.0000e-04 1.0000e+00
2   u      0.0000e+00 7.0711e-01

When b is estimated:

index   name   value      std dev
1   b      0.0000e+00 1.4142e+00
2   dummy  2.0000e-04 1.0000e+00
3   u      0.0000e+00 1.0000e+00

so that there is an inflation in the std of the random effect (u)
caused by the second term in (3.1) in the re-manual.

Hans

On Wed, Mar 16, 2011 at 6:28 AM, Chris Gast <cmgast at gmail.com> wrote:
> For anyone who may be following along, I believe I may have found the
> answer.  By combining my two previous posts (where I try to estimate var(p1)
> and var(g)), I was able to verify ADMB's estimate of var(p1), however this
> has lead me to discover that the printed estimated in the .std and .cor
> files is incorrect.
> This is what I did:
> While estimating var(g) [where g=exp(t1)], I used the delta method and
> obtained a number that disagreed with ADMB's estimate.  I hypothesized that
> this was because the printed estimate for SE(t1) in the .std and .cor file
> was only the first of the two components from eqn 3.1 of the ADMB-RE manual.
>  That is, I hypothesized that ADMB printed SE(t1), but it actually used eqn
> 3.1 when using the delta method to compute variances of functions involving
> t1.
> So SE(g) was estimated in ADMB as .057166 (see prior email), and my
> hypothesis was that
> SE(g) \approx sqrt{ exp(2*t1)*(var(t1) + X)}
> where X = du/dtheta^T * Cov(theta) * du/dtheta is the second term of eqn
> 3.1, which was not included in the SE column for t1 in the .cor and .std
> files.
> Given all of the other parts in the equation for SE(g) in the equation
> above, I solved for X, and found, in this case, that X=.003009676.  When I
> recompute
> Var(t1) = ADMB's printed var(t1) + .003009676, and then redo the delta
> method myself, I obtain the exact same variance estimate for p1 as ADMB
> does.
> Since I can't output the sensitivity matrix (du/dtheta) myself, I must
> assume that the second component of eqn 3.1 in the ADMB RE manual for Cov(u)
> is computed as prescribed (that is, that .003009676 is the correct number).
> What this probably-confusing discussion amounts to is that the printed
> variance estimates for RE terms appears to only account for the
> inverse-Hessian component of variability, and do not by default add the
> contribution of variability due to estimation of the fixed parameters.  This
> does not seem to be a problem when adding your own functions to ADMB code as
> sdreport_ variables, as ADMB still uses both components in the delta method
> variance approximation.  However, this is problematic when someone such as
> myself needs to specify a sampling distribution for portions of such a
> function.  My example with p1 is going to be extended to use a
> Horvitz-Thompson estimation approach to compute N1=x1/p1, where I can
> estimate the variability for p1 within ADMB, but I would like to specify the
> sampling distribution for x1 outside of ADMB (in some post-processing R
> code).
> I'm a modestly capable C/C++ programmer, and could probably spend some time
> digging around in the source to figure out how to print both components of
> the variance estimate for REs (I've dabbled in the source a bit already to
> change the precision of printed values), but I suspect there is a much more
> knowledgeable individual who might be able (and willing?) to accomplish
> this.
> Is there a specific reason that printed SE values for REs omit this
> component?  Would it not be better if, by default, the printed value for
> SE(t1) was equal to the square root of eqn 3.1 in the manual---or perhaps if
> both components were available in some output file?
>
>
> Chris Gast
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: simple_var.dat
Type: application/octet-stream
Size: 1 bytes
Desc: not available
-------------- next part --------------
A non-text attachment was scrubbed...
Name: simple_var.pin
Type: application/octet-stream
Size: 5 bytes
Desc: not available