Weighted cumulative joint longitudinal-survival models

The do file with all the code from this tutorial is available to download here

The main thing I like about merlin is the ability to quickly implement newly proposed methodology. This is going to be one of those examples.

A recently published paper in Statistics in Medicine by Mauff et al. (2017) extended the cumulative association structure to a weighted formulation, proposing two forms of the weight function, with estimation conducted in a Bayesian framework using MCMC, implemented in the JMBayes package in R (points for providing accessible software!). In this tutorial, I’ll show how to fit the first of their weight functions using merlin. I’ll use the standard Primary Biliary Cirrhosis (PBC) example, as they did in their paper (they also used a second example).

We load the dataset as usual

. use https://www.mjcrowther.co.uk/data/jm_example.dta, clear
(Example dataset for joint modelling, Michael Crowther 2011)

In the below models I’ll use splines to model the fixed effects of time and assume a random intercept and random linear trend, and a Weibull survival submodel. I’ll also use Monte-Carlo integration to integrate out the random effects. All of these things can of course be changed.

Current value

We can start with the standard current value parametrisation

. merlin  (logb rcs(time,df(3)) fp(time,pow(1))#M2[id]@1 M1[id]@1               ///
>                                        , family(gaussian) timevar(time))      ///
>         (stime trt EV[logb] , family(weibull, failure(died)) timevar(stime))  ///
>         , covariance(unstructured) restartvalues(M2 0.01)
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_3
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1

Fitting fixed effects model:

Fitting full model:

Iteration 0:   log likelihood = -3087.2721  (not concave)
Iteration 1:   log likelihood = -2235.4369  
Iteration 2:   log likelihood = -1938.8448  
Iteration 3:   log likelihood = -1918.2626  
Iteration 4:   log likelihood = -1909.4656  
Iteration 5:   log likelihood = -1909.0985  
Iteration 6:   log likelihood = -1909.0976  
Iteration 7:   log likelihood = -1909.0976  

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood = -1909.0976
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |            
     rcs():1 |   .1194282   .0247549     4.82   0.000     .0709095     .167947
     rcs():2 |  -.0100245   .0041173    -2.43   0.015    -.0180943   -.0019547
     rcs():3 |    .002303   .0012286     1.87   0.061    -.0001051    .0047111
 fp()#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5351209   .0584775     9.15   0.000      .420507    .6497348
  sd(resid.) |   .3435371   .0066191                      .3308058    .3567585
-------------+----------------------------------------------------------------
stime:       |            
         trt |   .0477638   .1797229     0.27   0.790    -.3044867    .4000143
        EV[] |   1.250361    .093301    13.40   0.000     1.067494    1.433227
       _cons |  -4.404265   .2721998   -16.18   0.000    -4.937767   -3.870763
  log(gamma) |   .0002902   .0827115     0.00   0.997    -.1618213    .1624016
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .9889317   .0422184                      .9095524    1.075239
      sd(M2) |   .1941417   .0134206                       .169542    .2223107
 corr(M2,M1) |   .4399148   .0700617                      .2930095    .5665281
------------------------------------------------------------------------------

which shows a log hazard ratio of 1.250 (95% CI: 1.067, 1.433) for a one-unit increase in log serum bilirubin at time $t$, showing very clearly the increase in the mortality rate with increasing log serum bilirubin.

Cumulative association structure

We can assess the cumulative effect of the biomarker directly to survival, for example

$$ h(t) = h_{0}(t) \exp \left( X_{1ij}\mathbf{\beta_{1}} + \alpha \int_{0}^{t} m_{i}(u) \partial u \right) $$

where the integral is approximated using numerical integration. This can be fitted with merlin using the iEV[] element. I’ll also use the estimates from the previous model as my starting values to help things along.

. merlin  (logb rcs(time,df(3)) fp(time,pow(1))#M2[id]@1 M1[id]@1                ///
>                                         , family(gaussian) timevar(time))      ///
>         (stime trt iEV[logb] , family(weibull, failure(died)) timevar(stime))  ///
>         , covariance(unstructured) restartvalues(M2 0.01)
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_3
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1

Fitting fixed effects model:

Fitting full model:

Iteration 0:   log likelihood = -3087.2721  
Iteration 1:   log likelihood = -2662.7076  (not concave)
Iteration 2:   log likelihood = -2535.5059  (not concave)
Iteration 3:   log likelihood = -2515.3547  
Iteration 4:   log likelihood = -2351.1225  (not concave)
Iteration 5:   log likelihood = -2010.1795  
Iteration 6:   log likelihood = -1983.8663  
Iteration 7:   log likelihood = -1977.6287  
Iteration 8:   log likelihood = -1977.5851  
Iteration 9:   log likelihood = -1977.5851  

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood = -1977.5851
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |            
     rcs():1 |   .1171953   .0247255     4.74   0.000     .0687343    .1656564
     rcs():2 |  -.0093604   .0041659    -2.25   0.025    -.0175254   -.0011954
     rcs():3 |   .0021762   .0012442     1.75   0.080    -.0002625    .0046148
 fp()#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .5332384   .0584448     9.12   0.000     .4186887    .6477882
  sd(resid.) |   .3453356   .0067133                      .3324253    .3587472
