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

[Mark Wolters]

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.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: kidneyA.tpl
Type: application/octet-stream
Size: 2256 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121129/0c9b6c4b/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: kidneyB.tpl
Type: application/octet-stream
Size: 2410 bytes
Desc: not available
URL: <http://lists.admb-project.org/pipermail/users/attachments/20121129/0c9b6c4b/attachment-0001.obj>