# Part 4 msm models with covariates

## 4.1 Reminder of covariates in multi-state models

The intensity for each \(r-s\) transition can be expressed as a log-linear function of covariates, denoted as a vector \(\mathbf{z}\).

\[q_{rs} = q_{rs}^{(0)} \exp(\boldsymbol{\beta}_{rs}^{\prime} \mathbf{z})\] Psoriatic arthritis example: investigate risk factors for joint damage

`psor`

data contains two covariates:

`hieffusn`

: Five or more effused joints`ollwsdrt`

: Low erythrocyte sedimentation rate (less than 15 mm/h)

Each row of the data contains value for that person at study entry (people just diagnosed), so covariate values are the same throughout each person’s follow-up.

`head(psor, 5)`

```
## ptnum months state hieffusn ollwsdrt
## 1 1 6.46 1 0 1
## 2 1 17.08 1 0 1
## 3 2 26.32 1 0 0
## 4 2 29.48 3 0 0
## 5 2 30.58 4 0 0
```

## 4.2 Choices for modelling covariates

There are three transitions in the model. Which transitions should be affected by which covariates?

Don’t just put all covariates on all transitions without thinking — you may not have enough data informing transitions out of particular states.

`statetable.msm(state, ptnum, data=psor)`

```
## to
## from 1 2 3 4
## 1 183 56 16 8
## 2 0 100 35 18
## 3 0 0 48 37
```

With 37 transitions, for example, unlikely to get usefully-precise effect estimates for 10 covariates on the State 3 - State 4 transition rate, but 3 covariates may be feasible.

No shrinkage / penalized likelihood methods implemented in `msm`

unfortunately, but standard variable selection methods, e.g. AIC, are simple to implement, as we will see.

## 4.3 msm syntax: covariates on rates of specific transitions

Use the `covariates`

argument to `msm`

to place covariates on transition intensities.

Set `covariates`

to a list of linear model formulae, named by the transition, to specify which intensities get which covariates.

```
<- msm(state ~ months,
psorc_spec.msm subject = ptnum,
data = psor,
qmatrix = Q,
covariates = list("1-2" = ~ hieffusn,
"2-3" = ~ hieffusn + ollwsdrt),
gen.inits=TRUE)
```

## 4.4 Covariates on all transition rates

To put the same covariates on *all* the transition rates, can set the `covariates`

argument to a single
formula (instead of a list of formulae).

```
<- msm(state ~ months,
psorc_all.msm subject = ptnum,
data = psor,
qmatrix = Q,
covariates = ~ hieffusn + ollwsdrt,
gen.inits=TRUE)
```

## 4.5 Identical effects of same covariate on different transitions

`constraint`

argument: a list with one named element for each covariate with constrained effects, so that the hazard ratio for that covariate is the same for different transitions.

Each list element is a vector of integers, of same length as the number of transition rates.

Which elements of this vector are the same determines which covariate effects are assumed to be equal.

In this model, the effect of `hieffusn`

is the same for the first two transition rates, i.e. the progression rates from State 1 to 2 and from State 2 to 3, but different for the State 3 - 4 rate

```
<- msm(state ~ months,
psorc_cons.msm subject = ptnum,
data = psor,
qmatrix = Q,
covariates = ~ hieffusn + ollwsdrt,
constraint = list(hieffusn = c(1,1,2)),
gen.inits=TRUE)
```

## 4.6 Covariate selection

Compare the three different models above using Akaike’s information criterion (AIC = -2\(\times\)log likelihood + 2*df), where “df” is the number of estimated parameters (degrees of freedom). AIC is an approximate measure of the predictive ability of a model.

There is not much difference, but the model with the lowest AIC is `psorc_spec.msm`

with different predictors on specific transitions.

```
<- c(logLik(psorc_spec.msm), logLik(psorc_cons.msm), logLik(psorc_all.msm))
logliks <- AIC(psorc_spec.msm, psorc_cons.msm, psorc_all.msm)
AICs cbind(logliks, AICs)
```

```
## logliks df AIC
## psorc_spec.msm -558 6 1128
## psorc_cons.msm -557 8 1129
## psorc_all.msm -556 9 1131
```

