In Xa4P/pacheck: Probabilistic Analysis Check Package

devtools::load_all()
library(knitr)
library(ggplot2)
data(df_pa)

Introduction

This article describes an example workflow for the use of the pacheck package.
This vignette is still in development.

1. Upload of original health economic model inputs and outputs

df_pa <- calculate_nb(df = df_pa,
                      e_int = "t_qaly_d_int",
                      e_comp = "t_qaly_d_comp",
                      c_int = "t_costs_d_int",
                      c_comp = "t_costs_d_comp",
                      wtp = 80000)# calculate net benefits

2. Investigate model inputs and outputs

a. Summary statistics of (user-selected) model inputs and outputs

To examine whether cost inputs are always positive for instance

df <- generate_sum_stats(df_pa)
kable(df)
rm(df)

b. Correlation matrix inputs (and outputs)

generate_cor(df_pa)

c. Inspect types of variables

Perform quick checks on input values

Are utility values/probabilities between 0-1, costs positive etc...

do_quick_check(df = df_pa,
               v_utilities = c("u_pfs", "u_pd"),
               v_costs = c("c_pfs", "c_pd", "c_thx"),
               v_rr = "rr"
               )

Perform quick comparison discounted and undiscounted results

Are utility values/probabilities between 0-1, costs positive etc...

do_discount_check(df = df_pa,
                   v_outcomes = c("t_qaly_comp", "t_qaly_int"),
                   v_outcomes_d = c("t_qaly_d_comp", "t_qaly_d_int")
                   )

Positive variables

Are these variables strictly positive?

check_positive("c_pfs", "c_pd", df = df_pa)

Variables between 0-1

Are these variables strictly positive?

check_binary("u_pfs", "p_pfspd", df = df_pa)

Sum of probabilities

To check whether the sum is lower than or equal to, or equal to 1

check_sum_probs("p_pfspd", "p_pfsd", df = df_pa, check = "lower") # output is a text

d. Histogram, density distribution of (user-selected) model inputs and outputs

To visually investigate the parameter distributions

Single parameter

p_1 <- vis_1_param(df = df_pa,
            param = "u_pfs",
            binwidth = NULL,
            type = "histogram",
            dist = c("norm", "beta", "gamma", "lnorm"))
p_1

paste("Probability to be in user-defined range: ",
      check_range(df_pa,
                  param = "u_pfs",
                  min_val = 0.77,
                  max_val = 0.80
                  ), "%") # add this to plot? + lines of min/ max on plot?

Two parameters

p_2p <- vis_2_params(df = df_pa,
                     param_1 = "u_pfs",
                     param_2 = "u_pd",
                     slope = 1,
                     check = "param_2 > param_1")
p_2p

e. Possibility to fit user-selected distributions (beta, gamma, lognormal, normal) + visualisation + parameters of the fitted distribution.

To cross check with parameters reported in documentation/report/ journal article, as an implementation check
- To do XP: Add probabilistic mean value per distribution

p_2 <- vis_1_param(df = df_pa,
            param = "u_pfs",
            binwidth = NULL,
            type = "density",
            dist = c("norm", "beta", "gamma", "lnorm"),
            user_dist = "beta",
            user_param_1 = 0.8,
            user_param_2 = 0.2,
            user_mean = 0.75)
p_2

Distributions' parameters & statistical fit

l_dist <- fit_dist(df = df_pa,
                   param = "u_pfs",
                   dist = c("norm", "beta", "gamma", "lnorm"))
l_dist[[1]]
l_dist[[2]]

3. Investigate model outputs

a. Incremental cost-effectiveness plane & summary

p_3 <- plot_ice(
  df = df_pa,
  e_int = "t_qaly_d_int",
  e_comp = "t_qaly_d_comp",
  c_int = "t_costs_d_int",
  c_comp = "t_costs_d_comp",
  wtp = 8000
)
p_3

summary_ice(df_pa,
            e_int = "t_qaly_d_int",
            e_comp = "t_qaly_d_comp",
            c_int = "t_costs_d_int",
            c_comp = "t_costs_d_comp")

b. Cost-effectiveness acceptability curve.

df_ceac <- calculate_ceac(df = df_pa,
                     e_int = "t_qaly_d_int",
                     e_comp = "t_qaly_d_comp",
                     c_int = "t_costs_d_int",
                     c_comp = "t_costs_d_comp")

plot_ceac(df = df_ceac,
          wtp = "WTP_threshold")
df_ceac

c. Histogram and density distribution of total and incremental costs and effects.

Can use the function above!

d. Convergence graph of outcomes

plot_convergence(df = df_pa,
                 param = "iNMB"
                 )

e. Check whether the mean quality of life outcome of each iteration remain between the maximum and minimum utility values of the specific iteration.

check_mean_qol(df = df_pa,
               t_ly = "t_ly_comp",
               t_qaly = "t_qaly_comp",
               u_values = c("u_pfs", "u_pd")
               )

4. Investigate relation between inputs and outputs

Single

• To do XP: Add R^2^

lm_rr <- fit_lm_metamodel(df = df_pa,
                          x_vars = "rr",
                          y_var = "inc_qaly")
lm_pred <- unlist(predict(lm_rr$fit, data.frame(rr = df_pa$rr)))
df_obs_pred <- data.frame(
  Values = df_pa$rr,
  Observed = df_pa$inc_qaly,
  Predicted = lm_pred
)
ggplot(data = df_obs_pred, aes(x = Values, y = Observed)) +
  geom_point(shape = 1, colour = "lightgrey") +
  geom_smooth(method = "lm") +
  theme_bw()
plot(lm_rr$fit)

Multiple

lm_full <- fit_lm_metamodel(df = df_pa,
                            x_vars = c("rr",
                                       "u_pfs",
                                       "u_pd",
                                       "c_pfs",
                                       "c_pd",
                                       "c_thx",
                                       "p_pfspd",
                                       "p_pfsd",
                                       "p_pdd"),
                            y_var = "inc_qaly")
lm_pred_full <- predict(lm_full$fit, data.frame(rr = df_pa$rr,
                            u_pfs = mean(df_pa$u_pfs),
                            u_pd = mean(df_pa$u_pd),
                            c_pfs = mean(df_pa$c_pfs),
                            c_pd = mean(df_pa$c_pd),
                            c_thx = mean(df_pa$c_thx),
                            p_pfspd = mean(df_pa$p_pfspd),
                            p_pfsd = mean(df_pa$p_pfsd),
                            p_pdd = mean(df_pa$p_pdd)))

df_obs_pred_full <- data.frame(
  Values = df_pa$rr,
  Observed = df_pa$inc_qaly,
  Predicted = lm_pred_full
)

ggplot(data = df_obs_pred_full, aes(x = Values, y = Observed)) +
  geom_point(shape = 1, colour = "lightgrey") +
  geom_line(data = df_obs_pred_full, aes(x = Values, y = Predicted), colour = "blue") +
  theme_bw()

plot(lm_full$fit)

5. Predictions based on metamodel

# Predicting using metamodel fitted above
predict(
  lm_full$fit,
  data.frame(
    rr = 0.9,
    u_pfs = 0.75,
    u_pd = 0.4,
    c_pfs = 1000,
    c_pd = 2000,
    c_thx = 0,
    p_pfspd = 0.5,
    p_pfsd = 0.3,
    p_pdd = 0.4
  )
)

Xa4P/pacheck documentation built on April 14, 2025, 1:51 p.m.