# 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

. 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.449
Iteration 2:   log likelihood =  -1938.343
Iteration 3:   log likelihood = -1918.0422
Iteration 4:   log likelihood = -1909.4608
Iteration 5:   log likelihood = -1909.0986
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.) |   .3435372   .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.870764
log(gamma) |   .0002902   .0827114     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.219 (95% CI: 1.041, 1.398) for a one-unit increase in log serum bilirubin at time $t$, which equates to a hazard ratio of 3.385 (95% CI: 2.832, 4.047) 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.8943  (not concave)
Iteration 2:   log likelihood = -2535.6354  (not concave)
Iteration 3:   log likelihood = -2514.7156
Iteration 4:   log likelihood = -2428.0628  (not concave)
Iteration 5:   log likelihood = -2029.9276
Iteration 6:   log likelihood = -1980.4806
Iteration 7:   log likelihood = -1977.5937
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    .1656563
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[] |   .1582421   .0143939    10.99   0.000     .1300305    .1864536
_cons |  -2.643812   .2009622   -13.16   0.000     -3.03769   -2.249933
log(gamma) |   -.226594   .1006117    -2.25   0.024    -.4237893   -.0293987
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   .9872938   .0421947                      .9079627    1.073556
sd(M2) |   .1839214   .0127884                      .1604896    .2107744
corr(M2,M1) |   .4319422   .0737857                      .2770592    .5649472
------------------------------------------------------------------------------



Our association parameter alpha shows a log hazard ratio of 0.160 (95% CI: 0.130, 0.189), which equates to a hazard ratio of 1.174 (95% CI: 1.139, 1.208) 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:
------------------------------------------------- mata (type end to exit) ---------------
: //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, alpha, 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 - sta
> ndard 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)
> }
note: variable alpha unused

: 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.2856  (not concave)
Iteration 1:   log likelihood = -2286.6158  (not concave)
Iteration 2:   log likelihood = -2282.9909  (not concave)
Iteration 3:   log likelihood = -2261.5794  (not concave)
Iteration 4:   log likelihood = -2111.0234  (not concave)
Iteration 5:   log likelihood = -2092.8191
Iteration 6:   log likelihood = -2050.0109  (not concave)
Iteration 7:   log likelihood =  -2031.439
Iteration 8:   log likelihood = -2029.0544
Iteration 9:   log likelihood = -2028.0068
Iteration 10:  log likelihood = -2026.7388
Iteration 11:  log likelihood = -2026.6998
Iteration 12:  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 |   .1289967   .0253711     5.08   0.000     .0792701    .1787232
rcs():2 |  -.0048557   .0046057    -1.05   0.292    -.0138826    .0041713
rcs():3 |   .0010047   .0013755     0.73   0.465    -.0016913    .0037007
fp()#M2[id] |          1          .        .       .            .           .
M1[id] |          1          .        .       .            .           .
_cons |   .4603147   .0256491    17.95   0.000     .4100434    .5105859
sd(resid.) |   .4022483   .0070663                      .3886343    .4163393
-------------+----------------------------------------------------------------
stime:       |
trt |   .0202497   .1766732     0.11   0.909    -.3260234    .3665227
_cons |   -4.24704   .2663488   -15.95   0.000    -4.769074   -3.725006
ap:1 |  -.0327652   .0893214    -0.37   0.714    -.2078318    .1423015
ap:2 |   1.259937   .0967406    13.02   0.000     1.070329    1.449545
ap:3 |  -3.033374   .2071093   -14.65   0.000    -3.439301   -2.627447
-------------+----------------------------------------------------------------
id:          |
sd(M1) |   .8331727   .0162519                      .8019208    .8656426
sd(M2) |   .1046533   .0028968                      .0991269    .1104878
corr(M2,M1) |   .3488691   .0464445                       .254774    .4364196
------------------------------------------------------------------------------



We obtain an estimate of alpha = 1.236 (95% CI: 1.055, 1.417), which gives us a hazard ratio of 3.443 (95% CI: 2.873, 4.126) 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.050 (95% CI: 0.031, 0.078), 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.

This example used merlin version 2.0.2.