[ADMB Users] SE Estimates including process error in ADMB

[Chris Gast]

Thanks for your reply, Hans. I see that it does appear to work in the simple
example you attached, but I still cannot find an explanation for why my
example is not working.  I have tried to "break" your example in the same
manner in which mine is "broken", but it seems to hold up.

I'll send you the relevant files for my example--I would put them on the
list, but they will not satisfy the 40kb size limit on attachments.
 Basically, I use ADMB to estimate var(exp(u)), and I try to approximate it
myself with the delta method, using exp(2*u)*var(u) (which is what I am led
to believe ADMB is using internally), and I do not get the same result. My
initial post (the p1 function, in which I am actually interested) was more
complicated, but suffers from the same problem.

I hope you can take a look at it, and discover my mistake.


Thanks,

Chris Gast



-----------------------------
Chris Gast
cmgast at gmail.com


On Wed, Mar 16, 2011 at 2:33 AM, H. Skaug <hskaug at gmail.com> wrote:

> 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 --------------
An HTML attachment was scrubbed...
URL: <http://lists.admb-project.org/pipermail/users/attachments/20110316/60106500/attachment.html>