-------------+----------------------------------------------------------------
stime:       |            
         trt |  -.0981252   .1730614    -0.57   0.571    -.4373193    .2410689
       iEV[] |   .1582423   .0143939    10.99   0.000     .1300307    .1864538
       _cons |  -2.643807   .2009621   -13.16   0.000    -3.037686   -2.249929
  log(gamma) |  -.2265967   .1006119    -2.25   0.024    -.4237924   -.0294009
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .9872938   .0421947                      .9079626    1.073556
      sd(M2) |   .1839214   .0127884                      .1604896    .2107744
 corr(M2,M1) |   .4319423   .0737857                      .2770592    .5649472
------------------------------------------------------------------------------

Our association parameter alpha shows a log hazard ratio of 0.158 (95% CI: 0.130, 0.186), which equates to a hazard ratio of 1.171 (95% CI: 1.139, 1.204) for a one unit increase in the cumulative value of log serum bilirubin at time $t$.

Weighted cumulative association structure

We can now extend the above model by incorporating a weight function within the integral, for example to assign higher weights to more recent values of the biomarker. I’ll implement the first proposal of Mauff et al. (2017),

$$ h(t) = h_{0}(t) \exp \left( X_{1ij}\mathbf{\beta_{1}} + \alpha \int_{0}^{t} w(t-u) m_{i}(u) \partial u \right) $$

with the weight function defined as,

$$ w(t-u) = \phi((t-u)|0,\sigma^{2}) / \Phi(0,\text{max}(T)) $$

where $\phi((t-u)|0,\sigma^{2})$ is the normal density function evaluated at $t-u$, with mean 0 and variance $\sigma^{2}$, and $\Phi(0,\text{max}(T))$ the standard normal cumulative distribution function between the limits of 0 and the maximum observed survival time, $T$.

To fit such a bespoke model, we can use merlin’s family(user) functionality to define the survival submodel. Here my function will return the log hazard under a Weibull baseline hazard function. It’s actually a pretty simple function, which should hopefully be easy to follow with the comments all the way through:

. mata mata clear

