Initial commit

time-dependent-9-tunnel-states.R

@@ -0,0 +1,67 @@
+# We will use cycles from 0 to 100
+cycles <- 0:100
+
+# We use 10 states for Diseased (9 tunnel and 1 with self-loop)
+n_tunnel <- 9
+states <- c("Healthy", paste0("Diseased", sprintf("%02i", 1:(n_tunnel+1))), "Dead")
+
+# Create an array for the state membership
+state_membership <- array(NA_real_, c(length(cycles), length(states)), list(cycle = cycles, state = states))
+state_membership[1, ] <- c(1, rep(0, length(states) - 1))
+
+# Define the survival model parameters
+# - Healthy to Diseased
+lambda_inc  <- 2e-5
+gamma_inc   <- 2.5
+# - Healthy to Dead
+a_general   <- 0.2
+b_general   <- 1e-7
+# - Diseased to Dead
+lambda_surv <- 0.2
+gamma_surv  <- 0.5
+
+# Calculate exponential distribution with same mean survival as conditional
+# Weibull distribution for Diseased to Dead given survived to 10th state
+upper_incomplete_gamma <- function(a, x) pgamma(x, a, lower = FALSE) * gamma(a)
+mean_surv <- exp(lambda_surv * n_tunnel ^ gamma_surv) * lambda_surv ^ (-1 / gamma_surv) * upper_incomplete_gamma(1 + 1 / gamma_surv, lambda_surv * n_tunnel ^ gamma_surv)
+rate_surv <- 1 / mean_surv
+
+# Calculate the hazard rates for the transitions out of the Healthy state
+healthy_haz <- array(NA_real_, c(length(cycles), 2), list(cycle = cycles, to = c("Diseased", "Dead")))
+healthy_haz[, 1] <- lambda_inc * ((cycles + 1) ^ gamma_inc - cycles ^ gamma_inc)
+healthy_haz[, 2] <- b_general / a_general * (exp(a_general * (cycles + 1)) - exp(a_general * cycles))
+
+# Create an array for transition probabilities in each cycle
+transition_probs <- array(
+  NA_real_,
+  c(length(cycles), 2 + n_tunnel + 1),
+  list(
+    cycle = cycles,
+    transition = c(
+      "Healthy->Diseased01",
+      "Healthy->Dead",
+      paste0("Diseased", sprintf("%02i", 1:(n_tunnel)), "->Diseased", sprintf("%02i", 2:(n_tunnel + 1))),
+      paste0("Diseased", sprintf("%02i", n_tunnel + 1), "->Diseased", sprintf("%02i", n_tunnel + 1)))))
+transition_probs[, 1:2] <- healthy_haz / rowSums(healthy_haz) * (1 - exp(-rowSums(healthy_haz)))
+for (t in 1:n_tunnel) {
+  transition_probs[, t + 2] <- exp(-lambda_surv * (t ^ gamma_surv - (t - 1) ^ gamma_surv))
+}
+transition_probs[, n_tunnel + 3] <- exp(-rate_surv)
+
+# Run cohort simulation
+for (i in 2:length(cycles)) {
+  # Update state membership
+  state_membership[i, "Healthy"] <- state_membership[i - 1, "Healthy"] * (1 - transition_probs[i - 1, "Healthy->Diseased01"] - transition_probs[i - 1, "Healthy->Dead"])
+  state_membership[i, "Diseased01"] <- state_membership[i - 1, "Healthy"] * transition_probs[i - 1, "Healthy->Diseased01"]
+  for (t in 1:n_tunnel) {
+    from <- paste0("Diseased", sprintf("%02i", t))
+    to   <- paste0("Diseased", sprintf("%02i", t + 1))
+    state_membership[i, to] <- state_membership[i - 1, from] * transition_probs[i - 1, paste0(from, "->", to)]
+  }
+  last <- paste0("Diseased", sprintf("%02i", n_tunnel + 1))
+  state_membership[i, last] <- state_membership[i - 1, last] * transition_probs[i - 1, paste0(last, "->", last)] + state_membership[i, last]
+  state_membership[i, "Dead"] <- 1 - sum(state_membership[i, 1:(n_tunnel+2)])
+}
+
+summarised_state_membership <- cbind(state_membership[, "Healthy"], rowSums(state_membership[, 2:(n_tunnel + 2)]), state_membership[, "Dead"])
+dimnames(summarised_state_membership) <- list(cycle = cycles, state = c("Healthy", "Diseased", "Dead"))

