```{rdiagram,echo=FALSE,fig.cap="Structural diagram of allowed transitions between GOSe states for the proposed multi-state model.",fig.width=8,fig.height=4}
Mixed effects models are a widely used approach in longitudinal
data analysis and model individual deviations from a population mean trajectory
(@verbeke2009linear).
To account for the fact that the GOSe outcome is an ordered factor,
we employ a cumulative link function model with flexible intercepts
(@Agresti2003).
The population mean is modeled as cubic spline function to allow a
non-linear population mean trajectory.
Patient-individual deviations from this population mean are modeled
as quadratic polynomials to allow sufficient flexibility (random effects).
Baseline covariates are added as linear fixed effects to to the
population mean.
The model was fitted using Bayesian statistics via the BRMS
package (@brms2017; @brms2018) for the R environment for statistical
computing (@R2016) and the Stan modelling language (@stan2017).
Non-informative priors were used for the model parameters.
A potential drawback of the proposed longitudinal mixed effects model
is the fact that the individual deviations from the population mean
are modeled globally using polynomials.
Since linear and quadratic terms are not identifiable for patients with
only one observed GOSE value, this implies large uncertainty over the
patient-specific effects for individuals with one or two observations
only by falling back on the non-informative priors on these model
parameters.
Thus, the overall uncertainty associated with model-based imputations
at exactly 180 days may become relatively large.
### Gaussian process model
Gaussian process regression allows flexible modelling of
both the individual GOSE trajectories as well as the population mean
in a Bayesian non-parametric way [@rasmussen2006].
This non-parametric paradigm leads to low model-uncertainty in the
vicinity of actually observed GOSe outcomes.
To account for the discreteness of the GOSe outcome,
the continuous output of the Gaussian process model is rounded to the
nearest integer in 1 to 8 (GOSe categories).
The squared exponential covariance function with shared
length scale for all individuals is used to model intra-individual
dependency of GOSe outcomes.
The population mean trajectory of the Gaussian process is modeled as
mean of an independent Gaussian process with pseudo observations at
`r config$gp$t_knots %>% as.numeric()` days post-injury
(also using a squared exponential covariance function).
This approach maintains flexibility of the population mean function
similar to the spline-based approach for the mixed effects model
while avoiding the the computational complexity of a fully
hierarchical Gaussian process model.
Again, the impact of baseline covariates is modeled via linear effects
on the population mean of the Gaussian process.
All parameters are estimated in a fully Bayesian fashion using the Stan
modelling language @stan2017 and non-informative priors except for the
length scale of the squared exponential kernel.
Due to the sparseness of the data, the estimated length scale will naturally
tend towards extremely large values implying unrealistically
long-range dependency between observations.
We therefore limit the length scale to a maximum of 120 days (4 months) and
impose a Gaussian prior with a mean of 60 days post injury and and a
standard deviation of 14.
### Multi-state model
Both the mixed effects model as well as the Gaussian process regression
model are essentially non-linear regression techniques for longitudinal
data.
While they are both powerful tools to model longitudinal trajectories,
they do not explicitly model the probability of transitions between
GOSe states.
Since the number of observations per individual is very limited in our
data set (1 to 4 GOSe observations per individual),
an approach explicitly modelling transition probabilities might be
more suitable to capture the dynamics of the GOSe trajectories.
To explore this further, a Markov multi-state model is considered (@meira2009).
This model class assumes that the transitions between adjacent GOSe
states can be modeled as a Markov process and the transition
intensities between adjacent states are fitted to the observed data.
```{r diagram, echo=FALSE, fig.cap="Structural diagram of allowed transitions between GOSe states for the proposed multi-state model.", fig.width=8, fig.height=4}
# define transition matrix
Q <- matrix(0, nrow = 8, ncol = 8)
for (i in 1:8) {
for (j in 1:8) {
if (i == j - 1 & i != 1) {
Q[i, j] <- 1
}
if (i == j + 1) {
Q[i, j] <- 1
}
}
}
Q[2:8, 1] <- 1 # allow instantaneous deaths
Q <- t(Q)
curvature <- Q / 5
curvature[1, ] <- 0
diagram::plotmat(
Q,
relsize = 1,
shadow.size = 0,
pos = c(1, 7),
box.size = .033,
arr.pos = .55,
arr.length = .33,
arr.type = "triangle",
cex.txt = 0,
curve = curvature
)
```
To account for the fact that state-transitions might be more frequent
in the early post-injury phase,
piece wise constant transition intensities were fitted to the intervals
[0, 90), [90, 270), and 270+ days post-injury.
The model was fit using the msm @msm2011 package for the R environment
for statistical computing @R2016.
Due to the relatively large number of 19 transition intensities in the proposed model
(cf. arrows in Figure ???, structure of transition graph),
inclusion of all baseline covariates turned out to be numerically unstable.
For the MSM model, instead of including all covariates,
only a model adjusting for age at injury via a proportional hazard approach was fit.