Most complex model `psorc_all.msm`

has highest log-likelihood, but the extra parameters do not lead to better predictive performance

It seems sufficient to assume that `hieffusn`

is a predictor of the 1-2 and 2-3 transitions, and `ollwsdrt`

is a predictor of the 2-3 transition

Check the coefficients from the fitted models to verify that this is plausible…

## 4.7 Interpretation of the hazard ratio

Hazard ratios represent how the *instantaneous risk* of making a particular transition is modified by the covariate. Here, someone in state 1, with 5 or more effused joints (`hieffusn=1`

), is at over twice the risk of progression to state 2, at any time, compared to someone with 4 or fewer (`hieffusn=0`

).

`print(psorc_all.msm, digits=2)`

```
##
## Call:
## msm(formula = state ~ months, subject = ptnum, data = psor, qmatrix = Q, gen.inits = TRUE, covariates = ~hieffusn + ollwsdrt)
##
## Maximum likelihood estimates
## Baselines are with covariates set to their means
##
## Transition intensities with hazard ratios for each covariate
## Baseline hieffusn ollwsdrt
## State 1 - State 1 -0.10 (-0.127,-0.079)
## State 1 - State 2 0.10 ( 0.079, 0.127) 2.3 (1.09,5.0) 0.73 (0.43,1.26)
## State 2 - State 2 -0.16 (-0.206,-0.128)
## State 2 - State 3 0.16 ( 0.128, 0.206) 1.7 (0.95,3.0) 0.46 (0.26,0.79)
## State 3 - State 3 -0.26 (-0.350,-0.195)
## State 3 - State 4 0.26 ( 0.195, 0.350) 1.4 (0.77,2.5) 1.58 (0.78,3.19)
##
## -2 * log-likelihood: 1113
## [Note, to obtain old print format, use "printold.msm"]
```

Can also extract with

`hazard.msm(psorc_all.msm)`

```
## $hieffusn
## HR L U
## State 1 - State 2 2.34 1.094 5.00
## State 2 - State 3 1.68 0.950 2.98
## State 3 - State 4 1.39 0.774 2.51
##
## $ollwsdrt
## HR L U
## State 1 - State 2 0.732 0.426 1.258
## State 2 - State 3 0.458 0.264 0.793
## State 3 - State 4 1.576 0.778 3.193
```

## 4.8 Outputs for different covariate values

Comparing **absolute outcomes** between groups can be easier to interpret than hazard ratios.

Most of `msm`

’s prediction functions (`totlos.msm`

, `ppass.msm`

, etc.) take an argument called `covariates`

, a list that defines the covariate values to calculate the output for.

Example: probability of transition from state 1 to severest state 4 by 10 years, compared between combinations of risk factors

`covariates`

argument indicates the covariate values \(\mathbf{z}\) to calculate the transition probabilities for.

`msm`

will calculate \(P(t | \hat{q}, \hat\beta, \mathbf{z})\) where \(\hat{q},\hat\beta\) are the fitted baseline intensities and covariate effects.

```
<- pmatrix.msm(psorc_spec.msm, t=10, covariates=list(hieffusn=0,ollwsdrt=0),
p00 ci="normal", B=100)
<- pmatrix.msm(psorc_spec.msm, t=10, covariates=list(hieffusn=0,ollwsdrt=1),
p01 ci="normal", B=100)
<- pmatrix.msm(psorc_spec.msm, t=10, covariates=list(hieffusn=1,ollwsdrt=0),
p10 ci="normal", B=100)
<- pmatrix.msm(psorc_spec.msm, t=10, covariates=list(hieffusn=1,ollwsdrt=1),
p11 ci="normal", B=100)
<- rbind(p00[1,4], p01[1,4], p10[1,4], p11[1,4]) p_14
```

Low sedimentation rate is protective (lower transition probabilities for `ollwsdrt=1`

), while high effusion is a risk factor for progression (higher transition probabilities for `hieffusn=1`

)

`cbind(hieffusn = c(0,0,1,1), ollwsdrt=c(0,1,0,1), p_14)`