time-dependent-98-tunnel-states.R

@@ -0,0 +1,67 @@
+# We will use cycles from 0 to 100
+cycles <- 0:100
+
+# We use 10 states for Diseased (9 tunnel and 1 with self-loop)
+n_tunnel <- 98
+states <- c("Healthy", paste0("Diseased", sprintf("%02i", 1:(n_tunnel+1))), "Dead")
+
+# Create an array for the state membership
+state_membership <- array(NA_real_, c(length(cycles), length(states)), list(cycle = cycles, state = states))
+state_membership[1, ] <- c(1, rep(0, length(states) - 1))
+
+# Define the survival model parameters
+# - Healthy to Diseased
+lambda_inc  <- 2e-5
+gamma_inc   <- 2.5
+# - Healthy to Dead
+a_general   <- 0.2
+b_general   <- 1e-7
+# - Diseased to Dead
+lambda_surv <- 0.2
+gamma_surv  <- 0.5
+
+# Calculate exponential distribution with same mean survival as conditional
+# Weibull distribution for Diseased to Dead given survived to 10th state
+upper_incomplete_gamma <- function(a, x) pgamma(x, a, lower = FALSE) * gamma(a)
+mean_surv <- exp(lambda_surv * n_tunnel ^ gamma_surv) * lambda_surv ^ (-1 / gamma_surv) * upper_incomplete_gamma(1 + 1 / gamma_surv, lambda_surv * n_tunnel ^ gamma_surv)
+rate_surv <- 1 / mean_surv
+
+# Calculate the hazard rates for the transitions out of the Healthy state
+healthy_haz <- array(NA_real_, c(length(cycles), 2), list(cycle = cycles, to = c("Diseased", "Dead")))
+healthy_haz[, 1] <- lambda_inc * ((cycles + 1) ^ gamma_inc - cycles ^ gamma_inc)
+healthy_haz[, 2] <- b_general / a_general * (exp(a_general * (cycles + 1)) - exp(a_general * cycles))
+
+# Create an array for transition probabilities in each cycle
+transition_probs <- array(
+  NA_real_,
+  c(length(cycles), 2 + n_tunnel + 1),
+  list(
+    cycle = cycles,
+    transition = c(
+      "Healthy->Diseased01",
+      "Healthy->Dead",
+      paste0("Diseased", sprintf("%02i", 1:(n_tunnel)), "->Diseased", sprintf("%02i", 2:(n_tunnel + 1))),
+      paste0("Diseased", sprintf("%02i", n_tunnel + 1), "->Diseased", sprintf("%02i", n_tunnel + 1)))))
+transition_probs[, 1:2] <- healthy_haz / rowSums(healthy_haz) * (1 - exp(-rowSums(healthy_haz)))
+for (t in 1:n_tunnel) {
+  transition_probs[, t + 2] <- exp(-lambda_surv * (t ^ gamma_surv - (t - 1) ^ gamma_surv))
+}
+transition_probs[, n_tunnel + 3] <- exp(-rate_surv)
+
+# Run cohort simulation
+for (i in 2:length(cycles)) {
+  # Update state membership
+  state_membership[i, "Healthy"] <- state_membership[i - 1, "Healthy"] * (1 - transition_probs[i - 1, "Healthy->Diseased01"] - transition_probs[i - 1, "Healthy->Dead"])
+  state_membership[i, "Diseased01"] <- state_membership[i - 1, "Healthy"] * transition_probs[i - 1, "Healthy->Diseased01"]
+  for (t in 1:n_tunnel) {
+    from <- paste0("Diseased", sprintf("%02i", t))
+    to   <- paste0("Diseased", sprintf("%02i", t + 1))
+    state_membership[i, to] <- state_membership[i - 1, from] * transition_probs[i - 1, paste0(from, "->", to)]
+  }
+  last <- paste0("Diseased", sprintf("%02i", n_tunnel + 1))
+  state_membership[i, last] <- state_membership[i - 1, last] * transition_probs[i - 1, paste0(last, "->", last)] + state_membership[i, last]
+  state_membership[i, "Dead"] <- 1 - sum(state_membership[i, 1:(n_tunnel+2)])
+}
+
+summarised_state_membership <- cbind(state_membership[, "Healthy"], rowSums(state_membership[, 2:(n_tunnel + 2)]), state_membership[, "Dead"])
+dimnames(summarised_state_membership) <- list(cycle = cycles, state = c("Healthy", "Diseased", "Dead"))

