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

[H. Skaug]

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
>