. mata:
: //define my function, which needs the second input of time and we will be integrating
: //it out to get the cumulative hazard required in the likelihood                        
: real matrix userhaz(gml, t)
> {
>         //declare everything we need
>         real matrix             linpred, logh, gq, nodes, weights, intres
>         real scalar     gam1, alph, sd, maxtime, nq, Nobs
> 
>         //evaluate the main complex linear predictor log(lambda) of the Weibull model
>         linpred = merlin_util_xzb(gml,t)
>         
>         //evaluate the anciliary parameters, ensuring correct scale
>         gam1    = exp(merlin_util_ap(gml,1))            //gamma of the Weibull model
>         alph    = merlin_util_ap(gml,2)                 //association parameter
>         sd              = exp(merlin_util_ap(gml,3))    //scale parameter - standard deviation
>         sd2             = sd^2                          //variance
>         
>         //setup for numerical integration
>         nq              = 30
>         gq              = merlin_gq(nq,"legendre")      //extract basis nodes and weights
>         Nobs    = rows(t)
>         nodes   = t :* J(Nobs,1,gq[,1]'):/2 :+ t:/2
>         weights = t :* J(Nobs,1,gq[,2]'):/2
>         
>         //evaluate denominator of weight function -> constant 
>         maxtime = max(merlin_util_depvar(gml)[,1])              
>         mvnmaxt = mvnormalcv(0,maxtime,0,sd2)                   
>         
>         //conduct integration
>         //-> merlin_util_xzb_mod() evaluates the linear predictor of the biomarker, 
>         //at specified timepoints
>         intres  = J(Nobs,1,0)
>         for (q=1;q<=nq;q++) {
>                 intres = intres :+ weights[,q] :* normalden(t:-nodes[,q],sd) :* 
>                                                   merlin_util_xzb_mod(gml,1,nodes[,q]) 
>         }
>         intres = intres :/ mvnmaxt
>         
>         //build and return hazard function
>         logh = linpred :+ log(gam1) :+ (gam1:-1):*log(t)  :+ alph :* intres
>         return(logh)
> }
: end

Note that this implementation relies on the Stata 15 Mata function mvnormalcv(), so won’t work in 14.2, the minimum version required to use merlin.

Now let’s put the function into context by examining the merlin command call. The main things to note are we have 3 anciliary parameters, the first being the shape parameter of the Weibull baseline (on the log scale), the association parameter alpha, and the standard deviation of the normal density in the weight function (modelled on the log standard deviation scale). Once more we can use the starting values from the previous model fit, but this time I need to add in an extra 0 in the right place, as I have an extra parameter.

Note that specifying starting values in the right place is rather challenging unless you wrote the command…I’m working on improvements to the results table, and starting value specification.

Ok, let’s fit our model,

. mat svs = e(b)

. mat svs2 = svs[1,1..7]          //logb

. mat svs2 = svs2,svs[1,8]        //trt

. mat svs2 = svs2,svs[1,10]       //intercept

. mat svs2 = svs2,svs[1,9]        //alpha

. mat svs2 = svs2,svs[1,11]       //loggamma

. mat svs2 = svs2,-2                      //sd

. mat svs2 = svs2,svs[1,12..14]

. merlin  (logb rcs(time,df(3)) fp(time,pow(1))#M2[id]@1 M1[id]@1                           ///
>                                   , family(gaussian) timevar(time))                       ///
>         (stime trt, family(user, loghfunc(userhaz) nap(3) failure(died)) timevar(stime))  ///
>         , covariance(unstructured) from(svs2) intmethod(gh)
variables created for model 1, component 1: _cmp_1_1_1 to _cmp_1_1_3
variables created for model 1, component 2: _cmp_1_2_1 to _cmp_1_2_1

Fitting full model:

Iteration 0:   log likelihood = -2521.2863  (not concave)
Iteration 1:   log likelihood = -2286.6161  (not concave)
Iteration 2:   log likelihood =  -2283.007  (not concave)
Iteration 3:   log likelihood = -2261.6611  (not concave)
Iteration 4:   log likelihood = -2110.9218  (not concave)
Iteration 5:   log likelihood = -2092.7007  
Iteration 6:   log likelihood = -2053.8711  (not concave)
Iteration 7:   log likelihood = -2034.0065  
Iteration 8:   log likelihood = -2027.5344  
Iteration 9:   log likelihood =   -2026.74  
Iteration 10:  log likelihood = -2026.6996  
Iteration 11:  log likelihood = -2026.6996  

Mixed effects regression model                  Number of obs     =      1,945
Log likelihood = -2026.6996
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
logb:        |            
     rcs():1 |   .1289953   .0253389     5.09   0.000      .079332    .1786587
     rcs():2 |  -.0048559   .0045983    -1.06   0.291    -.0138683    .0041565
     rcs():3 |   .0010048   .0013733     0.73   0.464    -.0016869    .0036964
 fp()#M2[id] |          1          .        .       .            .           .
      M1[id] |          1          .        .       .            .           .
       _cons |   .4603155   .0256413    17.95   0.000     .4100594    .5105715
  sd(resid.) |   .4022483   .0070663                      .3886343    .4163393
-------------+----------------------------------------------------------------
stime:       |            
         trt |   .0202493   .1766731     0.11   0.909    -.3260236    .3665223
       _cons |  -4.247035   .2663483   -15.95   0.000    -4.769068   -3.725002
        ap:1 |   -.032767   .0893215    -0.37   0.714     -.207834    .1422999
        ap:2 |   1.259936   .0967404    13.02   0.000     1.070328    1.449543
        ap:3 |  -3.033368   .2071124   -14.65   0.000    -3.439301   -2.627435
-------------+----------------------------------------------------------------
id:          |            
      sd(M1) |   .8331726   .0162518                       .801921    .8656422
      sd(M2) |   .1046533   .0028968                      .0991269    .1104878
 corr(M2,M1) |   .3488695   .0464411                      .2547815     .436414
------------------------------------------------------------------------------

We obtain an estimate of alpha = 1.260 (95% CI: 1.070, 1.450), which gives us a hazard ratio of 3.525 (95% CI: 2.915, 4.263) for a one-unit increase in the weighted cumulative log serum bilirubin, at time $t$. For the standard deviation of the weight function, we get $\sigma =$ 0.048 (95% CI: 0.032, 0.072), indicating much of the weight is placed close to the timepoint $t$ of interest, with a sharp decline, indicating only the most recent values of the biomarker appear important.

There’s a lot of ways to extend this further, with more complex trajectory functions, random effects design matrices, further linking between outcomes, non-continuous longitudinal outcomes, multiple longitudinal outcomes…the list goes on. I’ll add more to this later, including the second weight function.

Eventually I will write an element type which does this weight function, such as wiEV[] which will make life a lot easier.

Related

comments powered by Disqus