time-dependent-CSSTD.R

@@ -0,0 +1,109 @@
+# We will use cycles from 0 to 100
+cycles <- 0:100
+
+# Create an array for the state membership
+state_membership <- array(NA_real_, c(length(cycles), 3), list(cycle = cycles, state = c("Healthy", "Diseased", "Dead")))
+state_membership[1, ] <- c(1, 0, 0)
+
+# Define the survival model parameters
+# - Healthy to Diseased
+lambda_inc  <- 2e-5
+gamma_inc   <- 2.5
+# - Healthy to Dead
+a_general   <- 0.2
+b_general   <- 1e-7
+# - Diseased to Dead
+lambda_surv <- 0.2
+gamma_surv  <- 0.5
+
+# Calculate the hazard rates for the transitions out of the Healthy state
+healthy_haz <- array(NA_real_, c(length(cycles), 2), list(cycle = cycles, to = c("Diseased", "Dead")))
+healthy_haz[, 1] <- lambda_inc * ((cycles + 1) ^ gamma_inc - cycles ^ gamma_inc)
+healthy_haz[, 2] <- b_general / a_general * (exp(a_general * (cycles + 1)) - exp(a_general * cycles))
+
+# Create an array for transition probabilities in each cycle
+transition_probs <- array(NA_real_, c(length(cycles), 3), list(cycle = cycles, transition = c("Healthy->Diseased", "Healthy->Dead", "Diseased->Dead")))
+
+# The transition probabilities for the Healthy state can be calculated already
+transition_probs[, 1:2] <- healthy_haz / rowSums(healthy_haz) * (1 - exp(-rowSums(healthy_haz)))
+
+# Create an array for the moments of the overall CSSTD in each cycle
+csstd_moments <- array(NA_real_, c(length(cycles), 5), list(cycle = cycles, moment = c("1st raw (mean)", "2nd raw", "3rd raw", "2nd central (variance)", "3rd central")))
+
+# Create an array for the weights of the two components in the CSSTD
+component_weights <- array(NA_real_, c(length(cycles), 2), list(cycle = cycles, component = c("Remaining in state", "Entering state")))
+
+# Create an array for the moments for the first component (remaining in the state)
+remaining_moments <- array(NA_real_, c(length(cycles), 6), list(cycle = cycles, moment = c("1st raw (mean)", "1st around 1", "2nd around 1", "3rd around 1", "2nd central (variance)", "3rd central")))
+
+# Create a vector for the probability for remaining in the Diseased state
+ephi <- rep(NA_real_, length(cycles))
+
+# Create an array for rho
+rho <- array(NA_real_, c(length(cycles), 3), list(cycle = cycles, derivative = 0:2))
+
+# Create an expression for phi(s), the conditional probability of staying in the state one more cycle
+# Note: We aren't evaluating this yet - it is just symbolic
+phi <- expression( exp(-lambda_surv * ((s + 1) ^ gamma_surv - s ^ gamma_surv)) )
+
+# Create expressions for the first two derivatives of phi(s)
+dphi  <- D(  phi, "s")
+d2phi <- D( dphi, "s")
+d3phi <- D(d2phi, "s")
+
+# Initialise the model
+csstd_moments[1, ] <- 0
+ephi[1] <- eval(phi, list(s = 0))
+transition_probs[1, 3] <- 1 - ephi[1]
+rho[1, 1] <- 1
+rho[1, 2:3] <- 0
+
+# Iterate through the cycles
+for (i in 2:length(cycles)) {
+  
+  # Update state membership
+  state_membership[i, 1] <- state_membership[i - 1, 1] * (1 - sum(transition_probs[i - 1, 1:2]))
+  state_membership[i, 2] <- state_membership[i - 1, 1] * transition_probs[i - 1, 1] +
+    state_membership[i - 1, 2] * (1 - transition_probs[i - 1, 3])
+  state_membership[i, 3] <- state_membership[i - 1, 1] * transition_probs[i - 1, 2] +
+    state_membership[i - 1, 2] * transition_probs[i - 1, 3] +
+    state_membership[i - 1, 3]
+  
+  # Calculate weights
+  component_weights[i, 2] <- state_membership[i - 1, 1] * transition_probs[i - 1, 1] / state_membership[i, 2]
+  component_weights[i, 1] <- 1 - component_weights[i, 2]
+  
+  # Calculate moments for those remaining in the state (since the previous cycle)
+  remaining_moments[i, 1] <- 1 + csstd_moments[i - 1, "1st raw (mean)"] * rho[i - 1, 1] +
+    (csstd_moments[i - 1, "2nd raw"] - csstd_moments[i - 1, "1st raw (mean)"] ^ 2) * rho[i - 1, 2] +
+    (csstd_moments[i - 1, "3rd raw"] - 2 * csstd_moments[i - 1, "1st raw (mean)"] * csstd_moments[i - 1, "2nd raw"] + csstd_moments[i - 1, "1st raw (mean)"] ^ 3) * rho[i - 1, 3] / 2
+  remaining_moments[i, 2] <- remaining_moments[i, 1] - 1
+  remaining_moments[i, 3] <- csstd_moments[i - 1, "2nd raw"] * rho[i - 1, 1] +
+    (csstd_moments[i - 1, "3rd raw"] - csstd_moments[i - 1, "1st raw (mean)"] * csstd_moments[i - 1, "2nd raw"]) * rho[i - 1, 2]
+  remaining_moments[i, 4] <- csstd_moments[i - 1, "3rd raw"] * rho[i - 1, 1]
+  remaining_moments[i, 5] <- remaining_moments[i, 3] - (1 - remaining_moments[i, 1]) ^ 2
+  remaining_moments[i, 6] <- remaining_moments[i, 4] + 3 * remaining_moments[i, 3] * (1 - remaining_moments[i, 1]) - 2 * (1 - remaining_moments[i, 1]) ^ 3
+  
+  # Calculate overall CSSTD moments
+  csstd_moments[i, 1] <- component_weights[i, 1] * remaining_moments[i, 1] + component_weights[i, 2] * 0.5
+  temp <- remaining_moments[i, 1] - csstd_moments[i, 1]
+  csstd_moments[i, 4] <- component_weights[i, 1] * (temp^2 + remaining_moments[i, 5]) + component_weights[i, 2] * (0.5 - csstd_moments[i, 1]) ^ 2
+  csstd_moments[i, 5] <- component_weights[i, 1] * (temp^3 + 3 * remaining_moments[i, 5] * temp + remaining_moments[i, 6]) + component_weights[i, 2] * (0.5 - csstd_moments[i, 1])^3
+  csstd_moments[i, 2] <- csstd_moments[i, 1] ^ 2 + csstd_moments[i, 4]
+  csstd_moments[i, 3] <- csstd_moments[i, 1] ^ 3 + 3 * csstd_moments[i, 4] * csstd_moments[i, 1] + csstd_moments[i, 5]
+  
+  # Calculate E[phi(s)]
+  mu <- csstd_moments[i, 1]
+  phi_eval   <- eval(  phi, list(s = mu))
+  dphi_eval  <- eval( dphi, list(s = mu))
+  d2phi_eval <- eval(d2phi, list(s = mu))
+  d3phi_eval <- eval(d3phi, list(s = mu))
+  
+  ephi[i] <- phi_eval + d2phi_eval / 2 * csstd_moments[i, "2nd central (variance)"] + d3phi_eval / 6 * csstd_moments[i, "3rd central"]
+  transition_probs[i, 3] <- 1 - ephi[i]
+  
+  # Calculate rho and derivatives
+  rho[i, 1] <- phi_eval / ephi[i]
+  rho[i, 2] <- dphi_eval / ephi[i]
+  rho[i, 3] <- d2phi_eval / ephi[i]
+}

