[ADMB Users] Problems with RE model, using FUNCTIONs for integration

Fri Nov 30 07:48:42 PST 2012

Hi Hans,

Again, thank you for taking the time to respond to me.  Your
suggestion is a good one, and is exactly what I was trying to do in my
second .tpl file I mentioned (kidneyB.tpl).  In that file, I had
written my own small function trapz to do the integration (I believe
the one included with ADMB is called trapzd, sorry if this caused
confusion).  And it did not work fine!  That's why I was hoping to get
feedback, whether the problem is something simple I did wrong, or if
it is a bigger problem.

If you have a moment I would really appreciate if you (or anyone with
expertise) could try to build and run my kidneyB.tpl file and let me
know if you have any insights.  I've attached the .tpl, .dat, and .pin
files to this email just in case.  I put a few extra comments in the
.tpl file to make it clear what I've done.

Regards,

Mark

P.S. I got your follow-up message about getting details of the model.
I will send you something shortly.

On Thu, Nov 29, 2012 at 11:34 PM, Hans J. Skaug <Hans.Skaug at math.uib.no> wrote:
> Hi,
>
>
>
> OK, so you are using the Weibull, but trying to use numerical integration
> instead of the analytical formula.
>
> The point is that both adromb() and trapz() are both under development, and
> hence unstable. You should
>
> start without them, and switch when they are mature. IN the mean time you
> can code up
>
> your own crude numerical integration rule at the point where you make the
> call “Ht = adromb(&model_”.
>
> That should work fine, at least for a start.
>
>
>
> However, I would recommend that your code up at joint model with the
> analytical Weibull Ht.
>
> You should make sure that this model runs well for the amount of data you
> imagine. There are choices
>
> to be made. You must probably read about SEPARABLE_FUNCTION in the user
> manual.
>
>
>
> Hans
>
>
>
> From: Mark Wolters [mailto:kramsretlow at gmail.com]
> Sent: Friday, November 30, 2012 4:32 AM
> To: H. Skaug
> Subject: Re: [ADMB Users] Problems with RE model, using FUNCTIONs for
> integration
>
>
>
> Hi Hans,
>
> Glad to see there is some interest in working on such a model.  ADMB seems
> very suitable as a general tool for fitting this class of models, which is
> why I got started on this in the first place.
>
> I think there might be a slight misunderstanding about the hazard function
> in my test example--unless it's me that's  misunderstanding your comment.
> That kidney example is already very simple, it's just a Cox model with a
> Weibull baseline hazard.  You can get the analytical form of the cumulative
> hazard function, and this is what the original kidney.tpl had in it.  I was
> trying to get the cumulative hazard by integration in order to verify that
> it works. When it comes to those joint models, the hazard will usually have
> time-dependent terms in the exponential part, making numerical integration a
> necessity.
>
> Thanks for your help,
>
> Mark
>
> On 2012-11-29 6:48 PM, "H. Skaug" <hskaug at gmail.com> wrote:
>
> Mark,
>
> It would be great to have examples of a "Joint longitudinal and
> survival models" in
> ADMB. It has been on my mind for a long time.
>
> My main suggestion: why don't you start out with a cumulative hazard
> with a known
> parametric form such as the Weibull?
>
> http://en.wikipedia.org/wiki/Weibull_distribution
>
> This is just to get going. When everything is running smoothly you can
> switch
> to whathever hazard you like.
>
> (sorry if I am misunderstanding here)
>
> Hans
>
>
>
> On Thu, Nov 29, 2012 at 8:53 PM, Mark Wolters <kramsretlow at gmail.com> wrote:
>> Hi all,
>>
>> I am new to ADMB and have spent the last 6 weeks or so trying to get
>> up to speed. My ultimate goal is to use ADMB to fit joint longitudinal
>> and survival models, where a longitudinal model (typically a linear
>> mixed effect model) is combined with a survival (usually proportional
>> hazards) model.  There are different ways to link the two parts of the
>> model, for example by including the same random effects in both parts,
>> or by using the predicted value of the longitudinal response as a
>> covariate in the survival model.
>>
>> The likelihoods in question will generally involve an integral (the
>> cumulative hazard function) that needs to be evaluated numerically.
>> The integrand will usually involve one or more random effects.  And it
>> would be desirable/natural to break up the likelihood calculation into
>> a few sub functions, e.g. one to calculate the hazard function, one to
>> compute the longitudinal response, one to do the integration, etc.
>>
>> I've been trying to work my way up to something like the situation I'm
>> interested in by starting with the kidney data example
>>
>> (http://www.admb-project.org/examples/survival-analysis/weibull-regression-with-censoring)
>> and adding features a little at a time.  Unfortunately I haven't been
>> able to get far without some problems, so I would be grateful if
>> anyone could provide some guidance. Below I describe two failed
>> attempts; TPL files for these two cases are attached (the .pin and
>> .dat files can be found at the above web link).
>>
>> I am hoping I have just made some rookie mistakes (also I'm not much
>> of a C++ programmer so forgive me if I'm missing something obvious).
>> Thanks in advance,
>>
>> Mark Wolters
>>
>>
>> ## Problem 1: using adromb( ) ## (file kidneyA.tpl)
>>
>> For my first test, I decided to use adromb( ) to evaluate the
>> cumulative hazard numerically, rather than using the exact form of the
>> cumulative hazard.  My modifications to the .tpl look like this:
>>
>> PROCEDURE_SECTION
>>   // some unchanged code here
>>   // replace the line assigning Ht by the following:
>>       Ht = adromb(&model_parameters::integrand,0,t(i,j),8);
>>       Ht *= r*mfexp(eta);
>>   // more lines
>>
>> FUNCTION dvariable integrand(const dvariable& s)
>>   dvariable tmp;
>>   tmp = pow(s,r-1);
>>   return tmp;
>>
>> When I try to build from this tpl file I get the following error message:
>>
>> *** adcomp -r kidneyA
>> g++ -c -O3 -Wno-deprecated -D__GNUDOS__ -Dlinux -DOPT_LIB
>> -DUSE_LAPLACE -fpermissive -I. -I"C:\admb\ADMB-11\include"
>> -I"C:\admb\ADMB-11\contrib" -o kidneyA.obj kidneyA.cpp
>> kidneyA.cpp: In member function 'virtual void
>> df1b2_parameters::user_function()':
>> kidneyA.cpp:166:55: error: no match for 'operator=' in 'Ht =
>> function_minimizer::adromb(PMF, double, double,
>> int)(&model_parameters::integrand, 0.0,
>>
>> ((df1b2_parameters*)this)->df1b2_parameters::<anonymous>.df1b2_pre_parameters::<anonymous>.model_parameters::<anonymous>.model_data::t.data_matrix::<anonymous>.named_dmatrix::<anonymous>.dmatrix::operator()(i,
>> j), 8)'
>>
>>
>> ## Problem 2: using a trapezoidal rule ## (file kidneyB.tpl)
>>
>> In an attempt to make progress, I decided to drop adromb( ) and just
>> use a function to do trapezoidal integration with the hazard function
>> hard-coded in, to see if it would work. My tpl file now includes the
>> following:
>>
>> PROCEDURE_SECTION
>>   // some unchanged code here
>>   // replace the line assigning Ht by the following:
>>       Ht = trapz(0,t(i,j),100); // arguments are integration limits
>> and number of intervals.
>>       Ht *= r*mfexp(eta);
>>   // more lines
>>
>> FUNCTION dvariable trapz(double a, double b, int N)
>>   int j;
>>   dvariable out;
>>   double h = (b-a)/(N-1);
>>   out = pow(a,r-1) + pow(b,r-1);
>>   for(j=1;j<=N-2;j++)
>>   {
>>     out += 2*pow(a+j*h,r-1);
>>   }
>>   out *= h/2;
>>   return out;
>>
>> This version builds successfully, but when I run it I get "nan" values
>> for my random effects, and #QNAN values in my gradient output.
>>
>> _______________________________________________
>> Users mailing list
>> Users at admb-project.org
>> http://lists.admb-project.org/mailman/listinfo/users
>>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: kidneyB.dat
Type: application/octet-stream
Size: 821 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121130/14d02590/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: kidneyB.pin
Type: application/octet-stream
Size: 261 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121130/14d02590/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: kidneyB.tpl
Type: application/octet-stream
Size: 2901 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121130/14d02590/attachment-0002.obj>