```
## hieffusn ollwsdrt estimate lower upper
## [1,] 0 0 0.205 0.1500 0.270
## [2,] 0 1 0.120 0.0855 0.178
## [3,] 1 0 0.484 0.2922 0.688
## [4,] 1 1 0.324 0.1866 0.499
```

## 4.9 Exercises

We want to investigate whether heart transplant recipients who have previously had ischaemic heart disease are more at risk of CAV, which is a similar condition. The “primary diagnosis” (reason for having a heart transplant) is included in the CAV data in the categorical variable `pdiag`

, which takes the value `IHD`

for ischaemic heart disease.

Convert this to a binary variable, and develop a multi-state model like the one you built in the previous chapters, including this binary covariate on all transitions.

Adapt the model appropriately so that the covariate is only included on transitions which the covariate is predictive of.

Investigate and compare with a model where the effect of this covariate is constrained to be the same between different transitions.

Compare the clinical prognosis for people with and without pre-transplant IHD, using outcomes similar to those calculated in the previous chapter.

(Note: if calculating the observed and expected prevalences for a subset of patients, then a

`subset`

argument to`plot.prevalence.msm`

should be specified. This should be set to a vector of the unique patient IDs included in the appropriate subset.)

## 4.10 Solutions

We convert the categorical variable

`pdiag`

to a binary 0/1 covariate indicating primary diagnosis of IHD. First we fit the full model where the binary covariate affects all transitions, and examine the hazard ratios under this model.`$IHD <- as.numeric(cav$pdiag == "IHD") cav<- msm(state ~ years, subject=PTNUM, data=cav, covariates = ~IHD, cavcov.msm deathexact = TRUE, gen.inits = TRUE, qmatrix=Q_cav) cavcov.msm`

`## ## Call: ## msm(formula = state ~ years, subject = PTNUM, data = cav, qmatrix = Q_cav, gen.inits = TRUE, covariates = ~IHD, deathexact = TRUE) ## ## Maximum likelihood estimates ## Baselines are with covariates set to their means ## ## Transition intensities with hazard ratios for each covariate ## Baseline IHD ## State 1 - State 1 -0.17057 (-0.19084,-0.15246) ## State 1 - State 2 0.12845 ( 0.11159, 0.14785) 1.5279 (1.15331, 2.024) ## State 1 - State 4 0.04212 ( 0.03365, 0.05274) 1.4870 (0.94878, 2.331) ## State 2 - State 1 0.22974 ( 0.17089, 0.30885) 0.9270 (0.51291, 1.675) ## State 2 - State 2 -0.61327 (-0.71676,-0.52472) ## State 2 - State 3 0.34288 ( 0.27220, 0.43191) 0.9633 (0.60716, 1.528) ## State 2 - State 4 0.04065 ( 0.01130, 0.14623) 0.9333 (0.07213,12.076) ## State 3 - State 2 0.13300 ( 0.08069, 0.21922) 0.7380 (0.27167, 2.005) ## State 3 - State 3 -0.44500 (-0.56314,-0.35165) ## State 3 - State 4 0.31200 ( 0.24218, 0.40196) 0.7412 (0.44660, 1.230) ## ## -2 * log-likelihood: 3925 ## [Note, to obtain old print format, use "printold.msm"]`

We see that the pre-transplant IHD is associated with an increased risk of CAV onset (state 1 - state 2), and possibly also with death before CAV onset (state 1 - state 4). The confidence intervals for the hazard ratios of IHD on the other transition rates are wide and include the null 1.

Hence we build a simplified model, in which IHD only affects the transition rates out of State 1. This model has an improved AIC compared to the bigger model, despite having 5 fewer parameters (

`df`

).`<- msm(state ~ years, subject=PTNUM, data=cav, cavcov1.msm covariates = list("1-2"=~IHD, "1-4"=~IHD), deathexact = TRUE, gen.inits = TRUE, qmatrix=Q_cav) AIC(cavcov.msm, cavcov1.msm)`

`## df AIC ## cavcov.msm 14 3953 ## cavcov1.msm 9 3945`