time-dependent-comparison.R

@@ -0,0 +1,45 @@
+library(tidyverse)
+
+csstd      <- new.env()
+no_tunnel  <- new.env()
+tunnel9    <- new.env()
+tunnel98   <- new.env()
+phase_type <- new.env()
+
+system.time(source("time-dependent-CSSTD.R",            echo = TRUE, local = csstd))
+system.time(source("time-dependent-no-tunnel-states.R", echo = TRUE, local = no_tunnel))
+system.time(source("time-dependent-9-tunnel-states.R",  echo = TRUE, local = tunnel9))
+system.time(source("time-dependent-98-tunnel-states.R", echo = TRUE, local = tunnel98))
+system.time(source("time-dependent-phase-type.R",       echo = TRUE, local = phase_type))
+
+res <- 
+  tibble(
+    cycle               = csstd$cycles,
+    diseased.csstd      = csstd$state_membership[, 2],
+    diseased.no_tunnel  = no_tunnel$state_membership[, 2],
+    diseased.tunnel9    = tunnel9$summarised_state_membership[, 2],
+    diseased.tunnel98   = tunnel98$summarised_state_membership[, 2],
+    diseased.phase_type = phase_type$summarised_state_membership[, 2],
+    survival.csstd      = 1 - csstd$state_membership[, 3],
+    survival.no_tunnel  = 1 - no_tunnel$state_membership[, 3],
+    survival.tunnel9    = 1 - tunnel9$summarised_state_membership[, 3],
+    survival.tunnel98   = 1 - tunnel98$summarised_state_membership[, 3],
+    survival.phase_type = 1 - phase_type$summarised_state_membership[, 3]) %>%
+    pivot_longer(cols = -cycle, names_to = c(".value", "method"), names_sep = "\\.") %>%
+  left_join(
+    tibble(
+      method = c("csstd", "no_tunnel", "tunnel9", "tunnel98", "phase_type"),
+      description = c("CSSTD", "No tunnel states", "Nine tunnel states", "98 tunnel states", "Phase-type")))
+
+ggplot(res, aes(x = cycle, y = diseased, colour = description, linetype = description)) +
+  geom_line() +
+  labs(x = "Cycle", y = "Proportion of cohort in Diseased state", colour = NULL, linetype = NULL) +
+  theme_bw()
+ggsave("time-dependent-diseased.png", width = 15.92, height = 15.92 * 3/4, units = "cm", dpi = 600)
+
+ggplot(res, aes(x = cycle, y = survival, colour = description, linetype = description)) +
+  geom_line() +
+  labs(x = "Cycle", y = "Proportion of cohort alive", colour = NULL, linetype = NULL) +
+  theme_bw()
+ggsave("time-dependent-survival.png", width = 15.92, height = 15.92 * 3/4, units = "cm", dpi = 600)
+

