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.<div>
<br></div><div>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.</div>
<div><br></div><div>I hope you can take a look at it, and discover my mistake.</div><div><br></div><div><br></div><div>Thanks,</div><div><br></div><div>Chris Gast</div><div><br></div><div><br><div><br>-----------------------------<br>
Chris Gast<br><a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a><br>
<br><br><div class="gmail_quote">On Wed, Mar 16, 2011 at 2:33 AM, H. Skaug <span dir="ltr"><<a href="mailto:hskaug@gmail.com">hskaug@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Chris,<br>
<br>
I think everything is OK with the variance formula in ADMB.<br>
<br>
Informal proof:<br>
I prefer to work with the much simpler example:<br>
<br>
datum = b + u, where b=parameter, u=random effect and datum=0 is<br>
observed (see attached files).<br>
<br>
When b is fixed at its MLE b=0 I get the std-file (ignore dummy)<br>
<br>
index name value std dev<br>
1 dummy 2.0000e-04 1.0000e+00<br>
2 u 0.0000e+00 7.0711e-01<br>
<br>
When b is estimated:<br>
<br>
index name value std dev<br>
1 b 0.0000e+00 1.4142e+00<br>
2 dummy 2.0000e-04 1.0000e+00<br>
3 u 0.0000e+00 1.0000e+00<br>
<br>
so that there is an inflation in the std of the random effect (u)<br>
caused by the second term in (3.1) in the re-manual.<br>
<br>
<br>
Hans<br>
<div><div></div><div class="h5"><br>
<br>
<br>
<br>
<br>
On Wed, Mar 16, 2011 at 6:28 AM, Chris Gast <<a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a>> wrote:<br>
> For anyone who may be following along, I believe I may have found the<br>
> answer. By combining my two previous posts (where I try to estimate var(p1)<br>
> and var(g)), I was able to verify ADMB's estimate of var(p1), however this<br>
> has lead me to discover that the printed estimated in the .std and .cor<br>
> files is incorrect.<br>
> This is what I did:<br>
> While estimating var(g) [where g=exp(t1)], I used the delta method and<br>
> obtained a number that disagreed with ADMB's estimate. I hypothesized that<br>
> this was because the printed estimate for SE(t1) in the .std and .cor file<br>
> was only the first of the two components from eqn 3.1 of the ADMB-RE manual.<br>
> That is, I hypothesized that ADMB printed SE(t1), but it actually used eqn<br>
> 3.1 when using the delta method to compute variances of functions involving<br>
> t1.<br>
> So SE(g) was estimated in ADMB as .057166 (see prior email), and my<br>
> hypothesis was that<br>
> SE(g) \approx sqrt{ exp(2*t1)*(var(t1) + X)}<br>
> where X = du/dtheta^T * Cov(theta) * du/dtheta is the second term of eqn<br>
> 3.1, which was not included in the SE column for t1 in the .cor and .std<br>
> files.<br>
> Given all of the other parts in the equation for SE(g) in the equation<br>
> above, I solved for X, and found, in this case, that X=.003009676. When I<br>
> recompute<br>
> Var(t1) = ADMB's printed var(t1) + .003009676, and then redo the delta<br>
> method myself, I obtain the exact same variance estimate for p1 as ADMB<br>
> does.<br>
> Since I can't output the sensitivity matrix (du/dtheta) myself, I must<br>
> assume that the second component of eqn 3.1 in the ADMB RE manual for Cov(u)<br>
> is computed as prescribed (that is, that .003009676 is the correct number).<br>
> What this probably-confusing discussion amounts to is that the printed<br>
> variance estimates for RE terms appears to only account for the<br>
> inverse-Hessian component of variability, and do not by default add the<br>
> contribution of variability due to estimation of the fixed parameters. This<br>
> does not seem to be a problem when adding your own functions to ADMB code as<br>
> sdreport_ variables, as ADMB still uses both components in the delta method<br>
> variance approximation. However, this is problematic when someone such as<br>
> myself needs to specify a sampling distribution for portions of such a<br>
> function. My example with p1 is going to be extended to use a<br>
> Horvitz-Thompson estimation approach to compute N1=x1/p1, where I can<br>
> estimate the variability for p1 within ADMB, but I would like to specify the<br>
> sampling distribution for x1 outside of ADMB (in some post-processing R<br>
> code).<br>
> I'm a modestly capable C/C++ programmer, and could probably spend some time<br>
> digging around in the source to figure out how to print both components of<br>
> the variance estimate for REs (I've dabbled in the source a bit already to<br>
> change the precision of printed values), but I suspect there is a much more<br>
> knowledgeable individual who might be able (and willing?) to accomplish<br>
> this.<br>
> Is there a specific reason that printed SE values for REs omit this<br>
> component? Would it not be better if, by default, the printed value for<br>
> SE(t1) was equal to the square root of eqn 3.1 in the manual---or perhaps if<br>
> both components were available in some output file?<br>
><br>
><br>
> Chris Gast<br>
><br>
</div></div></blockquote></div><br></div></div>