Correct discounting and EVPI calculation

THR-MCPSA.R

@@ -117,8 +117,8 @@ runPSA <- function(...) {
   V_base[6, 4] <- theta[12]
   V_base[6, 6] <- theta[12]
   
-  for (t in 1:60) {
-    V[, , t] <- V_base %*% diag(rep(c(1.06^(-t), 1.015^(-t)), 3), 6, 6)
+  for (t in 0:59) {
+    V[, , t+1] <- V_base %*% diag(rep(c(1.06^(-t), 1.015^(-t)), 3), 6, 6)
   }
   
   # STANDARD MARKOV COHORT SIMULATION

THR-deltaPSA.R

@@ -125,8 +125,8 @@ V_base[2, 6] <- -expit(theta_mean[12])
 V_base[6, 4] <- expit(theta_mean[12])
 V_base[6, 6] <- expit(theta_mean[12])
 
-for (t in 1:60) {
-  V[, , t] <- V_base %*% diag(rep(c(1.06^(-t), 1.015^(-t)), 3), 6, 6)
+for (t in 0:59) {
+  V[, , t+1] <- V_base %*% diag(rep(c(1.06^(-t), 1.015^(-t)), 3), 6, 6)
 }
 
 # STANDARD MARKOV COHORT SIMULATION
@@ -306,7 +306,7 @@ p_inc_QALYs <- ggplot(g_psa, aes(x = Incremental_QALYs)) +
 (p_standard_costs + p_NP1_costs + p_inc_costs) / (p_standard_QALYs + p_NP1_QALYs + p_inc_QALYs)
 ggsave(here("Delta-PSA comparison.png"), width = 20.28, height = 10.62, units = "cm")
 
-wtp <- seq(0, 10000, by = 100)
+wtp <- seq(0, 10000, by = 50)
 
 ceac_psa <- sapply(wtp, function(lambda) mean(lambda * g_psa$Incremental_QALYs >= g_psa$Incremental_Costs))
 
@@ -324,10 +324,10 @@ ggsave(here("Delta-PSA CEAC.png"), width = 20.28, height = 10.62, units = "cm",
 
 evpi_psa <- sapply(wtp, function(lambda) {
   inmb <- lambda * g_psa$Incremental_QALYs - g_psa$Incremental_Costs
-  evpi <- mean(abs(inmb)) - abs(mean(inmb))
+  evpi <- mean(pmax(0, inmb)) - pmax(0, mean(inmb))
   evpi
 })
-evpi_delta <- (s * sqrt(2/pi) * exp(-m^2 / (2 * s^2)) - m * (1 - 2 * pnorm(m/s))) - abs(m)
+evpi_delta <- m * pnorm(m/s) + s * dnorm(m/s) - pmax(0, m)
 
 ggplot(tibble(x = wtp, evpi.delta = evpi_delta, evpi.psa = evpi_psa), aes(x)) +
   geom_line(aes(y = evpi.delta, linetype = "Delta-PSA", colour = "Delta-PSA")) +