Thanks again, Hans. The mixed order of tau and mu in the pin file was definitely an error. Fortunately, the effects do not appear to be too great, as they were of the same sign and similar magnitude.<div><br></div><div>As for the commented-out code:</div>
<div><br></div><div>totL += -.5*norm2(u);</div><div>e=tau*u;</div><div><br></div><div>I found it a bit confusing as well, but this is how it is listed in the ADMB-RE .pdf manual that I downloaded. I don't think the order of the transformation matters---it's simply a matter of evaluating from the N(0,1) density, and then reporting the "correct" values of e, following the appropriate transformation after also estimating tau. Especially since my tau is likely to be small (at least in my simulations), it is good to avoid using the "norm2(e/tau)" term for numerical stability. I think the code as written would be equivalent to the reverse-order of the statements.</div>
<div><br></div><div><div>e=tau*u;</div><div>totL += -.5*norm2(u);</div><div><br></div><div>Regardless, I actually seem to find more problems with successful model fit when employing that transformation than without it, for unknown reasons.</div>
<div><br></div><div>I haven't yet finished investigating the issue about ADMB estimating tau=0. Instead of using p_i = 1-exp(-c * f_i), I employed a logistic transformation p_i = exp(-c f_i)/(1+exp(-c * f_i)) in search of numerical stability. This, along with some perturbed simulation parameter values seems to have some mixed success.</div>
<div><br></div><div>I am still encountering a problem while fitting similar models (I am investigating a random effect for survival, capture probability, and both simultaneously) with ADMB trying to find appropriate values for the random effects (such as the "e" vector above) that don't seem to make sense. By using frequent cout statements, I can see that ADMB first guesses relatively small values that would make sense given the assumed mean (0), and estimated SD, tau. But when an "erroneous" (for lack of a better word) value is guessed for an entry in the "e" vector, loglikelihood then takes the value of "nan." Subsequent searches for appropriate values of "e" become increasingly extreme in their distance from 0 (values in the hundreds, and then thousands), which exacerbate the problem (loglikelihood continues to be "nan", and never recovers). Thus, my question persists. Is there a way to control ADMB's evaluations of the random effects vectors such that I can limit its search for estimates of the "e" vector?</div>
<div><br></div><div>I will try your proposal of fixing tau at my simulated value, and see what happens with the random effects estimates.</div><div><br></div><div>Thanks for your help,</div><div><br></div><div>Chris Gast</div>
<div><br></div><div><br></div><div><br clear="all"><br>-----------------------------<br>Chris Gast<br><a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a><br>
<br><br><div class="gmail_quote">On Mon, May 3, 2010 at 12:39 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;">
Hi,<br>
<br>
I tried to run you model, and find that the log-likelihood is 20 units higher<br>
when tau=0 than tau=5 (ca). Hence there seems to be information<br>
about tau in your data, contrary to what I suggested without<br>
having run your model. You must find out why the model<br>
is saying so clearly that tau=0. Maybe something is wrong.<br>
When you fix tau at the value used to simulate data,<br>
and plot true(e) against parfile(e). Do scatter around a straight line?<br>
<br>
Further, order of tau and mu was mixed in your pin file.<br>
And, the code you had commented out<br>
<div class="im"><br>
/*totL += -.5*norm2(u);<br>
e=tau*u;*/<br>
<br>
</div>Here e is being assigned a value after it is being use in the code.<br>
<br>
Hans<br>
<div><div></div><div class="h5"><br>
<br>
<br>
<br>
On Mon, May 3, 2010 at 12:34 AM, Chris Gast <<a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a>> wrote:<br>
> Thanks, Hans, for your thoughtful response.<br>
> I have tried a number of ways to increase the information about tau,<br>
> including those you describe. Success has been mixed--but this does bring<br>
> up a related issue that has caused me great difficulty, and has backed me<br>
> into the corner in which I find myself. When estimating random effects and<br>
> variance components (I typically follow the manual's advice and estimate<br>
> both of these in a second phase), I frequently encounter loglikelihood=NaN<br>
> in this second phase. It seems to be occurring when values for the random<br>
> effects that are too large (in absolute value) are tested. The definition<br>
> of "too large", of course, depends on the parameter estimates and the<br>
> resulting distribution based on the structure (typically,<br>
> normally-distributed) that I have assumed/simulated. Frequently, my vector<br>
> of random effects contains values in the hundreds or even thousands, when I<br>
> have assumed (and simulated) a normal distribution centered around 0, and my<br>
> estimated variance is ~5 (I simulated it to be about 2.5).<br>
> Is anyone aware of a way to exert some control over this search, so I can<br>
> prevent NaN's from occurring, even if at the cost of a less-than-full search<br>
> of the random effects' distribution?<br>
><br>
><br>
><br>
> -----------------------------<br>
> Chris Gast<br>
> <a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a><br>
><br>
><br>
> On Sat, May 1, 2010 at 10:40 PM, H. Skaug <<a href="mailto:hskaug@gmail.com">hskaug@gmail.com</a>> wrote:<br>
>><br>
>> Chris,<br>
>><br>
>> What you are indirectly doing when using<br>
>><br>
>> totL += (.5)*(-(nyears-1)*log(tau)-.5*norm2(e/tau));<br>
>><br>
>> is to place a penalty on tau, favoring large values of tau. (Somewhat<br>
>> like a Bayesian prior on tau.) The penalty you have added is<br>
>><br>
>> totL += (.5)*(nyears-1)*log(tau);<br>
>><br>
>><br>
>> But, your MLE of tau is 0. Since this is happending repeatedly<br>
>> on simulated datasets where you used tau>0 to generate the<br>
>> data, the conclusion is that tau is weakly identified under your<br>
>> model (given parameters and amount of data). This is maybe not<br>
>> good news for you, but you should appreciate ADMB telling you this.<br>
>> To investigate my claim you can increase the information about<br>
>> tau in the simulations, by increasing value of tau and the amount of data.<br>
>><br>
>><br>
>> Hans<br>
>><br>
>><br>
>><br>
>><br>
>> On Sat, May 1, 2010 at 9:44 PM, Chris Gast <<a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a>> wrote:<br>
>> > Hello,<br>
>> > I'm still relatively new to ADMB--that information might be helpful in<br>
>> > diagnosing the problem I seem to be encountering.<br>
>> > The setup is this: I have a relatively simple product-multinomial model<br>
>> > that I can fit well (with simulated data) when all effects are<br>
>> > considered<br>
>> > fixed. The data are age-specific harvest numbers of a fictional animal.<br>
>> > I<br>
>> > am simulating stochasticity in both survival ( logistic transformation<br>
>> > of<br>
>> > \beta+e_i, e_i~N(0,tau^2)), and in my vulnerability parameter. I have a<br>
>> > single vulnerability coefficient c, such that the probability of harvest<br>
>> > is<br>
>> > 1-exp(-(c+e_i)*f_i), for hunter-effort level f_i in year i. I am<br>
>> > currently<br>
>> > working with 12 simulated years of data and three age classes. For the<br>
>> > discussion below, I am limiting myself to fitting the single random<br>
>> > effect<br>
>> > for survival (rather than both random effects for survival and<br>
>> > vulnerability). (I am also simulating 4 years of binomial telemetry<br>
>> > data to<br>
>> > provide information on the probability of harvest.)<br>
>> > The problem is this: ADMB seems to have some trouble estimating the<br>
>> > single<br>
>> > variance component. My resulting report either spits out my initial<br>
>> > value<br>
>> > (I'm using true simulated values as initial estimates), or ADMB forces<br>
>> > the<br>
>> > parameter estimates as close to zero as I will allow it (of course, it's<br>
>> > constrained to be positive). I have tried many avenues of attack from<br>
>> > the<br>
>> > manual and in other examples, and none of them seem to alleviate the<br>
>> > problem. Having run out of options, I tried replacing the "prior"<br>
>> > distribution portion of the log-likelihood with<br>
>> > totL += (.5)*(-(nyears-1)*log(tau)-.5*norm2(e/tau));<br>
>> > instead of the usual<br>
>> > totL += (-(nyears-1)*log(tau)-.5*norm2(e/tau));<br>
>> > and the variance component estimates are actually OK--they are at least<br>
>> > on<br>
>> > the correct order of magnitude, and not precisely equal to my initial<br>
>> > parameter estimate. I don't really understand why this down-weighting<br>
>> > is<br>
>> > successful (or if it is truly successful or I have just gotten lucky<br>
>> > somehow). I would also really like to avoid this seemingly ad hoc<br>
>> > procedure, since I think ADMB is powerful enough to fit this model<br>
>> > without<br>
>> > resorting to such things--it seems to fit more complex models relatively<br>
>> > easily.<br>
>> > For the record, I have tried the transformation procedure of using a<br>
>> > standard N(0,1) density for the random effect, and then scaling up to<br>
>> > match<br>
>> > the true N(0,tau^2) distribution, and in this case, this "trick" seems<br>
>> > to<br>
>> > hasten ADMBs result of loglikelihood=nan, from which it does not<br>
>> > recover. I<br>
>> > have also tried playing with boundary values on constrained parameters,<br>
>> > and<br>
>> > met with no success.<br>
>> > I have pasted a stripped-down version of my .tpl file (the working<br>
>> > version,<br>
>> > with the down-weighting intact), with accompanying .dat and .pin files,<br>
>> > in<br>
>> > the hope that someone would be kind enough to examine them. I am aware<br>
>> > that<br>
>> > there are probably some inefficiencies here, but I'm not yet concerned<br>
>> > about<br>
>> > optimizing the program for run-time. Of course, I'm only posting one<br>
>> > simulated dataset---this phenomenon occurs for hundreds of simulations<br>
>> > in a<br>
>> > row.<br>
>> > Thanks very much,<br>
>> > Chris Gast<br>
>> ><br>
>> ><br>
>> > /********************* start tpl file ******************************/<br>
>> > DATA_SECTION<br>
>> > init_int nyears;<br>
>> > init_int nages;<br>
>> > init_int m;<br>
>> > init_vector effdat(0,nyears-1);<br>
>> > init_int telnyears;<br>
>> > init_vector numtagged(0,telnyears-1);<br>
>> > init_vector numrecovered(0,telnyears-1);<br>
>> > init_imatrix dat(0,nyears-1,0,nages-1);<br>
>> > imatrix fcohort(0,nyears-1,0,nages-1);<br>
>> > imatrix pcohort(0,nages-2,0,nages-1);<br>
>> > imatrix fdat(0,nyears+nages-2,0,nages-1);<br>
>> > int i;<br>
>> > int j;<br>
>> > PARAMETER_SECTION<br>
>> > init_vector logparms(0,m-1,1);<br>
>> > init_bounded_number beta(-2,1.7,1);<br>
>> > init_bounded_number logcmu(-6,1,1);<br>
>> > init_bounded_number logtau(-10,1.5,2);<br>
>> > random_effects_vector e(0,nyears-2,2);<br>
>> > number tau;<br>
>> > vector parms(0,m-1);<br>
>> > number prob;<br>
>> > number prevq;<br>
>> > number prod;<br>
>> > number datL;<br>
>> > number telL;<br>
>> > number auxL;<br>
>> > number extraL;<br>
>> > number s;<br>
>> > number c;<br>
>> > number sprod;<br>
>> > objective_function_value totL;<br>
>> > PROCEDURE_SECTION<br>
>> > const double pi=3.14159265358979323846264338327950288;<br>
>> > double datsum;<br>
>> > int colrank;<br>
>> > datL=0.0; auxL=0.0; telL=0.0; tau=mfexp(logtau);<br>
>> > parms=mfexp(logparms);<br>
>> > c=mfexp(logcmu);<br>
>> > cout << "tau=" << tau << ", c=" << c << "\n";<br>
>> > for(i=0;i<(nyears+nages-1);i++){<br>
>> > datsum=0;<br>
>> > datL+=gammln(parms[i]+1);<br>
>> > prevq=1;<br>
>> > prod=0;<br>
>> > colrank=0;<br>
>> > sprod=1;<br>
>> > for(j=0;j<nages;j++){<br>
>> > if(fdat[i][j]!=-1){<br>
>> > if(i<nages) prob=1-exp(-c*effdat[colrank]);<br>
>> > else if(i>=nages) prob=1-exp(-c*effdat[i-nages+1+j]);<br>
>> > if(colrank==0) s=1;<br>
>> > else if(i<nages)<br>
>> > s=exp(beta+e[colrank-1])/(1+exp(beta+e[colrank-1]));<br>
>> > else if(i>=nages)<br>
>> > s=exp(beta+e[i-nages+1+j-1])/(1+exp(beta+e[i-nages+1+j-1]));<br>
>> > sprod*=s;<br>
>> > datsum+=fdat[i][j];<br>
>> > datL-=gammln(fdat[i][j]+1);<br>
>> > if(colrank>0) prod+=prevq*sprod*prob;<br>
>> > else prod=prob;<br>
>> > //cout << "prod=" << prod << ", prob=" << prob << ", s=" << s <<<br>
>> > ",<br>
>> > sprod=" << sprod << ", prevq=" << prevq << "\n";<br>
>> > datL+=fdat[i][j]*log(prevq*sprod*prob);<br>
>> > colrank++;<br>
>> > prevq*=(1-prob);<br>
>> > }<br>
>> > }<br>
>> > datL+=(parms[i]-datsum)*log(1-prod) - gammln(parms[i]-datsum+1);<br>
>> > }<br>
>> > //now compute radio-telemetry likelihood<br>
>> > for(i=0;i<telnyears;i++){<br>
>> > telL+=gammln(numtagged[i]+1) - gammln(numrecovered[i]+1) -<br>
>> > gammln(numtagged[i]-numrecovered[i]+1) +<br>
>> > numrecovered[i]*log((1-exp(-c*effdat[nyears-1]))) +<br>
>> > (numtagged[i]-numrecovered[i])*log(exp(-c*effdat[nyears-1]));<br>
>> > }<br>
>> > for(i=0;i<(nyears-1);i++){<br>
>> > cout << e[i] << "\t";<br>
>> > }<br>
>> > cout << "\n\ndatL=" << datL << ", telL=" << telL << ", extraL=" <<<br>
>> > extraL<br>
>> > << ", totL=" << -(datL+telL) << "\n";<br>
>> > cout << "cest=" << c << ", betaest=" << beta << ", tauest=" << tau <<<br>
>> > "\n\n\n\n\n\n\n\n\n";<br>
>> > totL=(datL+telL);<br>
>> > totL += (.5)*(-(nyears-1)*log(tau)-.5*norm2(e/tau));<br>
>> > /*totL += -.5*norm2(u);<br>
>> > e=tau*u;*/<br>
>> > totL *= -1;<br>
>> > REPORT_SECTION<br>
>> > report << "s\n" << exp(beta)/(1+exp(beta)) << "\n";<br>
>> > report << "tau\n" << tau << "\n";<br>
>> > report << "cmu\n" << c << "\n";<br>
>> > report << "parms\n" << parms << "\n";<br>
>> > PRELIMINARY_CALCS_SECTION<br>
>> > for(i=0;i<nyears;i++){<br>
>> > effdat[i]/=1000;<br>
>> > }<br>
>> > for(i=0;i<(nages-1);i++){<br>
>> > for(j=0;j<nages;j++){<br>
>> > pcohort[i][j]=-1;<br>
>> > }<br>
>> > }<br>
>> > for(i=0;i<nyears;i++){<br>
>> > for(j=0;j<nages;j++){<br>
>> > fcohort[i][j]=-1;<br>
>> > }<br>
>> > }<br>
>> > //get the partial cohorts first<br>
>> > for(i=0;i<(nages-1);i++){<br>
>> > for(j=0;j<nages;j++){<br>
>> > if(j>i){<br>
>> > pcohort[nages-j-1+i][j] = dat[i][j];<br>
>> > }<br>
>> > }<br>
>> > }<br>
>> > //now get the "full" cohorts<br>
>> > for(i=0;i<nyears;i++){<br>
>> > for(j=0;j<nages;j++){<br>
>> > if(i<=(nyears-nages)) fcohort[i][j]=dat[i+j][j];<br>
>> > else if((j+i)<nyears) fcohort[i][j]=dat[i+j][j];<br>
>> > }<br>
>> > }<br>
>> > //now put all the cohort data together in a single matrix<br>
>> > for(i=0;i<(nyears+nages-1);i++){<br>
>> > for(j=0;j<nages;j++){<br>
>> > if(i<(nages-1)) fdat[i][j]=pcohort[i][j];<br>
>> > else fdat[i][j]=fcohort[i-(nages-1)][j];<br>
>> > }<br>
>> > }<br>
>> ><br>
>> > /********************* end tpl file ******************************/<br>
>> ><br>
>> > /********************* start dat file ******************************/<br>
>> > #nyears<br>
>> > 12<br>
>> > #nages<br>
>> > 3<br>
>> > #parameters<br>
>> > 14<br>
>> > #effort<br>
>> > 379 233 829 1012 807 698 859 137 779 574 460 684<br>
>> > #telnyears<br>
>> > 4<br>
>> > #numtagged<br>
>> > 100 100 100 100<br>
>> > #numrecovered<br>
>> > 13 19 17 11<br>
>> > #harvest data<br>
>> > 28 13 9<br>
>> > 15 8 4<br>
>> > 49 24 13<br>
>> > 53 27 13<br>
>> > 42 24 12<br>
>> > 36 16 9<br>
>> > 33 18 8<br>
>> > 5 3 1<br>
>> > 29 13 8<br>
>> > 18 12 5<br>
>> > 16 7 4<br>
>> > 23 12 5<br>
>> ><br>
>> > /********************* end dat file ******************************/<br>
>> > /********************* start pin file ******************************/<br>
>> > #logparms<br>
>> > 5.45959 5.86647 5.92426 5.85793 5.81114 5.69709 5.5835 5.60212 5.39816<br>
>> > 5.30827 5.34711 5.17615 5.32301 5.11799<br>
>> > #beta<br>
>> > 0.322773<br>
>> > #logtau<br>
>> > -1.20397<br>
>> > #logcmu<br>
>> > -1.60944<br>
>> > #e<br>
>> > 0 0 0 0 0 0 0 0 0 0 0<br>
>> ><br>
>> > /********************* end pin file ******************************/<br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> ><br>
>> > -----------------------------<br>
>> > Chris Gast<br>
>> > <a href="mailto:cmgast@gmail.com">cmgast@gmail.com</a><br>
>> ><br>
>> > _______________________________________________<br>
>> > Users mailing list<br>
>> > <a href="mailto:Users@admb-project.org">Users@admb-project.org</a><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>
><br>
><br>
</div></div></blockquote></div><br></div></div>