time-dependent-no-tunnel-states.R

@@ -0,0 +1,43 @@
+# We will use cycles from 0 to 100
+cycles <- 0:100
+
+# Create an array for the state membership
+state_membership <- array(NA_real_, c(length(cycles), 3), list(cycle = cycles, state = c("Healthy", "Diseased", "Dead")))
+state_membership[1, ] <- c(1, 0, 0)
+
+# Define the survival model parameters
+# - Healthy to Diseased
+lambda_inc  <- 2e-5
+gamma_inc   <- 2.5
+# - Healthy to Dead
+a_general   <- 0.2
+b_general   <- 1e-7
+# - Diseased to Dead
+lambda_surv <- 0.2
+gamma_surv  <- 0.5
+
+# Calculate exponential distribution with same mean survival as Weibull
+# distribution for Diseased to Dead
+mean_surv <- lambda_surv ^ (-1 / gamma_surv) * gamma(1 + 1 / gamma_surv)
+rate_surv <- 1 / mean_surv
+
+# Calculate the hazard rates for the transitions out of the Healthy state
+healthy_haz <- array(NA_real_, c(length(cycles), 2), list(cycle = cycles, to = c("Diseased", "Dead")))
+healthy_haz[, 1] <- lambda_inc * ((cycles + 1) ^ gamma_inc - cycles ^ gamma_inc)
+healthy_haz[, 2] <- b_general / a_general * (exp(a_general * (cycles + 1)) - exp(a_general * cycles))
+
+# Create an array for transition probabilities in each cycle
+transition_probs <- array(NA_real_, c(length(cycles), 3), list(cycle = cycles, transition = c("Healthy->Diseased", "Healthy->Dead", "Diseased->Dead")))
+transition_probs[, 1:2] <- healthy_haz / rowSums(healthy_haz) * (1 - exp(-rowSums(healthy_haz)))
+transition_probs[, 3] <- 1 - exp(-rate_surv)
+
+# Run cohort simulation
+for (i in 2:length(cycles)) {
+  # Update state membership
+  state_membership[i, 1] <- state_membership[i - 1, 1] * (1 - sum(transition_probs[i - 1, 1:2]))
+  state_membership[i, 2] <- state_membership[i - 1, 1] * transition_probs[i - 1, 1] +
+    state_membership[i - 1, 2] * (1 - transition_probs[i - 1, 3])
+  state_membership[i, 3] <- state_membership[i - 1, 1] * transition_probs[i - 1, 2] +
+    state_membership[i - 1, 2] * transition_probs[i - 1, 3] +
+    state_membership[i - 1, 3]
+}

