Exam Topics in Economic Evaluation: Ruth Edigin
HEVAL5200, 13 December 2024
Set up packages, data and working directory
rm(list=ls(all=TRUE)) # clear workspacelist_of_packages <- c("dplyr", "haven", "foreign", "ggplot2", "gmodels", "ggpubr", "knitr", "formattable", "coin", "broom", "margins", "sandwich", "lmtest", "gtsummary", "brms", "MASS", "boot", "openxlsx")
new_packages <- list_of_packages[!(list_of_packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
lapply(list_of_packages, library, character.only = TRUE)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: survival## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: Rcpp## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:margins':
##
## marginal_effects## The following object is masked from 'package:survival':
##
## kidney## The following object is masked from 'package:stats':
##
## ar##
## Attaching package: 'MASS'## The following object is masked from 'package:gtsummary':
##
## select## The following object is masked from 'package:formattable':
##
## area## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'boot'## The following object is masked from 'package:survival':
##
## aml## [[1]]
## [1] "dplyr" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "haven" "dplyr" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[3]]
## [1] "foreign" "haven" "dplyr" "stats" "graphics" "grDevices"
## [7] "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "ggplot2" "foreign" "haven" "dplyr" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "gmodels" "ggplot2" "foreign" "haven" "dplyr" "stats"
## [7] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[6]]
## [1] "ggpubr" "gmodels" "ggplot2" "foreign" "haven" "dplyr"
## [7] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [13] "base"
##
## [[7]]
## [1] "knitr" "ggpubr" "gmodels" "ggplot2" "foreign" "haven"
## [7] "dplyr" "stats" "graphics" "grDevices" "utils" "datasets"
## [13] "methods" "base"
##
## [[8]]
## [1] "formattable" "knitr" "ggpubr" "gmodels" "ggplot2"
## [6] "foreign" "haven" "dplyr" "stats" "graphics"
## [11] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "coin" "survival" "formattable" "knitr" "ggpubr"
## [6] "gmodels" "ggplot2" "foreign" "haven" "dplyr"
## [11] "stats" "graphics" "grDevices" "utils" "datasets"
## [16] "methods" "base"
##
## [[10]]
## [1] "broom" "coin" "survival" "formattable" "knitr"
## [6] "ggpubr" "gmodels" "ggplot2" "foreign" "haven"
## [11] "dplyr" "stats" "graphics" "grDevices" "utils"
## [16] "datasets" "methods" "base"
##
## [[11]]
## [1] "margins" "broom" "coin" "survival" "formattable"
## [6] "knitr" "ggpubr" "gmodels" "ggplot2" "foreign"
## [11] "haven" "dplyr" "stats" "graphics" "grDevices"
## [16] "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "sandwich" "margins" "broom" "coin" "survival"
## [6] "formattable" "knitr" "ggpubr" "gmodels" "ggplot2"
## [11] "foreign" "haven" "dplyr" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "lmtest" "zoo" "sandwich" "margins" "broom"
## [6] "coin" "survival" "formattable" "knitr" "ggpubr"
## [11] "gmodels" "ggplot2" "foreign" "haven" "dplyr"
## [16] "stats" "graphics" "grDevices" "utils" "datasets"
## [21] "methods" "base"
##
## [[14]]
## [1] "gtsummary" "lmtest" "zoo" "sandwich" "margins"
## [6] "broom" "coin" "survival" "formattable" "knitr"
## [11] "ggpubr" "gmodels" "ggplot2" "foreign" "haven"
## [16] "dplyr" "stats" "graphics" "grDevices" "utils"
## [21] "datasets" "methods" "base"
##
## [[15]]
## [1] "brms" "Rcpp" "gtsummary" "lmtest" "zoo"
## [6] "sandwich" "margins" "broom" "coin" "survival"
## [11] "formattable" "knitr" "ggpubr" "gmodels" "ggplot2"
## [16] "foreign" "haven" "dplyr" "stats" "graphics"
## [21] "grDevices" "utils" "datasets" "methods" "base"
##
## [[16]]
## [1] "MASS" "brms" "Rcpp" "gtsummary" "lmtest"
## [6] "zoo" "sandwich" "margins" "broom" "coin"
## [11] "survival" "formattable" "knitr" "ggpubr" "gmodels"
## [16] "ggplot2" "foreign" "haven" "dplyr" "stats"
## [21] "graphics" "grDevices" "utils" "datasets" "methods"
## [26] "base"
##
## [[17]]
## [1] "boot" "MASS" "brms" "Rcpp" "gtsummary"
## [6] "lmtest" "zoo" "sandwich" "margins" "broom"
## [11] "coin" "survival" "formattable" "knitr" "ggpubr"
## [16] "gmodels" "ggplot2" "foreign" "haven" "dplyr"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
##
## [[18]]
## [1] "openxlsx" "boot" "MASS" "brms" "Rcpp"
## [6] "gtsummary" "lmtest" "zoo" "sandwich" "margins"
## [11] "broom" "coin" "survival" "formattable" "knitr"
## [16] "ggpubr" "gmodels" "ggplot2" "foreign" "haven"
## [21] "dplyr" "stats" "graphics" "grDevices" "utils"
## [26] "datasets" "methods" "base"#Set Working directory to where file is saved
setwd("C:/Users/ruthe/OneDrive/Documents/Exam")
load("Exam_data_2024-2.Rdata")
View(Exam_data_2024)Learning about Dataset
group_labels <- c("Surgery", "Control")
Exam_data_2024$Group <- factor(Exam_data_2024$Group, labels = group_labels)table(Exam_data_2024$Group)##
## Surgery Control
## 25 25freq_table_group <- table(Exam_data_2024$Group)
proption_table_group <- prop.table(freq_table_group)
print(proption_table_group)##
## Surgery Control
## 0.5 0.5Gender_labels <- c("Female", "Male")
Exam_data_2024$Gender <- factor(Exam_data_2024$Gender, labels = Gender_labels)
table(Exam_data_2024$Gender)##
## Female Male
## 44 6freq_table_gender <- table(Exam_data_2024$Gender)
proption_table_gender <- prop.table(freq_table_gender)
print(proption_table_gender)##
## Female Male
## 0.88 0.12# Summary of Groups
table(Exam_data_2024$Group)##
## Surgery Control
## 25 25CrossTable(Exam_data_2024$Gender, Exam_data_2024$Group, prop.r = TRUE)##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 50
##
##
## | Exam_data_2024$Group
## Exam_data_2024$Gender | Surgery | Control | Row Total |
## ----------------------|-----------|-----------|-----------|
## Female | 20 | 24 | 44 |
## | 0.182 | 0.182 | |
## | 0.455 | 0.545 | 0.880 |
## | 0.800 | 0.960 | |
## | 0.400 | 0.480 | |
## ----------------------|-----------|-----------|-----------|
## Male | 5 | 1 | 6 |
## | 1.333 | 1.333 | |
## | 0.833 | 0.167 | 0.120 |
## | 0.200 | 0.040 | |
## | 0.100 | 0.020 | |
## ----------------------|-----------|-----------|-----------|
## Column Total | 25 | 25 | 50 |
## | 0.500 | 0.500 | |
## ----------------------|-----------|-----------|-----------|
##
## # Testing continuous variables - AGE
summary(Exam_data_2024$Age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 59.00 65.25 72.00 72.64 78.75 88.00ggplot(Exam_data_2024, aes(x = Age)) + geom_histogram(binwidth = 5) + facet_wrap(~ Group)
# T-test for Age (assuming normal distribution)
t.test(Age ~ Group, data = Exam_data_2024)##
## Welch Two Sample t-test
##
## data: Age by Group
## t = -0.35907, df = 44.968, p-value = 0.7212
## alternative hypothesis: true difference in means between group Surgery and group Control is not equal to 0
## 95 percent confidence interval:
## -5.816254 4.056254
## sample estimates:
## mean in group Surgery mean in group Control
## 72.20 73.08#rank-sum test (not normal distribution)
wilcox_test(Age ~ Group, data = Exam_data_2024)##
## Asymptotic Wilcoxon-Mann-Whitney Test
##
## data: Age by Group (Surgery, Control)
## Z = -0.35951, p-value = 0.7192
## alternative hypothesis: true mu is not equal to 0# Generate continuous variable into categories - AGE
Exam_data_2024 <- Exam_data_2024 %>%
mutate(age_group = case_when(
Age < 65 ~ 1,
Age >= 65 & Age < 70 ~ 2,
Age >= 70 & Age < 75 ~ 3,
Age >= 75 & Age < 80 ~ 4,
Age >= 80 & Age < 85 ~ 5,
Age >= 85 ~ 6
))
# Create age group labels
age_labels <- c("<65", "65-69", "70-74", "75-79", "80-84", "85+")
Exam_data_2024$age_group <- factor(Exam_data_2024$age_group, labels = age_labels)
# Tabulate age groups
Exam_data_2024$age_group <- factor(
Exam_data_2024$age_group,
labels = age_labels # Assign labels
)
# Assuming 'Group' has levels that you want to label (e.g., 1, 2, 3)
Exam_data_2024$Group <- factor(
Exam_data_2024$Group,
labels = group_labels # Assign labels
)
# Create the CrossTable with labeled factors
CrossTable(
Exam_data_2024$age_group,
Exam_data_2024$Group,
prop.r = TRUE,
chisq = TRUE
)## Warning in chisq.test(t, correct = FALSE, ...): Chi-squared approximation may
## be incorrect##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 50
##
##
## | Exam_data_2024$Group
## Exam_data_2024$age_group | Surgery | Control | Row Total |
## -------------------------|-----------|-----------|-----------|
## <65 | 4 | 6 | 10 |
## | 0.200 | 0.200 | |
## | 0.400 | 0.600 | 0.200 |
## | 0.160 | 0.240 | |
## | 0.080 | 0.120 | |
## -------------------------|-----------|-----------|-----------|
## 65-69 | 4 | 5 | 9 |
## | 0.056 | 0.056 | |
## | 0.444 | 0.556 | 0.180 |
## | 0.160 | 0.200 | |
## | 0.080 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 70-74 | 6 | 2 | 8 |
## | 1.000 | 1.000 | |
## | 0.750 | 0.250 | 0.160 |
## | 0.240 | 0.080 | |
## | 0.120 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 75-79 | 6 | 5 | 11 |
## | 0.045 | 0.045 | |
## | 0.545 | 0.455 | 0.220 |
## | 0.240 | 0.200 | |
## | 0.120 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 80-84 | 4 | 2 | 6 |
## | 0.333 | 0.333 | |
## | 0.667 | 0.333 | 0.120 |
## | 0.160 | 0.080 | |
## | 0.080 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 85+ | 1 | 5 | 6 |
## | 1.333 | 1.333 | |
## | 0.167 | 0.833 | 0.120 |
## | 0.040 | 0.200 | |
## | 0.020 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## Column Total | 25 | 25 | 50 |
## | 0.500 | 0.500 | |
## -------------------------|-----------|-----------|-----------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 5.935354 d.f. = 5 p = 0.3125602
##
##
## table(Exam_data_2024$age_group)##
## <65 65-69 70-74 75-79 80-84 85+
## 10 9 8 11 6 6table(Exam_data_2024$age_group)##
## <65 65-69 70-74 75-79 80-84 85+
## 10 9 8 11 6 6Exam_data_2024$age_group <- factor(
Exam_data_2024$age_group,
labels = age_labels # Assign labels
)
# Create the CrossTable with labeled factors
CrossTable(
Exam_data_2024$age_group,
Exam_data_2024$Group,
prop.r = TRUE,
chisq = TRUE
)## Warning in chisq.test(t, correct = FALSE, ...): Chi-squared approximation may
## be incorrect##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 50
##
##
## | Exam_data_2024$Group
## Exam_data_2024$age_group | Surgery | Control | Row Total |
## -------------------------|-----------|-----------|-----------|
## <65 | 4 | 6 | 10 |
## | 0.200 | 0.200 | |
## | 0.400 | 0.600 | 0.200 |
## | 0.160 | 0.240 | |
## | 0.080 | 0.120 | |
## -------------------------|-----------|-----------|-----------|
## 65-69 | 4 | 5 | 9 |
## | 0.056 | 0.056 | |
## | 0.444 | 0.556 | 0.180 |
## | 0.160 | 0.200 | |
## | 0.080 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 70-74 | 6 | 2 | 8 |
## | 1.000 | 1.000 | |
## | 0.750 | 0.250 | 0.160 |
## | 0.240 | 0.080 | |
## | 0.120 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 75-79 | 6 | 5 | 11 |
## | 0.045 | 0.045 | |
## | 0.545 | 0.455 | 0.220 |
## | 0.240 | 0.200 | |
## | 0.120 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 80-84 | 4 | 2 | 6 |
## | 0.333 | 0.333 | |
## | 0.667 | 0.333 | 0.120 |
## | 0.160 | 0.080 | |
## | 0.080 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 85+ | 1 | 5 | 6 |
## | 1.333 | 1.333 | |
## | 0.167 | 0.833 | 0.120 |
## | 0.040 | 0.200 | |
## | 0.020 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## Column Total | 25 | 25 | 50 |
## | 0.500 | 0.500 | |
## -------------------------|-----------|-----------|-----------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 5.935354 d.f. = 5 p = 0.3125602
##
##
## # Sort by Group
Exam_data_2024 %>% group_by(Group) %>% summarise(age_group_counts = list(table(age_group)))## # A tibble: 2 × 2
## Group age_group_counts
## <fct> <list>
## 1 Surgery <table [6]>
## 2 Control <table [6]># Test for differences (categorical variable)
CrossTable(Exam_data_2024$age_group, Exam_data_2024$Group, prop.c = TRUE, chisq = TRUE)## Warning in chisq.test(t, correct = FALSE, ...): Chi-squared approximation may
## be incorrect##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 50
##
##
## | Exam_data_2024$Group
## Exam_data_2024$age_group | Surgery | Control | Row Total |
## -------------------------|-----------|-----------|-----------|
## <65 | 4 | 6 | 10 |
## | 0.200 | 0.200 | |
## | 0.400 | 0.600 | 0.200 |
## | 0.160 | 0.240 | |
## | 0.080 | 0.120 | |
## -------------------------|-----------|-----------|-----------|
## 65-69 | 4 | 5 | 9 |
## | 0.056 | 0.056 | |
## | 0.444 | 0.556 | 0.180 |
## | 0.160 | 0.200 | |
## | 0.080 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 70-74 | 6 | 2 | 8 |
## | 1.000 | 1.000 | |
## | 0.750 | 0.250 | 0.160 |
## | 0.240 | 0.080 | |
## | 0.120 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 75-79 | 6 | 5 | 11 |
## | 0.045 | 0.045 | |
## | 0.545 | 0.455 | 0.220 |
## | 0.240 | 0.200 | |
## | 0.120 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 80-84 | 4 | 2 | 6 |
## | 0.333 | 0.333 | |
## | 0.667 | 0.333 | 0.120 |
## | 0.160 | 0.080 | |
## | 0.080 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 85+ | 1 | 5 | 6 |
## | 1.333 | 1.333 | |
## | 0.167 | 0.833 | 0.120 |
## | 0.040 | 0.200 | |
## | 0.020 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## Column Total | 25 | 25 | 50 |
## | 0.500 | 0.500 | |
## -------------------------|-----------|-----------|-----------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 5.935354 d.f. = 5 p = 0.3125602
##
##
## BASELINE TABEL CHARACTERISTICS AND TESTSContinuous variables: t-test (normal distribution), not normal - wilcox ranksum
Categorical variables: chisq
summary_table <- Exam_data_2024 %>%
dplyr::select(Gender, age_group, Time0, Group) %>%
tbl_summary(
by =Group,
statistic = list(
all_continuous() ~ "{mean} ({sd})", # For continuous variables
all_categorical() ~ "{n} ({p}%)" # For categorical variables
)
)
summary_table| Characteristic | Surgery N = 251 |
Control N = 251 |
|---|---|---|
| Gender | ||
| Female | 20 (80%) | 24 (96%) |
| Male | 5 (20%) | 1 (4.0%) |
| age_group | ||
| <65 | 4 (16%) | 6 (24%) |
| 65-69 | 4 (16%) | 5 (20%) |
| 70-74 | 6 (24%) | 2 (8.0%) |
| 75-79 | 6 (24%) | 5 (20%) |
| 80-84 | 4 (16%) | 2 (8.0%) |
| 85+ | 1 (4.0%) | 5 (20%) |
| 15D at Baseline | 0.87 (0.07) | 0.83 (0.11) |
| 1 n (%); Mean (SD) | ||
# Test for differences (categorical variable)
CrossTable(Exam_data_2024$Gender, Exam_data_2024$Group, prop.c = TRUE, chisq = TRUE)## Warning in chisq.test(t, correct = TRUE, ...): Chi-squared approximation may be
## incorrect## Warning in chisq.test(t, correct = FALSE, ...): Chi-squared approximation may
## be incorrect##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 50
##
##
## | Exam_data_2024$Group
## Exam_data_2024$Gender | Surgery | Control | Row Total |
## ----------------------|-----------|-----------|-----------|
## Female | 20 | 24 | 44 |
## | 0.182 | 0.182 | |
## | 0.455 | 0.545 | 0.880 |
## | 0.800 | 0.960 | |
## | 0.400 | 0.480 | |
## ----------------------|-----------|-----------|-----------|
## Male | 5 | 1 | 6 |
## | 1.333 | 1.333 | |
## | 0.833 | 0.167 | 0.120 |
## | 0.200 | 0.040 | |
## | 0.100 | 0.020 | |
## ----------------------|-----------|-----------|-----------|
## Column Total | 25 | 25 | 50 |
## | 0.500 | 0.500 | |
## ----------------------|-----------|-----------|-----------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 3.030303 d.f. = 1 p = 0.08172275
##
## Pearson's Chi-squared test with Yates' continuity correction
## ------------------------------------------------------------
## Chi^2 = 1.704545 d.f. = 1 p = 0.1916946
##
## CrossTable(Exam_data_2024$age_group, Exam_data_2024$Group, prop.c = TRUE, chisq = TRUE)## Warning in chisq.test(t, correct = FALSE, ...): Chi-squared approximation may
## be incorrect##
##
## Cell Contents
## |-------------------------|
## | N |
## | Chi-square contribution |
## | N / Row Total |
## | N / Col Total |
## | N / Table Total |
## |-------------------------|
##
##
## Total Observations in Table: 50
##
##
## | Exam_data_2024$Group
## Exam_data_2024$age_group | Surgery | Control | Row Total |
## -------------------------|-----------|-----------|-----------|
## <65 | 4 | 6 | 10 |
## | 0.200 | 0.200 | |
## | 0.400 | 0.600 | 0.200 |
## | 0.160 | 0.240 | |
## | 0.080 | 0.120 | |
## -------------------------|-----------|-----------|-----------|
## 65-69 | 4 | 5 | 9 |
## | 0.056 | 0.056 | |
## | 0.444 | 0.556 | 0.180 |
## | 0.160 | 0.200 | |
## | 0.080 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 70-74 | 6 | 2 | 8 |
## | 1.000 | 1.000 | |
## | 0.750 | 0.250 | 0.160 |
## | 0.240 | 0.080 | |
## | 0.120 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 75-79 | 6 | 5 | 11 |
## | 0.045 | 0.045 | |
## | 0.545 | 0.455 | 0.220 |
## | 0.240 | 0.200 | |
## | 0.120 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## 80-84 | 4 | 2 | 6 |
## | 0.333 | 0.333 | |
## | 0.667 | 0.333 | 0.120 |
## | 0.160 | 0.080 | |
## | 0.080 | 0.040 | |
## -------------------------|-----------|-----------|-----------|
## 85+ | 1 | 5 | 6 |
## | 1.333 | 1.333 | |
## | 0.167 | 0.833 | 0.120 |
## | 0.040 | 0.200 | |
## | 0.020 | 0.100 | |
## -------------------------|-----------|-----------|-----------|
## Column Total | 25 | 25 | 50 |
## | 0.500 | 0.500 | |
## -------------------------|-----------|-----------|-----------|
##
##
## Statistics for All Table Factors
##
##
## Pearson's Chi-squared test
## ------------------------------------------------------------
## Chi^2 = 5.935354 d.f. = 5 p = 0.3125602
##
##
## COSTS# Checking summary statistics of AGE
summary(Exam_data_2024$Age)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 59.00 65.25 72.00 72.64 78.75 88.00# Checking summary statistics of AGE_Group
summary(Exam_data_2024$age_group)## <65 65-69 70-74 75-79 80-84 85+
## 10 9 8 11 6 6# Checking summary statistics of gender
summary(Exam_data_2024$Gender)## Female Male
## 44 6# Checking summary statistics of Group
summary(Exam_data_2024$Group)## Surgery Control
## 25 25# Initial hospital stay
Exam_data_2024%>%
group_by(Group) %>%
summarize(LOS_baseline= list(summary(LOS_baseline)),
LOS_baseline= list(LOS_baseline)) %>%
ungroup()## # A tibble: 2 × 2
## Group LOS_baseline
## <fct> <list>
## 1 Surgery <list [1]>
## 2 Control <list [1]># Histogram of LOS by Group
ggplot(Exam_data_2024, aes(x = LOS_baseline)) +
geom_histogram(bins = 30) +
facet_wrap(~ Group) +
labs(title = "Histogram of LOS_baseline by Group")
# Mann-Whitney test for LOS by Group
wilcox.test(LOS12 ~ Group, data = Exam_data_2024)## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: LOS12 by Group
## W = 374, p-value = 0.01591
## alternative hypothesis: true location shift is not equal to 0wilcox.test(LOS6 ~ Group, data = Exam_data_2024)## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: LOS6 by Group
## W = 465.5, p-value = 0.002911
## alternative hypothesis: true location shift is not equal to 0wilcox.test(LOS_baseline~ Group, data = Exam_data_2024)## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties##
## Wilcoxon rank sum test with continuity correction
##
## data: LOS_baseline by Group
## W = 465.5, p-value = 0.002911
## alternative hypothesis: true location shift is not equal to 0Cost healthcare perspective
# Cost length of stay baseline, initial
LOS_baseline_day <-500
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_LOS12 = if_else(
is.na(LOS12),
if_else(!is.na(LOS6), LOS6,
if_else(!is.na(LOS_baseline), LOS_baseline,0)), LOS12),
Cost_hospital_YEAR = Adjusted_LOS12 * LOS_baseline_day)
# Summarize the costs by group
Exam_data_2024_summary <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_hospital_YEAR, na.rm = TRUE),
Median_Cost = median(Cost_hospital_YEAR, na.rm = TRUE),
Total_Cost = sum(Cost_hospital_YEAR, na.rm = TRUE)
)
View(Exam_data_2024_summary)# Cost length of stay hospital stay
LOS_ARG_COST_YEAR <- 500
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_LOSarg12 = if_else(
is.na(LOSARG12),
if_else(!is.na(LOSARG6), LOSARG6,
if_else(!is.na(LOSARG3), LOSARG3,0)), LOSARG12),
Cost_hospital_arg12 = Adjusted_LOSarg12 * LOS_ARG_COST_YEAR)
# Summarize the costs by group
Exam_data_2024_summary <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_hospital_arg12, na.rm = TRUE),
Median_Cost = median(Cost_hospital_arg12, na.rm = TRUE),
Total_Cost = sum(Cost_hospital_arg12, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost length of stay, hospital stay, re-operation
LOS_ARG_REOP_YEAR <- 500
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_LOS_ARG12_reop = if_else(
is.na(LOS_ARG12_reop),
if_else(!is.na(LOS_ARG6_reop), LOS_ARG6_reop,
if_else(!is.na(LOS_ARG3_reop), LOS_ARG3_reop,0)), LOS_ARG12_reop),
Cost_ARG_12_months_reop = Adjusted_LOS_ARG12_reop * LOS_ARG_REOP_YEAR)
# Summarize the costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_ARG_12_months_reop, na.rm = TRUE),
Median_Cost = median(Cost_ARG_12_months_reop, na.rm = TRUE),
Total_Cost = sum(Cost_ARG_12_months_reop, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost rehabilitation
LOS_REHAB_YEAR <- 350
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_LOSrehab12 = if_else(
is.na(LOSrehab12),
if_else(!is.na(LOSrehab6), LOSrehab6,
if_else(!is.na(LOSrehab3), LOSrehab3,0)), LOSrehab12),
Cost_LOS_REHAB = Adjusted_LOSrehab12 * LOS_REHAB_YEAR)
# Summarize the costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_LOS_REHAB, na.rm = TRUE),
Median_Cost = median(Cost_LOS_REHAB, na.rm = TRUE),
Total_Cost = sum(Cost_LOS_REHAB, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of nursing home
LOS_nursing_day <- 250
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_LOSnursh12 = if_else(
is.na(LOSnursh12),
if_else(!is.na(LOSnursh6), LOSnursh6,
if_else(!is.na(LOSnursh3), LOSnursh3,0)),LOSnursh12),
Cost_LOS_nursing_year = Adjusted_LOSnursh12 * LOS_nursing_day)
# Summarize the nursing home costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_LOS_nursing_year, na.rm = TRUE),
Median_Cost = median(Cost_LOS_nursing_year, na.rm = TRUE),
Total_Cost = sum(Cost_LOS_nursing_year, na.rm = TRUE) )
View(Exam_data_2024_summary)# Cost of surgery material
Total_surg_material_costs <-1500
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_Material_surgery12 = if_else(
is.na(Material_surgery12),
if_else(!is.na(Material_surgery6), Material_surgery6,
if_else(!is.na(Material_surgery), Material_surgery,0)), Material_surgery12),
Cost_surgery_material_year = Adjusted_Material_surgery12 * Total_surg_material_costs)
# Summarize the surgery material by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_surgery_material_year, na.rm = TRUE),
Median_Cost = median(Cost_surgery_material_year, na.rm = TRUE),
Total_Cost = sum(Cost_surgery_material_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of CT scan
CT_scan <-120
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_CT12 = if_else(
is.na(CT12),
if_else(!is.na(CT6), CT6,
if_else(!is.na(CT_baseline), CT_baseline,0)), CT12),
Cost_CT_year = Adjusted_CT12 * CT_scan)
# Summarize the rehab costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_CT_year, na.rm = TRUE),
Median_Cost = median(Cost_CT_year, na.rm = TRUE),
Total_Cost = sum(Cost_CT_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of XRAY scan
XRAY_scan <-75
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_Xray12 = if_else(
is.na(Xray12),
if_else(!is.na(Xray6), Xray6,
if_else(!is.na(Xray3), Xray3,0)), Xray12),
Cost_XRAY_year = Adjusted_Xray12 * XRAY_scan)
# Summarize the XRAY costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_XRAY_year, na.rm = TRUE),
Median_Cost = median(Cost_XRAY_year, na.rm = TRUE),
Total_Cost = sum(Cost_XRAY_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of physio visit
Physio_visit <- 35
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_NOphysio12 = if_else(
is.na(NOphysio12),
if_else(!is.na(NOphysio6), NOphysio6,
if_else(!is.na(NOphysio3), NOphysio3,0)), NOphysio12),
Cost_physio_visit_year = Adjusted_NOphysio12 * Physio_visit)
# Summarize the rehab costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_physio_visit_year, na.rm = TRUE),
Median_Cost = median(Cost_physio_visit_year, na.rm = TRUE),
Total_Cost = sum(Cost_physio_visit_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of physio test
Physio_test <-40
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_Physiotest12 = if_else(
is.na(Physiotest12),
if_else(!is.na(Physiotest6), Physiotest6,
if_else(!is.na(Physiotest3), Physiotest3,0)), Physiotest12),
Cost_physio_test_year = Adjusted_Physiotest12 * Physio_test)
# Summarize the physio test costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_physio_test_year, na.rm = TRUE),
Median_Cost = median(Cost_physio_test_year, na.rm = TRUE),
Total_Cost = sum(Cost_physio_test_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of GP
GP_visit <-40
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_GP12 = if_else(
is.na(GP12),
if_else(!is.na(GP6), GP6,
if_else(!is.na(GP3), GP3,0)), GP12),
Cost_GP_year = Adjusted_GP12 * GP_visit)
# Summarize the GP costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_GP_year, na.rm = TRUE),
Median_Cost = median(Cost_GP_year, na.rm = TRUE),
Total_Cost = sum(Cost_GP_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of ER
ER_visit <-60
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_ER12 = if_else(
is.na(ER12),
if_else(!is.na(ER6), ER6,
if_else(!is.na(ER3), ER3,0)), ER12),
Cost_ER_year = Adjusted_ER12 * ER_visit)
# Summarize the ER costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_ER_year, na.rm = TRUE),
Median_Cost = median(Cost_ER_year, na.rm = TRUE),
Total_Cost = sum(Cost_ER_year, na.rm = TRUE) )
View(Exam_data_2024_summary)# Cost of outpatients
Outp_visit <-180
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_Outpatcli12 = if_else(
is.na(Outpatcli12),
if_else(!is.na(Outptcl6), Outptcl6,
if_else(!is.na(Outptcl3), Outptcl3,0)), Outpatcli12),
Cost_OUTP_year = Adjusted_Outpatcli12 * Outp_visit)
# Summarize the outpatients costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_OUTP_year, na.rm = TRUE),
Median_Cost = median(Cost_OUTP_year, na.rm = TRUE),
Total_Cost = sum(Cost_OUTP_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of home nurse visits
Nurse_visit <-40
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_HomeNurse12 = if_else(
is.na(HomeNurse12),
if_else(!is.na(HomeNurese6), HomeNurese6,
if_else(!is.na(HomeNurse3), HomeNurse3,0)), HomeNurse12),
Cost_HOME_NURSE_year = Adjusted_HomeNurse12 * Nurse_visit)
# Summarize the home nurse costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Cost_HOME_NURSE_year, na.rm = TRUE),
Median_Cost = median(Cost_HOME_NURSE_year, na.rm = TRUE),
Total_Cost = sum(Cost_HOME_NURSE_year, na.rm = TRUE) )
View(Exam_data_2024_summary) # Cost of help at home
Home_help <- 35
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_Homehelp12 = if_else(
is.na(Homehelp12),
if_else(!is.na(Homehelp6), Homehelp6,
if_else(!is.na(HomeHelp3), HomeHelp3,0)), Homehelp12),
Cost_HOME_HELP = Adjusted_Homehelp12 * Home_help )
# Summarize the home help costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize( Mean_Cost = mean(Cost_HOME_HELP, na.rm = TRUE),
Median_Cost = median(Cost_HOME_HELP, na.rm = TRUE),
Total_Cost = sum(Cost_HOME_HELP, na.rm = TRUE) )
View(Exam_data_2024_summary) Cost Societal Perspective
# Start cost SOCIETAL PERSPECTIVE
# Cost of having a day of sick leave
Sick_day <-225
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Adjusted_NOsickleave12 = if_else(
is.na(NOsickleave12),
if_else(!is.na(NOsickleav6), NOsickleav6,
if_else(!is.na(NOsickleave3), NOsickleave3,0)), NOsickleave12 ),
Cost_SICK_LEAVE_year = Adjusted_NOsickleave12 * Sick_day )
# Summarize the sick leave
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize( Mean_Cost = mean(Cost_SICK_LEAVE_year , na.rm = TRUE),
Median_Cost = median(Cost_SICK_LEAVE_year , na.rm = TRUE),
Total_Cost = sum(Cost_SICK_LEAVE_year , na.rm = TRUE) )
View(Exam_data_2024_summary)# Cost of Travelling
Travel_cost <-20
Exam_data_2024 <- Exam_data_2024 %>%
mutate(
Travel_cost_GP12 = Adjusted_GP12 * Travel_cost,
Travel_cost_Outpatcli12 = Adjusted_Outpatcli12 * Travel_cost,
Travel_cost_NOPhysio12 = Adjusted_NOphysio12 * Travel_cost,
Travel_cost_total = Travel_cost_GP12 + Travel_cost_Outpatcli12 + Travel_cost_NOPhysio12)
# Summarize the travelling costs by group
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Travel_cost_total , na.rm = TRUE),
Median_Cost = median(Travel_cost_total , na.rm = TRUE),
Total_Cost = sum(Travel_cost_total, na.rm = TRUE) )
View(Exam_data_2024_summary) Total cost healthcare perspective
#### TOTAL COST ####
Exam_data_2024 <- Exam_data_2024 %>%
mutate( Total_HealthcareP_Costs = Cost_HOME_NURSE_year + Cost_OUTP_year + Cost_ER_year + Cost_GP_year + Cost_physio_test_year + Cost_physio_visit_year + Cost_XRAY_year + Cost_CT_year + Cost_surgery_material_year + Cost_LOS_nursing_year + Cost_LOS_REHAB + Cost_ARG_12_months_reop + Cost_hospital_arg12 + Cost_hospital_YEAR + Cost_HOME_HELP)
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Total_HealthcareP_Costs , na.rm = TRUE),
Median_Cost = median(Total_HealthcareP_Costs , na.rm = TRUE),
Total_Cost = sum(Total_HealthcareP_Costs, na.rm = TRUE) )
View(Exam_data_2024_summary)Total cost Societal perspective
Exam_data_2024 <- Exam_data_2024 %>%
mutate( Total_SocietalP_Costs = Cost_HOME_NURSE_year + Cost_OUTP_year + Cost_ER_year + Cost_GP_year + Cost_physio_test_year + Cost_physio_visit_year + Cost_XRAY_year + Cost_CT_year + Cost_surgery_material_year + Cost_LOS_nursing_year + Cost_LOS_REHAB + Cost_ARG_12_months_reop + Cost_hospital_arg12 + Cost_hospital_YEAR + Cost_HOME_HELP + Cost_SICK_LEAVE_year +Travel_cost_total)
Exam_data_2024_summary <- Exam_data_2024 %>% group_by(Group) %>%
summarize(
Mean_Cost = mean(Total_SocietalP_Costs , na.rm = TRUE),
Median_Cost = median(Total_SocietalP_Costs , na.rm = TRUE),
Total_Cost = sum(Total_SocietalP_Costs, na.rm = TRUE) )
View(Exam_data_2024_summary) REGRESSIONS
Regression HEALTHCARE PERSPECTIVE
OLS
# Total_HealthcareP_Costs
# Regressions - OLS
model_OLS1 <- lm(Total_HealthcareP_Costs ~ factor(Group), data = Exam_data_2024)
summary(model_OLS1)##
## Call:
## lm(formula = Total_HealthcareP_Costs ~ factor(Group), data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10865 -8064 -4779 7124 32705
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12921 2136 6.049 2.11e-07 ***
## factor(Group)Control -1401 3021 -0.464 0.645
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10680 on 48 degrees of freedom
## Multiple R-squared: 0.004461, Adjusted R-squared: -0.01628
## F-statistic: 0.2151 on 1 and 48 DF, p-value: 0.6449# Check predictions
Exam_data_2024$ols1_Total_healthcareP_Costs <- predict(model_OLS1)
# Optional: Check model assumptions with diagnostic plots
par(mfrow = c(2, 2))
plot(model_OLS1)
linktest <- function(model_OLS1) {
pred <- predict(model_OLS1)
model2 <- lm(Total_SocietalP_Costs ~ pred + I(pred^2), data = Exam_data_2024)
return(summary(model_OLS1))
}
#linktest(model_OLS2)
resettest(model_OLS1)##
## RESET test
##
## data: model_OLS1
## RESET = 0, df1 = 2, df2 = 46, p-value = 1# Heteroscedasticity test
bptest(model_OLS1)##
## studentized Breusch-Pagan test
##
## data: model_OLS1
## BP = 1.0806, df = 1, p-value = 0.2986# HEALTHCARE INCLUDING COVARIATES
# Including covariates
model_OLS2 <- lm(Total_HealthcareP_Costs ~ factor(Group) + factor(age_group) + factor(Gender), data = Exam_data_2024)
summary(model_OLS2)##
## Call:
## lm(formula = Total_HealthcareP_Costs ~ factor(Group) + factor(age_group) +
## factor(Gender), data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10722 -6202 -2085 2736 26076
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4554 3674 1.239 0.222060
## factor(Group)Control -1530 2938 -0.521 0.605301
## factor(age_group)65-69 4452 4271 1.043 0.303139
## factor(age_group)70-74 8495 4648 1.828 0.074701 .
## factor(age_group)75-79 7332 4094 1.791 0.080551 .
## factor(age_group)80-84 15310 4954 3.090 0.003541 **
## factor(age_group)85+ 20593 4849 4.247 0.000117 ***
## factor(Gender)Male 2916 4593 0.635 0.528870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9208 on 42 degrees of freedom
## Multiple R-squared: 0.3525, Adjusted R-squared: 0.2446
## F-statistic: 3.267 on 7 and 42 DF, p-value: 0.00734# Check predictions for model 2
Exam_data_2024$ols2_Total_HealthcareP_Costs <- predict(model_OLS2)
par(mfrow = c(2, 2))
plot(model_OLS2)
# Linktest
linktest <- function(model_OLS2) {
pred <- predict(model_OLS2)
model2 <- lm(Total_HealthcareP_Costs ~ pred + I(pred^2), data = Exam_data_2024)
return(summary(model_OLS2))
}
# Linktest(model_OLS2)
resettest(model_OLS2)##
## RESET test
##
## data: model_OLS2
## RESET = 0.62332, df1 = 2, df2 = 40, p-value = 0.5413# Heteroscedasticity test
bptest(model_OLS2)##
## studentized Breusch-Pagan test
##
## data: model_OLS2
## BP = 7.6799, df = 7, p-value = 0.3617LOG
# Log transformation and regression
Exam_data_2024$lnTotal_HealthcareP_Costs <- log(Exam_data_2024$Total_HealthcareP_Costs)
model_log1 <- lm(lnTotal_HealthcareP_Costs ~ factor(Group) + factor(age_group) + factor(Gender), data = Exam_data_2024)
summary(model_log1)##
## Call:
## lm(formula = lnTotal_HealthcareP_Costs ~ factor(Group) + factor(age_group) +
## factor(Gender), data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.42265 -0.50108 0.02054 0.38321 1.78330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.3784 0.2963 28.276 < 2e-16 ***
## factor(Group)Control -0.4711 0.2370 -1.988 0.053372 .
## factor(age_group)65-69 0.2833 0.3444 0.822 0.415513
## factor(age_group)70-74 0.8440 0.3748 2.252 0.029630 *
## factor(age_group)75-79 0.9753 0.3302 2.954 0.005126 **
## factor(age_group)80-84 1.5671 0.3996 3.922 0.000319 ***
## factor(age_group)85+ 2.0292 0.3911 5.189 5.75e-06 ***
## factor(Gender)Male 0.2550 0.3704 0.688 0.494974
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7427 on 42 degrees of freedom
## Multiple R-squared: 0.5009, Adjusted R-squared: 0.4178
## F-statistic: 6.023 on 7 and 42 DF, p-value: 6.77e-05# Predicted log cost
Exam_data_2024$log_Total_HealthcareP_Costs <- predict(model_log1)
# Back-transforming log predictions
Exam_data_2024$log_Total_HealthcareP_Costs_raw <- exp(Exam_data_2024$log_Total_HealthcareP_Costs)
Exam_data_2024$smr <- exp(Exam_data_2024$lnTotal_HealthcareP_Costs - Exam_data_2024$log_Total_HealthcareP_Costs)
# Mean smearing factor
smear <- mean(Exam_data_2024$smr, na.rm = TRUE)
Exam_data_2024$mu <- exp(Exam_data_2024$log_Total_HealthcareP_Costs) * smear
# Comparing adjusted costs
Exam_data_2024 %>%
summarize(adjusted_vs_observed = list(c(mu = mean(mu), Total_HealthcareP_Costs = mean(Total_HealthcareP_Costs))))## # A tibble: 1 × 1
## adjusted_vs_observed
## <list>
## 1 <dbl [2]># Re-run regression with log costs
model_log2 <- lm(lnTotal_HealthcareP_Costs ~ factor(Group) + factor(age_group) + factor(Gender), data = Exam_data_2024)
summary(model_log2)##
## Call:
## lm(formula = lnTotal_HealthcareP_Costs ~ factor(Group) + factor(age_group) +
## factor(Gender), data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.42265 -0.50108 0.02054 0.38321 1.78330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.3784 0.2963 28.276 < 2e-16 ***
## factor(Group)Control -0.4711 0.2370 -1.988 0.053372 .
## factor(age_group)65-69 0.2833 0.3444 0.822 0.415513
## factor(age_group)70-74 0.8440 0.3748 2.252 0.029630 *
## factor(age_group)75-79 0.9753 0.3302 2.954 0.005126 **
## factor(age_group)80-84 1.5671 0.3996 3.922 0.000319 ***
## factor(age_group)85+ 2.0292 0.3911 5.189 5.75e-06 ***
## factor(Gender)Male 0.2550 0.3704 0.688 0.494974
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7427 on 42 degrees of freedom
## Multiple R-squared: 0.5009, Adjusted R-squared: 0.4178
## F-statistic: 6.023 on 7 and 42 DF, p-value: 6.77e-05par(mfrow = c(2, 2))
plot(model_log2)
# Linktest
linktest <- function(model_log2) {
pred <- predict(model_log2)
model3 <- lm(Total_HealthcareP_Costs ~ pred + I(pred^2), data = Exam_data_2024)
return(summary(model_log2))
}
#linktest(model_log2)
resettest(model_log2)##
## RESET test
##
## data: model_log2
## RESET = 3.2251, df1 = 2, df2 = 40, p-value = 0.05029# Heteroscedasticity test
bptest(model_log2)##
## studentized Breusch-Pagan test
##
## data: model_log2
## BP = 5.9296, df = 7, p-value = 0.548GLM
# Fit the generalized linear model
model_gamma_log <- glm(Total_HealthcareP_Costs ~ factor(Group) + factor(age_group) + factor(Gender), family = Gamma(link = "log"), data = Exam_data_2024)
summary(model_gamma_log)##
## Call:
## glm(formula = Total_HealthcareP_Costs ~ factor(Group) + factor(age_group) +
## factor(Gender), family = Gamma(link = "log"), data = Exam_data_2024)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.4453 0.3256 25.936 < 2e-16 ***
## factor(Group)Control -0.3753 0.2604 -1.441 0.15696
## factor(age_group)65-69 0.5606 0.3785 1.481 0.14606
## factor(age_group)70-74 1.1740 0.4119 2.850 0.00675 **
## factor(age_group)75-79 1.0817 0.3629 2.981 0.00476 **
## factor(age_group)80-84 1.5329 0.4391 3.491 0.00114 **
## factor(age_group)85+ 1.9785 0.4297 4.604 3.81e-05 ***
## factor(Gender)Male 0.3582 0.4071 0.880 0.38387
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.6660759)
##
## Null deviance: 40.532 on 49 degrees of freedom
## Residual deviance: 24.312 on 42 degrees of freedom
## AIC: 1028.4
##
## Number of Fisher Scoring iterations: 11# Get the linear predictors
Exam_data_2024$glmTotal_HealthcareP_Costs <- predict(model_gamma_log)
# Exponentiate the linear predictors to get the predicted values on the original scale
Exam_data_2024$glmTotal_HealthcareP_Costs_raw <- exp(Exam_data_2024$glmTotal_HealthcareP_Costs)
# Summarize glmTot_raw and TotalCostEURO
summary_results_glm <- Exam_data_2024 %>%
summarize(
mean_glmTotal_HealthcareP_Costs_raw = mean(glmTotal_HealthcareP_Costs_raw, na.rm = TRUE),
mean_Total_HealthcareP_Costs = mean(Total_HealthcareP_Costs, na.rm = TRUE),
sd_glmTotal_HealthcareP_Costs_raw = sd(glmTotal_HealthcareP_Costs_raw, na.rm = TRUE),
sd_Total_HealthcareP_Costs = sd(Total_HealthcareP_Costs, na.rm = TRUE)
)
print(summary_results_glm)## # A tibble: 1 × 4
## mean_glmTotal_HealthcareP_Cost…¹ mean_Total_Healthcar…² sd_glmTotal_Healthca…³
## <dbl> <dbl> <dbl>
## 1 12466. 12221. 7091.
## # ℹ abbreviated names: ¹mean_glmTotal_HealthcareP_Costs_raw,
## # ²mean_Total_HealthcareP_Costs, ³sd_glmTotal_HealthcareP_Costs_raw
## # ℹ 1 more variable: sd_Total_HealthcareP_Costs <dbl>linktest <- function(model_gamma_log) {
pred <- predict(model_gamma_log)
model3 <- lm(Total_HealthcareP_Costs ~ pred + I(pred^2), data = Exam_data_2024)
return(summary(model_gamma_log))
}
#linktest(model_gamma_log)
resettest(model_gamma_log)##
## RESET test
##
## data: model_gamma_log
## RESET = 0.62332, df1 = 2, df2 = 40, p-value = 0.5413# Heteroscedasticity test
bptest(model_gamma_log)##
## studentized Breusch-Pagan test
##
## data: model_gamma_log
## BP = 7.6799, df = 7, p-value = 0.3617Questionnare (15D, Time 0, OLS + summary)
model_OLS1_15D <- lm(Time0 ~ factor(Group), data = Exam_data_2024)
summary(model_OLS1_15D)##
## Call:
## lm(formula = Time0 ~ factor(Group), data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25955 -0.05603 0.02566 0.07417 0.11565
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.87003 0.01808 48.115 <2e-16 ***
## factor(Group)Control -0.04028 0.02557 -1.575 0.122
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09041 on 48 degrees of freedom
## Multiple R-squared: 0.04915, Adjusted R-squared: 0.02934
## F-statistic: 2.481 on 1 and 48 DF, p-value: 0.1218summary_15D_baseline <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
Mobilityb_mean = mean(Mobilityb, na.rm = TRUE),
Mobilityb_sd = sd(Mobilityb, na.rm = TRUE),
Visionb_mean = mean(Visionb, na.rm = TRUE),
Visionb_sd = sd(Visionb, na.rm = TRUE),
Hearingb_mean = mean(Hearingb, na.rm = TRUE),
Hearingb_sd = sd(Hearingb, na.rm = TRUE),
Breathingb_mean = mean(Breathingb, na.rm = TRUE),
Breathingb_sd = sd(Breathingb, na.rm = TRUE),
Sleepingb_mean = mean(Sleepingb, na.rm = TRUE),
Sleepingb_sd = sd(Sleepingb, na.rm = TRUE),
Eatingb_mean = mean(Eatingb, na.rm = TRUE),
Eatingb_sd = sd(Eatingb, na.rm = TRUE),
Speechb_mean = mean(Speechb, na.rm = TRUE),
Speechb_sd = sd(Speechb, na.rm = TRUE),
Excretionb_mean = mean(Excretionb, na.rm = TRUE),
Excretionb_sd = sd(Excretionb, na.rm = TRUE),
Usualactivitiesb_mean = mean(Usualactivitiesb, na.rm = TRUE),
Usualactivitiesb_sd = sd(Usualactivitiesb, na.rm = TRUE),
Mentalb_mean = mean(Mentalb, na.rm = TRUE),
Mentalb_sd = sd(Mentalb, na.rm = TRUE),
DiscomfortandSymptomsb_mean = mean(DiscomfortandSymptomsb, na.rm = TRUE),
DiscomfortandSymptomsb_sd = sd(DiscomfortandSymptomsb, na.rm = TRUE),
Depression1_mean = mean(Depression1, na.rm = TRUE),
Depression1_sd = sd(Depression1, na.rm = TRUE),
Distressb_mean = mean(Distressb, na.rm = TRUE),
Distressb_sd = sd(Distressb, na.rm = TRUE),
Vitalityb_mean = mean(Vitalityb, na.rm = TRUE),
Vitalityb_sd = sd(Vitalityb, na.rm = TRUE)
)
summary_15D_baseline## # A tibble: 2 × 29
## Group Mobilityb_mean Mobilityb_sd Visionb_mean Visionb_sd Hearingb_mean
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Surgery 1.28 0.614 1.08 0.277 1.28
## 2 Control 1.44 0.870 1.12 0.332 1.16
## # ℹ 23 more variables: Hearingb_sd <dbl>, Breathingb_mean <dbl>,
## # Breathingb_sd <dbl>, Sleepingb_mean <dbl>, Sleepingb_sd <dbl>,
## # Eatingb_mean <dbl>, Eatingb_sd <dbl>, Speechb_mean <dbl>, Speechb_sd <dbl>,
## # Excretionb_mean <dbl>, Excretionb_sd <dbl>, Usualactivitiesb_mean <dbl>,
## # Usualactivitiesb_sd <dbl>, Mentalb_mean <dbl>, Mentalb_sd <dbl>,
## # DiscomfortandSymptomsb_mean <dbl>, DiscomfortandSymptomsb_sd <dbl>,
## # Depression1_mean <dbl>, Depression1_sd <dbl>, Distressb_mean <dbl>, … QALYs
# Calculation for the QALYs without imputations
# QALY Calculation
Exam_data_2024$QALYs <- ((Exam_data_2024$Time0 + Exam_data_2024$Time3) / 2) * (3 / 12) + # From baseline to 3 months
((Exam_data_2024$Time3 + Exam_data_2024$Time6) / 2) * (3 / 12) + # From 3 to 6 months
((Exam_data_2024$Time6 + Exam_data_2024$Time12) / 2) * (6 / 12) # From 6 to 12 months
# Summary
library(dplyr)
# Summary of Time0 (Baseline)
Exam_data_2024_summary_Time0 <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(total_Time0 = mean(Time0, na.rm = TRUE))
print(Exam_data_2024_summary_Time0)## # A tibble: 2 × 2
## Group total_Time0
## <fct> <dbl>
## 1 Surgery 0.870
## 2 Control 0.830# Summary of Time3 (3 Months)
Exam_data_2024_summary_Time3 <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
total_Time3 = mean(Time3, na.rm = TRUE))
print(Exam_data_2024_summary_Time3)## # A tibble: 2 × 2
## Group total_Time3
## <fct> <dbl>
## 1 Surgery 0.814
## 2 Control 0.817# Summary of Time6 (6 Months)
Exam_data_2024_summary_Time6 <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(total_Time6 = mean(Time6, na.rm = TRUE))
print(Exam_data_2024_summary_Time6)## # A tibble: 2 × 2
## Group total_Time6
## <fct> <dbl>
## 1 Surgery 0.836
## 2 Control 0.817# Summary of Time12 (12 Months)
Exam_data_2024_summary_Time12 <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(total_Time12 = mean(Time12, na.rm = TRUE))
print(Exam_data_2024_summary_Time12)## # A tibble: 2 × 2
## Group total_Time12
## <fct> <dbl>
## 1 Surgery 0.841
## 2 Control 0.819# Summary of QALYs
Exam_data_2024_summary_QALYs <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
total_QALYs = mean(QALYs, na.rm = TRUE),
sd_QALYs = sd(QALYs, na.rm = TRUE)
)
print(Exam_data_2024_summary_QALYs)## # A tibble: 2 × 3
## Group total_QALYs sd_QALYs
## <fct> <dbl> <dbl>
## 1 Surgery 0.837 0.0738
## 2 Control 0.819 0.0762# Summarize QALYs by group
library(dplyr)
Exam_data_2024_summary_QALYs <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
total_QALYs = mean(QALYs, na.rm = TRUE),
sd_QALYs = sd(QALYs, na.rm = TRUE)
)
print(Exam_data_2024_summary_QALYs)## # A tibble: 2 × 3
## Group total_QALYs sd_QALYs
## <fct> <dbl> <dbl>
## 1 Surgery 0.837 0.0738
## 2 Control 0.819 0.0762# Manca method (et al) : adjusting for baseline differences
library(broom)
model_manca <- lm(QALYs ~ factor(Group) + Time0, data = Exam_data_2024, na.action = na.exclude)
summary(model_manca)##
## Call:
## lm(formula = QALYs ~ factor(Group) + Time0, data = Exam_data_2024,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.119229 -0.033020 0.004049 0.027713 0.130251
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.341260 0.080607 4.234 0.000112 ***
## factor(Group)Control 0.008794 0.016647 0.528 0.599919
## Time0 0.565326 0.090986 6.213 1.5e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05569 on 45 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.4695, Adjusted R-squared: 0.4459
## F-statistic: 19.91 on 2 and 45 DF, p-value: 6.394e-07tidy_results_manca <- broom::tidy(model_manca)
print(tidy_results_manca)## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.341 0.0806 4.23 0.000112
## 2 factor(Group)Control 0.00879 0.0166 0.528 0.600
## 3 Time0 0.565 0.0910 6.21 0.000000150# Imputation
# Missing information in several of the observation periods
# Total sample 50
# Intervention = Surgery = 1
# Control = Conservative = 2
# Surgery n = 25 and prothesis n = 25
# Missing values for 15D
# predicting missing variables - at 3 months for Time 3
# First regression model
model_Time3 <- lm(Time3 ~ factor(Group) + factor(Gender) + Time0, data = Exam_data_2024)
summary(model_Time3)##
## Call:
## lm(formula = Time3 ~ factor(Group) + factor(Gender) + Time0,
## data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.128627 -0.046584 -0.005661 0.044944 0.142213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.51745 0.08622 6.002 3.1e-07 ***
## factor(Group)Control 0.02069 0.01845 1.121 0.26808
## factor(Gender)Male 0.02486 0.02778 0.895 0.37562
## Time0 0.33446 0.09845 3.397 0.00143 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06122 on 45 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.2269, Adjusted R-squared: 0.1754
## F-statistic: 4.403 on 3 and 45 DF, p-value: 0.008458# Predict Time3 values
Exam_data_2024$Time3_p<- predict(model_Time3, newdata = Exam_data_2024)
summary(Exam_data_2024$Time3_p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.7289 0.7949 0.8208 0.8150 0.8389 0.8585# Create new variable Time_3_all and replace missing values
Exam_data_2024$Time3_all <- Exam_data_2024$Time3
Exam_data_2024$Time3_all[is.na(Exam_data_2024$Time3_all)] <-
Exam_data_2024$Time3_p[is.na(Exam_data_2024$Time3_all)]
# Summarize Time3_all by Group
summary_Time3_all <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(mean_Time3_all = mean(Time3_all, na.rm = TRUE),
sd_Time3_all = sd(Time3_all, na.rm = TRUE),
n = n())
print(summary_Time3_all)## # A tibble: 2 × 4
## Group mean_Time3_all sd_Time3_all n
## <fct> <dbl> <dbl> <int>
## 1 Surgery 0.813 0.0598 25
## 2 Control 0.817 0.0743 25# Summarize Time_3 by Group again (original values)# predicting missing variables - at 6 months for Time 6
# First regression model
model_Time6 <- lm(Time6 ~ factor(Group) + factor(Gender) + Time0, data = Exam_data_2024)
summary(model_Time6)##
## Call:
## lm(formula = Time6 ~ factor(Group) + factor(Gender) + Time0,
## data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.142281 -0.050581 0.004145 0.043291 0.172674
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.36260 0.10262 3.533 0.000978 ***
## factor(Group)Control 0.01279 0.02182 0.586 0.560606
## factor(Gender)Male 0.03740 0.03225 1.160 0.252394
## Time0 0.53057 0.11631 4.562 4.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07087 on 44 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3548, Adjusted R-squared: 0.3108
## F-statistic: 8.065 on 3 and 44 DF, p-value: 0.0002167# Predict Time6 values
Exam_data_2024$Time6_p<- predict(model_Time6, newdata = Exam_data_2024)
summary(Exam_data_2024$Time6_p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6779 0.7958 0.8371 0.8244 0.8642 0.9016# Create new variable Time_6_all and replace missing values
Exam_data_2024$Time6_all <- Exam_data_2024$Time6
Exam_data_2024$Time6_all[is.na(Exam_data_2024$Time6_all)] <-
Exam_data_2024$Time6_p[is.na(Exam_data_2024$Time6_all)]
# Summarize Time6_all by Group
summary_Time6_all <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(mean_Time6_all = mean(Time6_all, na.rm = TRUE),
sd_Time6_all = sd(Time6_all, na.rm = TRUE),
n = n())
print(summary_Time6_all)## # A tibble: 2 × 4
## Group mean_Time6_all sd_Time6_all n
## <fct> <dbl> <dbl> <int>
## 1 Surgery 0.832 0.0812 25
## 2 Control 0.817 0.0882 25# Summarize Time_6 by Group again (original values)# predicting missing variables - at 12 months for Time 12
# First regression model
model_Time12 <- lm(Time12 ~ factor(Group) + factor(Gender) + Time0, data = Exam_data_2024)
summary(model_Time12)##
## Call:
## lm(formula = Time12 ~ factor(Group) + factor(Gender) + Time0,
## data = Exam_data_2024)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.215210 -0.042438 0.004452 0.057386 0.165061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.32905 0.11063 2.974 0.00475 **
## factor(Group)Control 0.01660 0.02352 0.706 0.48411
## factor(Gender)Male 0.06528 0.03477 1.878 0.06704 .
## Time0 0.56770 0.12539 4.527 4.51e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07641 on 44 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.3786, Adjusted R-squared: 0.3363
## F-statistic: 8.937 on 3 and 44 DF, p-value: 9.737e-05# Predict Time6 values
Exam_data_2024$Time12_p<- predict(model_Time12, newdata = Exam_data_2024)
summary(Exam_data_2024$Time12_p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6694 0.7948 0.8392 0.8277 0.8658 0.9310# Create new variable Time_12_all and replace missing values
Exam_data_2024$Time12_all <- Exam_data_2024$Time12
Exam_data_2024$Time12_all[is.na(Exam_data_2024$Time12_all)] <-
Exam_data_2024$Time12_p[is.na(Exam_data_2024$Time12_all)]
# Summarize Time12_all by Group
summary_Time12_all <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(mean_Time12_all = mean(Time12_all, na.rm = TRUE),
sd_Time12_all = sd(Time12_all, na.rm = TRUE),
n = n())
print(summary_Time12_all)## # A tibble: 2 × 4
## Group mean_Time12_all sd_Time12_all n
## <fct> <dbl> <dbl> <int>
## 1 Surgery 0.836 0.103 25
## 2 Control 0.819 0.0826 25# Summarize Time_12 by Group again (original values)#Estimation of QALYs after imputation
Exam_data_2024$QALYs_all <- ((Exam_data_2024$Time0 + Exam_data_2024$Time3_all) / 2 )* (3 / 12) +
((Exam_data_2024$Time3_all + Exam_data_2024$Time6_all) / 2) * (3 / 12) +
((Exam_data_2024$Time6_all + Exam_data_2024$Time12_all) / 2) * (6 /12)
Exam_data_2024_summary_QALYs_all <- Exam_data_2024 %>%
group_by(Group) %>%
summarize(
total_QALYs_all = mean(QALYs_all, na.rm = TRUE),
sd_QALYs = sd(QALYs_all, na.rm = TRUE)
)
print(Exam_data_2024_summary_QALYs_all)## # A tibble: 2 × 3
## Group total_QALYs_all sd_QALYs
## <fct> <dbl> <dbl>
## 1 Surgery 0.833 0.0722
## 2 Control 0.819 0.0762# Adjusting for differences at baseline - Using method suggested by Manca et al
model_manca_all <- lm(QALYs_all ~ factor(Group) + Time0, data = Exam_data_2024, na.action = na.exclude)
summary(model_manca_all)##
## Call:
## lm(formula = QALYs_all ~ factor(Group) + Time0, data = Exam_data_2024,
## na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.119127 -0.030580 0.004092 0.026470 0.130154
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.341909 0.076484 4.470 4.91e-05 ***
## factor(Group)Control 0.008871 0.015809 0.561 0.577
## Time0 0.564452 0.087012 6.487 4.91e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0545 on 47 degrees of freedom
## Multiple R-squared: 0.4771, Adjusted R-squared: 0.4549
## F-statistic: 21.45 on 2 and 47 DF, p-value: 2.41e-07tidy_results_manca_all <- broom::tidy(model_manca_all)
print(tidy_results_manca_all)## # A tibble: 3 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 0.342 0.0765 4.47 0.0000491
## 2 factor(Group)Control 0.00887 0.0158 0.561 0.577
## 3 Time0 0.564 0.0870 6.49 0.0000000491 DETERMINISTIC ICER
# incremental costs / incremental QALYs
# ICER = INTERVENTION – CONTROL
# ICER = SURGERY – CONSERVATIVE
# ICER formula =
# (Cost intervention A - Cost intervention B) /
# (Effectiveness Intervention A - Effectiveness intervention B)
# HEALTHCARE
# COST
# Healthcare total cost mean surgery = 12921.2
# Healthcare total cost mean control = 11520.2
# 12921.2 - 11520.2 = 1,401
INC_COST_HP= 12921.2 - 11520.2
print(INC_COST_HP)## [1] 1401 # QALY
# QALY surgery = 0.8368793
# QALY control = 0.8191325
# 0.8368793 - 0.8191325 = 0.0177468
INC_QALY_HP= 0.8368793 - 0.8191325
print(INC_QALY_HP)## [1] 0.0177468 # ICER HEALTHCARE
ICER_HP= 1401 / 0.0177414
print(ICER_HP) ## [1] 78967.84 # 1401 / 0.0177414 = 78967.84 = ICER HEALTHCARE PERSPECTIVE# SOCIETY ICER CALCULATION
# Societal total cost mean surgery = 18735.8
# Societal total cost mean control = 16108.4
# 18735.8 - 16108.4 = 2,627.4
INC_COST_SP= 18735.8 - 16108.4
print(INC_COST_SP)## [1] 2627.4 # QALY
# QALY surgery = 0.8368793
# QALY control = 0.8191325
# 0.8368793 - 0.8191325 = 0.0177414
INC_QALY_SP= 0.8368793 - 0.8191325
print(INC_QALY_SP)## [1] 0.0177468 # ICER SOCIETY
ICER_SP= 2627.4 / 0.0177468
print(ICER_SP) ## [1] 148049.2 # 2,627.4 / 0.0177468 = 148049.2 = ICER SOCIETAL PERSPECTIVE UNCERTAINTY
BOOTSTRAP + AUC
# Bootstrap results of costs and QALYs
# Hemiprothesis = intervention = surgery
# Internal Fixation (IF) = control = conservative
# CostGroup1 = Total_HealthcareP_Costs
# Bootstrap for societal cost
# Bootstrap for Surgery (Qalys and intervention)
# Bootstrap for Conservative (Qalys and intervention)BOOTSTRAP HEALTHCARE PERSPECTIVE
# Bootstrap results of costs and QALYs
# HEALTHCARE PERSPECTIVE
# Create Total_HealthcareP_Costs_Surgery for Group == Intervention
Exam_data_2024$Total_HP_Costs_Surgery <- ifelse(Exam_data_2024$Group == "Surgery", Exam_data_2024$Total_HealthcareP_Costs,NA)
summary(Exam_data_2024$Total_HP_Costs_Surgery)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 4195 5775 9475 12921 19775 37985 25# Create totalcost_IF and CostGroup1_IF for Group == Control
Exam_data_2024$Total_HP_Costs_Control <- ifelse(Exam_data_2024$Group == "Control", Exam_data_2024$Total_HealthcareP_Costs,NA)
summary(Exam_data_2024$Total_HP_Costs_Control)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 655 2925 4320 11520 20495 44225 25# Define QALYs/AUC for the two interventions
Exam_data_2024$QALYs_HP_Surgery <- ifelse(Exam_data_2024$Group == "Surgery", Exam_data_2024$QALYs, NA)
summary(Exam_data_2024$QALYs_HP_Surgery)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.6598 0.8165 0.8610 0.8369 0.8929 0.9397 27Exam_data_2024$QALYs_HP_Control <- ifelse(Exam_data_2024$Group == "Control", Exam_data_2024$QALYs, NA)
summary(Exam_data_2024$QALYs_HP_Control)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.6284 0.7786 0.8104 0.8191 0.8869 0.9434 25# Bootstrap for Surgery
mean_fun <- function(data, indices) {
d <- data[indices, ] # Resample data with indices
return(c(mean(d$Total_HP_Costs_Surgery, na.rm = TRUE),
mean(d$QALYs_HP_Surgery, na.rm = TRUE)))
}
# Set the number of bootstrap replications
n_rep <- 1000
# Perform the bootstrap
boot_results_HP_Surgery <- boot(data = Exam_data_2024, statistic = mean_fun, R = n_rep)
# Extract bootstrap replicates (t) and convert it into a data frame
bootstrap_df_HP_Surgery <- as.data.frame(boot_results_HP_Surgery$t)
# Rename the columns for clarity
colnames(bootstrap_df_HP_Surgery) <- c("Total_HP_Costs_Surgery", "QALYs_HP_Surgery")
# Bootstrap for Control
mean_fun <- function(data, indices) {
d <- data[indices, ] # Resample data with indices
return(c(mean(d$Total_HP_Costs_Control, na.rm = TRUE),
mean(d$QALYs_HP_Control, na.rm = TRUE)))
}
# Set the number of bootstrap replications
n_rep <- 1000
# Perform the bootstrap
boot_results_HP_Control <- boot(data = Exam_data_2024, statistic = mean_fun, R = n_rep)
# Extract bootstrap replicates (t) and convert it into a data frame
bootstrap_df_HP_Control <- as.data.frame(boot_results_HP_Control$t)
# Rename the columns for clarity
colnames(bootstrap_df_HP_Control) <- c("Total_HP_Costs_Control", "QALYs_HP_Control")
##
write.xlsx(bootstrap_df_HP_Surgery, "/Users/ruthe/OneDrive/documents")
write.xlsx(bootstrap_df_HP_Control, "/Users/ruthe/OneDrive/documents")BOOTSTRAP SOCIETAL PERSPECTIVE
# Bootstrap results of costs and QALYs
# SOCIETAL PERSPECTIVE
# Create Total_SocietalP_Costs_Surgery for Group == Intervention
Exam_data_2024$Total_SP_Costs_Surgery <- ifelse(Exam_data_2024$Group == "Surgery", Exam_data_2024$Total_SocietalP_Costs,NA)
summary(Exam_data_2024$Total_SP_Costs_Surgery)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 4455 7055 12940 18736 22025 74820 25# Create totalcost_IF and CostGroup1_IF for Group == Control
Exam_data_2024$Total_SP_Costs_Control <- ifelse(Exam_data_2024$Group == "Control", Exam_data_2024$Total_SocietalP_Costs,NA)
summary(Exam_data_2024$Total_SP_Costs_Control)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1370 4835 16175 16108 22130 45025 25# Define QALYs/AUC for the two interventions
Exam_data_2024$QALYs_SP_Surgery <- ifelse(Exam_data_2024$Group == "Surgery", Exam_data_2024$QALYs, NA)
summary(Exam_data_2024$QALYs_SP_Surgery)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.6598 0.8165 0.8610 0.8369 0.8929 0.9397 27Exam_data_2024$QALYs_SP_Control <- ifelse(Exam_data_2024$Group == "Control", Exam_data_2024$QALYs, NA)
summary(Exam_data_2024$QALYs_SP_Control)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.6284 0.7786 0.8104 0.8191 0.8869 0.9434 25# Bootstrap for Surgery
mean_fun <- function(data, indices) {
d <- data[indices, ] # Resample data with indices
return(c(mean(d$Total_SP_Costs_Surgery, na.rm = TRUE),
mean(d$QALYs_SP_Surgery, na.rm = TRUE)))
}
# Set the number of bootstrap replications
n_rep <- 1000
# Perform the bootstrap
boot_results_SP_Surgery <- boot(data = Exam_data_2024, statistic = mean_fun, R = n_rep)
# Extract bootstrap replicates (t) and convert it into a data frame
bootstrap_df_SP_Surgery <- as.data.frame(boot_results_SP_Surgery$t)
# Rename the columns for clarity
colnames(bootstrap_df_SP_Surgery) <- c("Total_SP_Costs_Surgery", "QALYs_SP_Surgery")
# Bootstrap for Control
mean_fun <- function(data, indices) {
d <- data[indices, ] # Resample data with indices
return(c(mean(d$Total_SP_Costs_Control, na.rm = TRUE),
mean(d$QALYs_SP_Control, na.rm = TRUE)))
}
# Set the number of bootstrap replications
n_rep <- 1000
# Perform the bootstrap
boot_results_SP_Control <- boot(data = Exam_data_2024, statistic = mean_fun, R = n_rep)
# Extract bootstrap replicates (t) and convert it into a data frame
bootstrap_df_SP_Control <- as.data.frame(boot_results_SP_Control$t)
# Rename the columns for clarity
colnames(bootstrap_df_SP_Control) <- c("Total_SP_Costs_Control", "QALYs_SP_Control")
#####
write.xlsx(bootstrap_df_SP_Surgery, "/Users/ruthe/OneDrive/documents")
write.xlsx(bootstrap_df_SP_Control, "/Users/ruthe/OneDrive/documents")