# 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.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.