time-dependent-phase-type.R

@@ -0,0 +1,58 @@
+# We will use cycles from 0 to 100
+cycles <- 0:100
+
+# Create an array for the state membership
+state_membership <- array(NA_real_, c(length(cycles), 5), list(cycle = cycles, state = c("Healthy", "DiseasedA", "DiseasedB", "DiseasedC", "Dead")))
+state_membership[1, ] <- c(1, 0, 0, 0, 0)
+
+# Define the survival model parameters
+# - Healthy to Diseased
+lambda_inc  <- 2e-5
+gamma_inc   <- 2.5
+# - Healthy to Dead
+a_general   <- 0.2
+b_general   <- 1e-7
+# - Diseased to Dead
+lambda_surv <- 0.2
+gamma_surv  <- 0.5
+
+# Calculate the hazard rates for the transitions out of the Healthy state
+healthy_haz <- array(NA_real_, c(length(cycles), 2), list(cycle = cycles, to = c("Diseased", "Dead")))
+healthy_haz[, 1] <- lambda_inc * ((cycles + 1) ^ gamma_inc - cycles ^ gamma_inc)
+healthy_haz[, 2] <- b_general / a_general * (exp(a_general * (cycles + 1)) - exp(a_general * cycles))
+
+# Fit a phase-type distribution to the Weibull distribution for disease survival
+fit <- mapfit::phfit.density(
+  ph = mapfit::cf1(3),
+  f = function(x) lambda_surv * gamma_surv * x ^ (gamma_surv - 1) * exp(-lambda_surv * x ^ gamma_surv))
+Q <- Matrix::as.matrix(fit$model@Q)
+P <- expm::expm(Q)
+
+# Create an array for transition probabilities in each cycle
+transition_probs <- array(
+  NA_real_,
+  c(length(cycles), 10),
+  list(
+    cycle = cycles,
+    transition = c("Healthy->DiseasedA", "Healthy->DiseasedB", "Healthy->DiseasedC", "Healthy->Dead", "DiseasedA->DiseasedA", "DiseasedA->DiseasedB", "DiseasedB->DiseasedB", "DiseasedA->DiseasedC", "DiseasedB->DiseasedC", "DiseasedC->DiseasedC")))
+
+# The transition probabilities for the Healthy state can be calculated already
+healthy_to_diseased <- healthy_haz[, 1] / rowSums(healthy_haz) * (1 - exp(-rowSums(healthy_haz)))
+transition_probs[, 1:3] <- sapply(1:3, function(i) healthy_to_diseased * fit$model@alpha[i])
+transition_probs[, 4] <- healthy_haz[, 2] / rowSums(healthy_haz) * (1 - exp(-rowSums(healthy_haz)))
+transition_probs[, 5:10] <- t(sapply(cycles, function(i) P[upper.tri(P, diag = TRUE)]))
+
+# Run the cohort simulation
+for (i in 2:length(cycles)) {
+  TP <- matrix(0, 5, 5)
+  TP[1, 5] <- transition_probs[i - 1, 4]
+  TP[1, 2:4] <- transition_probs[i - 1, 1:3]
+  TP[1, 1] <- 1 - sum(TP[1, 2:5])
+  TP[(row(TP) %in% (2:4)) & (col(TP) %in% (2:4) & (row(TP) <= col(TP)))] <- transition_probs[i - 1, 5:10]
+  TP[2:5, 5] <- 1 - rowSums(TP[2:5, 1:4])
+  
+  state_membership[i, ] <- t(TP) %*% state_membership[i - 1, ]
+}
+
+summarised_state_membership <- cbind(state_membership[, 1], rowSums(state_membership[, 2:4]), state_membership[, 5])
+dimnames(summarised_state_membership) <- list(cycle = cycles, state = c("Healthy", "Diseased", "Dead"))