[ADMB Users] SE Estimates including process error in ADMB

[Chris Gast]

Hans was kind of enough to look at my model and detect my error.  I was
using random_effects_bounded_vectors instead of just random_effects_vectors;
the latter works as prescribed, but the former does not (which is why it is
not recommended to be used, nor mentioned in any documentation).

I encounter other optimization issues when not using bounded RE vectors that
has forced me to use them in the past, but I will now continue my search for
a way to avoid using them.

Thanks to Hans and any others who may have been thinking about this issue.


Chris




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


On Wed, Mar 16, 2011 at 7:50 AM, Chris Gast <cmgast at gmail.com> wrote:

> 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/b6722c1a/attachment.html>