added age adjusted msm

Snakefile

@@ -120,7 +120,7 @@ rule fit_model_validation_set:
 rule model_posteriors:
     input:
         ["output/v1.1/data/validation/posteriors/%s/df_posterior_mi_%i_fold_%i.rds" % (m, i, j)
-          for m in ("locf", "msm", "gp", "gp_nb", "mm", "mm_nb")
+          for m in ("locf", "msm", "msm_age", "gp", "gp_nb", "mm", "mm_nb")
           for i in range(1, config["mi_m"] + 1)
           for j in range(1, config["folds"] + 1)
         ]

models/msm_age/fit.R

+#!/usr/bin/env Rscript
+library(tidyverse)
+library(msm)
+
+args <- commandArgs(trailingOnly = TRUE)
+
+inputfile  <- args[[1]]
+outputfile <- args[[2]]
+
+config <- yaml::read_yaml("config.yml")
+
+# filter deaths - we use exact times
+# GOSE = 2 is too infreqeunt
+df_gose <- readRDS(inputfile) %>%
+  mutate(
+    Outcomes.DerivedCompositeGOSE = factor(
+      Outcomes.DerivedCompositeGOSE,
+      levels  = c(1, 3:8, 99), # 2 is too infrequent to be fitted
+      ordered = TRUE
+    ) %>% as.character %>% as.integer
+  ) %>% select(
+    gupi,
+    Outcomes.DerivedCompositeGOSE,
+    Outcomes.DerivedCompositeGOSEDaysPostInjury
+  ) %>%
+  arrange(
+    gupi,
+    Outcomes.DerivedCompositeGOSEDaysPostInjury
+  )
+
+# translate to msm formatting
+tmp <- df_gose %>%
+  rbind(tibble( # only predict GOSE at 180 days - takes too long otherwise
+    gupi   = df_gose$gupi %>% unique,
+    Outcomes.DerivedCompositeGOSEDaysPostInjury = config$t_out_msm + .5, # needed to offset
+    Outcomes.DerivedCompositeGOSE = 99
+  )) %>%
+  arrange(gupi, Outcomes.DerivedCompositeGOSEDaysPostInjury) %>%
+  mutate(
+    obstype = ifelse(Outcomes.DerivedCompositeGOSE == 1, 3, 1), # make death exactly observed transition
+    Outcomes.DerivedCompositeGOSE = factor(
+      Outcomes.DerivedCompositeGOSE,
+      levels = c(1, 3:8, 99),
+      labels = c(1:7, 99)
+    ) %>% as.character %>% as.numeric
+  ) %>%
+  group_by(gupi, Outcomes.DerivedCompositeGOSEDaysPostInjury) %>%
+  arrange(gupi, Outcomes.DerivedCompositeGOSEDaysPostInjury) %>%
+  # remove prediction states where we already have a true one
+  filter((n() == 1) | (Outcomes.DerivedCompositeGOSE != 99)) %>%
+  ungroup
+
+# define transition matrix
+Q <- matrix(0, nrow = 7, ncol = 7)
+for (i in 1:7) {
+  for (j in 1:7) {
+    if (i == j - 1 & i != 1) {
+      Q[i, j] <- 1
+    }
+    if (i == j + 1) {
+      Q[i, j] <- 1
+    }
+  }
+}
+Q[2:7, 1] <- 1 # allow instantaneous deaths
+
+fit <- msm(
+  Outcomes.DerivedCompositeGOSE ~ Outcomes.DerivedCompositeGOSEDaysPostInjury,
+  subject   = tmp$gupi,
+  data      = tmp,
+  obstype   = tmp$obstype,
+  gen.inits = TRUE,
+  qmatrix   = Q,
+  censor    = 99,
+  pci       = c(90, 270),
+  covariates = ~ Subject.Age,
+  censor.states = 1:7,
+  control   = list(
+    fnscale = config$fnscale_msm,
+    maxit   = config$maxiter_msm,
+    trace   = 2
+  )
+)
+
+# get fitted values at censored (prediction) states
+fitted <- viterbi.msm(fit)
+
+# store posteriors
+df_posteriors <- tibble(
+    gupi = fitted$subject,
+    t    = fitted$time
+  ) %>%
+  cbind(
+    fitted$pstate %>% {colnames(.) <- c("1", as.character(3:8)); .}
+  ) %>%
+  mutate("2" = 0) %>%
+  filter(t == config$t_out_msm + .5) %>%
+  mutate(t = config$t_out_msm) %>%
+  gather(GOSE, probability, as.character(1:8)) %>%
+  arrange(gupi, t, GOSE)
+
+saveRDS(df_posteriors, outputfile)