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

[Mark Maunder]

Mark,

You might want to make  h a dvariable so the derivatives are carried through.

Mark

-----Original Message-----
From: users-bounces at admb-project.org [mailto:users-bounces at admb-project.org] On Behalf Of Mark Wolters
Sent: Thursday, November 29, 2012 11:53 AM
To: users at admb-project.org
Subject: [ADMB Users] Problems with RE model, using FUNCTIONs for integration

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.