The hazard ratios for IHD on the risk of CAV onset and pre-CAV death appear to be similar, so we might constrain them to be equal, as follows. The AIC is slightly smaller than for the model where these effects are different.

`<- msm(state ~ years, subject=PTNUM, data=cav, cavcov2.msm covariates = list("1-2"=~IHD, "1-4"=~IHD), constraint = list(IHD = c(1,1)), deathexact = TRUE, gen.inits = TRUE, qmatrix=Q_cav) AIC(cavcov2.msm)`

`## [1] 3943`

However despite the statistical improvement, this model would not be clinically meaningful if pre-transplant IHD affects the risk of CAV onset and death after CAV through different mechanisms. We would then prefer to acknowledge that the true effects might be different.

We could compare the transition probabilities at a future time (e.g. 5 years after transplant) for people with and without pre-transplant IHD. The risk of occupying any of the CAV states at this time, and the risk of death by this time, is about 50% more for people with previous IHD, while the chance of being CAV-free (state 1) is correspondingly lower.

`pmatrix.msm(cavcov.msm, t=5, covariates=list(IHD=0))`

`## State 1 State 2 State 3 State 4 ## State 1 0.5884 0.1197 0.0703 0.222 ## State 2 0.2747 0.1273 0.1459 0.452 ## State 3 0.0716 0.0647 0.1264 0.737 ## State 4 0.0000 0.0000 0.0000 1.000`

`pmatrix.msm(cavcov.msm, t=5, covariates=list(IHD=1))`

`## State 1 State 2 State 3 State 4 ## State 1 0.4494 0.1571 0.114 0.280 ## State 2 0.2189 0.1483 0.210 0.423 ## State 3 0.0537 0.0712 0.207 0.668 ## State 4 0.0000 0.0000 0.000 1.000`

The passage probability to severe CAV (which we computed in the previous chapter), i.e. the chance of going through the “severe CAV” state before death, is similar between both groups of people (

`IHD=0`

) and (`IHD=1`

), though slightly higher for people who have had IHD.`ppass.msm(cavcov1.msm, tot=100, covariates=list(IHD=0))["State 1","State 3"]`

`## [1] 0.574`

`ppass.msm(cavcov1.msm, tot=100, covariates=list(IHD=1))["State 1","State 3"]`

`## [1] 0.595`

We might also compare the total length of stay in each state over a long time horizon (e.g. 20 years), finding that people who had IHD are expected to spend longer periods of time in the CAV states (2 and 3) and die sooner.

`totlos.msm(cavcov1.msm, tot=20, covariates=list(IHD=0))`

`## State 1 State 2 State 3 State 4 ## 8.80 1.66 1.21 8.33`

`totlos.msm(cavcov1.msm, tot=20, covariates=list(IHD=1))`

`## State 1 State 2 State 3 State 4 ## 6.47 1.96 1.45 10.11`

The expected times from state 1 to death (

`tostate=4`

) in the two groups can be compared as follows. These are shorter for people who have had IHD.`efpt.msm(cavcov1.msm, tostate=4, covariates=list(IHD=0), ci="normal")[,1]`

`## est 2.5% 97.5% ## 14.8 12.5 16.9`

`efpt.msm(cavcov1.msm, tostate=4, covariates=list(IHD=1), ci="normal")[,1]`

`## est 2.5% 97.5% ## 11.26 9.74 12.60`

Let’s also compare the observed and expected prevalences. Note the use of the

`subset`

argument to restrict the observed prevalences to those observed in the indicated set of patient IDs.`<- unique(cav$PTNUM[cav$IHD==0]) notihd_pts <- unique(cav$PTNUM[cav$IHD==1]) ihd_pts plot.prevalence.msm(cavcov.msm, covariates=list(IHD=0), subset=notihd_pts)`

`plot.prevalence.msm(cavcov.msm, covariates=list(IHD=1), subset=ihd_pts)`

If desired, the function

`prevalence.msm`

could be used to extract the numbers behind the plots shown by`plot.prevalence.msm`

, to enable multiple groups to be shown on the same plot, or to produce other customised plots.