Perhaps more simply, if I add the function (as a sdreport_number variable)<div><br></div><div>g = exp(t1)</div><div><br></div><div>where t1 is a RE term, ADMB computes</div><div><br></div><div>g = exp(0.032993) = 0.0057166</div>
<div><br></div><div>and estimates (in the .std and .cor files)</div><div><br></div><div>SE(g) = 0.057166</div><div><br></div><div>where a standard delta-method approximation would yield</div><div><br></div><div>Var(g) \approx (dg/dt1)^2 * var(t1) = exp(2*0.032993) * (0.0070425)^2 = 5.297957e-05 --> SE(g) = 0.007278707</div>
<div><br></div><div>So clearly, I must be doing something wrong, or have some misconception about what is reported in the .std and .cor files for RE terms. Any thoughts?</div><div><br></div><div><br></div><div><br></div>
<div>Thanks,</div><div><br></div><div>Chris Gast</div><div><br></div><div><br></div><div><br></div><div><br clear="all"><br>-----------------------------<br>Chris Gast<br><a href="mailto:cmgast@gmail.com" target="_blank">cmgast@gmail.com</a><br>
<br><br><div class="gmail_quote">On Tue, Mar 15, 2011 at 12:41 PM, Chris Gast <span dir="ltr"><<a href="mailto:cmgast@gmail.com" target="_blank">cmgast@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="gmail_quote"><div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote>Hello again,<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote>I'm still not sure what I think ADMB is doing and what it is actually doing are the same thing, and I think it's related to the second term on the rhs of eqn (3.1) in the ADMB-RE manual.<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote>As an example, I have a function <br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote>p1 = 1-exp((cmu+t1)*f1)/(1+exp((cmu+t1)*f1))<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote></div>where cmu is a parameter, f1 is a datum, and t1 is a random effect. I have added the function p1 as a sdreport_number in ADMB, and am trying to recompute the variance estimate via the estimates in the .cor file myself (I need to know what ADMB is doing, so I can use these estimates in additional approximations). I cannot seem to recreate ADMB's variance estimate for p1. I have attached a trimmed version of the .cor file (trimmed to relevant portions to keep file size down) and some commented R code wherein I'm trying to verify the SE estimate for p1. Would anyone who be willing to check my variance approximation, and see what I'm leaving out with respect to what ADMB is doing? Either I've made a dumb mistake, or ADMB is not doing what I think it's doing.<div>
<div></div><div><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote>Calculations like this seem to work (ADMB output can be verified) if RE terms are omitted, such as the meanp parameter function (noted in the R and .cor file), which is the same as p1 but with t1=0 and f1=1.4; in addition, my variance estimate for functions with REs tend to be lower than those from ADMB. These two facts may indicate that I probably don't handle the second term in eqn (3.1) of the ADMB-RE manual correctly. My understanding is that the SE value for a RE term (such as t1) is already computed as in eqn (3.1) in the ADMB-RE manual, and that is what is presented in the .std and .cor files. Perhaps this is not the case? Or perhaps there is something wrong with backwards-computing Cov(cmu,t1) = corr(cmu,t1)*SE(cmu)*SE(t1), where the three terms on the rhs are all found in the .cor file? I could do some testing of these theories if I could obtain the individual VCV matrix components from Dave Fournier's previous messages on this subject, but I don't see how that can be done (besides L_xx^{-1}, of course).<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote>Thanks very much,<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><font face="arial, helvetica, sans-serif">Chris Gast</font></div></div></div><div><div></div><div><div class="gmail_quote"><font face="arial, helvetica, sans-serif"><br></font></div><div class="gmail_quote">
<font face="arial, helvetica, sans-serif"><br></font><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote><font color="#888888"><br>
</font><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote><font color="#888888"><br></font><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote>-----------------------------<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"></blockquote>Chris Gast<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
</blockquote><a href="mailto:cmgast@gmail.com" target="_blank">cmgast@gmail.com</a><br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>
<div><div>
<br><br></div><div><div></div><div><div class="gmail_quote">On Tue, Nov 16, 2010 at 4:21 AM, dave fournier <span dir="ltr"><<a href="mailto:davef@otter-rsch.com" target="_blank">davef@otter-rsch.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Chris Gast wrote:<br>
<br>
I think this is discussed in the paper Hans and I wrote.<br>
anyway since uhat(x) is the value of uhat which maximizes L(x,u)<br>
It follows that<br>
<br>
L_u(x,uhat(x))=0 for all x so differentiating wrt x you get<br>
<br>
L_xu + L_uu uhat'(x)=0 so that<br>
<br>
uhat'(x) = - L_uu^{-1} * L_xu<br>
<br>
If there are n x's and m u's then uhat'(x) is an n x m matrix.<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>
Thanks again, Dave. Could you clarify what the general entry of the gradient uhat'(xhat) looks like? I had thought this would be dx_i/du_j, but I think these would all be 1's or 0's (since, for instance , Shat_1 = shat + uhat_1). It is likely that I am misunderstanding something.<br>
<br>
<br>
Thanks,<br>
<br>
Chris<br>
<br>
<br>
<br>
<br>
-----------------------------<br>
Chris Gast<br>
</div><a href="mailto:cmgast@gmail.com" target="_blank">cmgast@gmail.com</a> <mailto:<a href="mailto:cmgast@gmail.com" target="_blank">cmgast@gmail.com</a>><div><br>
<br>
<br>
On Mon, Nov 15, 2010 at 7:00 PM, dave fournier <<a href="mailto:davef@otter-rsch.com" target="_blank">davef@otter-rsch.com</a> <mailto:<a href="mailto:davef@otter-rsch.com" target="_blank">davef@otter-rsch.com</a>>> wrote:<br>
<br>
OK,<br>
<br>
there are inly ar efew thing to play with. for notation let<br>
<br>
F(x,u) be the likelihood for x and u.<br>
<br>
let<br>
<br>
L(x) = int F(x,u) du<br>
<br>
be what we get after integrating out u either exactly or by the<br>
laplace approx.<br>
<br>
and let xhat uhat(xhat) be the miximizing values.<br>
<br>
We have the Hessians L_xx(xhat), and F_uu(xhat,uhat(xhat)<br>
and the gradients uhat'(xhat)<br>
Then the covariance matrix for x,u is assumed to be<br>
<br>
L_xx^{-1} L_xx^{-1} * uhat'(x)<br>
<br>
uhat'(xhat)*L_xx^{-1} F_uu^{-1} +<br>
uhat'(xhat)*L_xx^{-1}*uhat'(x)<br>
<br>
The idea is that u-uhat(xhat) is independent of xhat i.e the only<br>
correlation is between<br>
xhat and uhat(x)<br>
<br>
Put in transposes where necessary.<br>
<br>
For any function of x,u the covariance is computed by the delta<br>
method.<br>
<br>
<br>
<br>
<br>
<br>
_______________________________________________<br>
Users mailing list<br></div>
<a href="mailto:Users@admb-project.org" target="_blank">Users@admb-project.org</a> <mailto:<a href="mailto:Users@admb-project.org" target="_blank">Users@admb-project.org</a>><div><br>
<a href="http://lists.admb-project.org/mailman/listinfo/users" target="_blank">http://lists.admb-project.org/mailman/listinfo/users</a><br>
<br>
<br>
</div></blockquote>
<br>
</blockquote></div><br></div></div></div></div>
</blockquote></div><br>
</div></div></blockquote></div><br></div>