Remove comments and use here()

THR-deltaPSA.R

@@ -1,19 +1,6 @@
 library(tidyverse)
 library(here)
 
-# optim(par = c(0, 1),
-#       fn  = function(params, beta_shape1, beta_shape2) {
-#         a <- beta_shape1
-#         b <- beta_shape2
-#         beta_mean <- a / (a + b)
-#         beta_var  <- a * b / ((a + b) ^ 2 * (a + b + 1))
-#         moments <- logitnorm::momentsLogitnorm(params[1], params[2])
-#         (moments[1] - beta_mean) ^ 2 + (moments[2] - beta_var) ^ 2
-#       },
-#       beta_shape1 = 69.7,
-#       beta_shape2 = 162.633333,
-#       method = "BFGS")
-
 theta_mean <- c(-5.490935,
                 -0.0367022,
                 0.768536,
@@ -40,21 +27,6 @@ theta_var[10,10] <- 0.235360^2
 theta_var[11,11] <- 0.2147719^2
 theta_var[12,12] <- 0.1433964^2
 
-# ggplot(tibble(x = c(0, 0.2)), aes(x)) +
-#   stat_function(aes(colour = "Beta(4, 96)"), fun = dbeta, args = list(shape1 = 4, shape2 = 96)) +
-#   stat_function(aes(colour = "LogitN(-3.29, 0.49)"), fun = logitnorm::dlogitnorm, args = list(mu = -3.2869827, sigma = 0.4889407)) +
-#   labs(x = "Re-revision risk", y = "Probability density function", colour = NULL) +
-#   theme_bw()
-# ggsave(filename = "logitnormal-approx.png", width = 20.28, height = 10.62, units = "cm", scale = 0.8)
-# 
-# 
-# ggplot(tibble(x = c(0, 20000)), aes(x)) +
-#   stat_function(aes(colour = "Gamma(12.67, 417.67)"), fun = dgamma, args = list(shape = 12.67, scale = 417.67)) +
-#   stat_function(aes(colour = "LogN(8.54, 0.28)"), fun = dlnorm, args = list(meanlog = 8.53636, sdlog = 0.2755688)) +
-#   labs(x = "Cost of revision THR", y = "Probability density function", colour = NULL) +
-#   theme_bw()
-# ggsave(filename = "lognormal-approx.png", width = 20.28, height = 10.62, units = "cm", scale = 0.8)
-
 state_labels  <- c("Standard_SuccessP", "Standard_Revision", "Standard_SuccessR", "Standard_Death",
                   "NP1_SuccessP", "NP1_Revision", "NP1_SuccessR", "NP1_Death")
 cycle_labels  <- 1:60
@@ -290,7 +262,7 @@ g_variance <- dg %*% theta_var %*% t(dg)
 g_sd <- sqrt(diag(g_variance))
 
 # COMPARE WITH PSA
-g_psa <- readRDS("THR-MCPSA.rds")
+g_psa <- readRDS(here("THR-MCPSA.rds"))
 
 library(patchwork)
 
@@ -332,7 +304,7 @@ p_inc_QALYs <- ggplot(g_psa, aes(x = Incremental_QALYs)) +
   theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank())
 
 (p_standard_costs + p_NP1_costs + p_inc_costs) / (p_standard_QALYs + p_NP1_QALYs + p_inc_QALYs)
-ggsave("Delta-PSA comparison.png", width = 20.28, height = 10.62, units = "cm")
+ggsave(here("Delta-PSA comparison.png"), width = 20.28, height = 10.62, units = "cm")
 
 wtp <- seq(0, 10000, by = 100)
 
@@ -348,7 +320,7 @@ ggplot(tibble(x = wtp, ceac.delta = pnorm(m/s), ceac.psa = ceac_psa), aes(x)) +
   scale_y_continuous("Probability NP1 is cost-effective", labels = scales::percent_format()) +
   scale_colour_discrete("Method") +
   scale_linetype_discrete("Method")
-ggsave("Delta-PSA CEAC.png", width = 20.28, height = 10.62, units = "cm", scale = 0.8)
+ggsave(here("Delta-PSA CEAC.png"), width = 20.28, height = 10.62, units = "cm", scale = 0.8)
 
 evpi_psa <- sapply(wtp, function(lambda) {
   inmb <- lambda * g_psa$Incremental_QALYs - g_psa$Incremental_Costs
@@ -364,4 +336,4 @@ ggplot(tibble(x = wtp, evpi.delta = evpi_delta, evpi.psa = evpi_psa), aes(x)) +
   scale_y_continuous("Expected value of perfect information", labels = scales::dollar_format(prefix = "£")) +
   scale_colour_discrete("Method") +
   scale_linetype_discrete("Method")
-ggsave("Delta-PSA EVPI.png", width = 20.28, height = 10.62, units = "cm", scale = 0.8)
+ggsave(here("Delta-PSA EVPI.png"), width = 20.28, height = 10.62, units = "cm", scale = 0.8)