What does this code do?

packages <- c("haven", "car", "broom", "Matching", "lmtest", "ggplot2",
              "sandwich", "AER")

#install any missing packages & load all packages
for (pkg in packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg)
      }
  library(pkg, character.only = TRUE)
  }

##set working directory
setwd("//cmvm.datastore.ed.ac.uk/cmvm/smgphs/users/gtramon2/Win7/Desktop/For Students/05-RWE materials")

#load data
data <- read.csv("RWE_data_2025.csv")

# COSTS
#note: "actual" treatment effect is 50 EUR!

# QUESTION 1 - summarise costs
summary(data$costs)





#some negative costs, =-10! need to be careful how we deal with these...
#this might be a missing data code
#this might say something about costs
#this might be a systematic typo + should be 10
#need to know more about our data, really...

#proceed under assumption here that -10 should just be zero
# NOTE: THIS WILL NOT NECESSARILY BE VALID


#replace negative costs with zero
data$costs[data$costs < 0] <- 0

#summarise costs again
summary(data$costs)

# QUESTION 2
#Regression with costs on height
summary(lm(costs ~ height, data = data))

# although this is a statistically significant relationship, we should think about potential 
# reasons for this:

# an artefact of using a large sample
# relationship between gender + outcomes
### e.g. check distributions of height by gender


#summarise costs by gender
summary(data$costs[data$male == 0])
summary(data$costs[data$male == 1])

#t-test of costs by gender
#again remember R very annoying here
t.test(data$costs ~ data$male)

#histogram of costs by gender
ggplot(data, aes(x = costs)) +
  geom_histogram(aes(fill = factor(male)), position = "identity", alpha = 0.5) +
  theme_minimal()

# QUESTION 3

#summarise costs by treatment & compare means
summary(data$costs[data$treat == 0])
summary(data$costs[data$treat == 1])

#first OLS Regression with costs on treatment
summary(lm(costs ~ treat, data = data))

#unlikely to be causal given what we know (& what we're about to show)

# QUESTION 4
summary(lm(costs ~ treat + male, data = data))

# note fall in estimated relationship here

# QUESTION 5

#We know that doctor's beliefs (although unobservable) affect the individual's
#frailty score and suspect they simultaneously affect both probability of 
#treatment and costs incurred - due to increased follow-up care.
#Frailty is therefore potentially an instrument for treatment.

#use ivreg command

summary(ivreg(costs ~ treat + male | frailty + male, data = data))

#Manually constructing IV 
#probit regression for treatment on frailty and gender
logit_model <- glm(treat ~ frailty + male, data = data, family = binomial(link = "logit"))
data$prob <- predict(logit_model, type = "response")
summary(lm(costs ~ prob + male, data = data))

#this is our causal effect (perhaps!) - approx 49 EUR (note logit not LPM)

# QUESTION 6

#regression of costs adjusted for treatment, male, and number of appointments
summary(lm(costs ~ treat + male + n_appts, data = data))

# QUESTION 7

#treatment Effects with nearest neighbour matching
#perform nearest neighbor matching using 'Match' from the 'Matching' package.
#Specifies Y outcome variable of costs, 'Tr' is the treatment indicator (treat),
#& 'X' contains the covariates to match on (male, n_appts).
#'M=1' specifies one-to-one matching, and 'replace=TRUE' allows for replacement.
m_out <- Match(Y = data$costs, Tr = data$treat,
               X = cbind(data$male, data$n_appts),
               M = 1, replace = TRUE)

# Estimate the average treatment effect on the treated (ATT)
r_out <- summary(m_out)

# OUTCOMES

# QUESTION 1

## First we need to generate QALYs from these. Area under curve...
data$qalys <- (data$eq5d_baseline + data$eq5d_final)/2

#regression of QALYs with treatment as only covariate
summary(lm(qalys ~ treat, data = data))

# QUESTION 2
#regression of QALYs with treatment and male as covariates
summary(lm(qalys ~ treat + male, data = data))

# QUESTION 3

# iv regression with QALYs on treatment and male
data$prob_out <- predict(logit_model, type = "response")
summary(lm(qalys ~ prob_out + male, data = data))

# deciles for frailty??

#cut by deciles
data$fraildec <- cut(data$frailty, breaks = 10, labels = FALSE)

#prepare blank df to store the coefficients and standard errors
treat_effects <- data.frame(decile = integer(0), coefficient = numeric(0), std_error = numeric(0))

#run regression by frailty decile and store results
for (i in 1:10) {
  model <- lm(qalys ~ treat + male, data = data, subset = (fraildec == i))
  coef_summary <- summary(model)$coefficients
  
  #as 'treat' is second coefficient in model
  treat_coef <- coef_summary["treat", "Estimate"]
  treat_se <- coef_summary["treat", "Std. Error"]  

#add the results to the data frame
  treat_effects <- rbind(treat_effects, c(i, treat_coef, treat_se))
}

#rename columns for clarity
names(treat_effects) <- c("Decile", "Coefficient", "Standard Error")

#output the summary of treatment effects by decile
print(treat_effects)
##these look pretty similar

#many options for testing here but for simplicity do decile 1 vs decile 10
coef1 <- treat_effects$Coefficient[treat_effects$Decile == 1]
se1 <- treat_effects$`Standard Error`[treat_effects$Decile == 1]
coef10 <- treat_effects$Coefficient[treat_effects$Decile == 10]
se10 <- treat_effects$`Standard Error`[treat_effects$Decile == 10]

#calculate the t-statistic & p-value
t_statistic <- (coef1 - coef10) / sqrt(se1^2 + se10^2)
p_value <- 2 * pnorm(-abs(t_statistic))

p_value
#not even close to significant at conventional levels

# //question 4
# /*
#   These estimates do not differ much - seems like this selection process isn't
# affecting estimates much (possibly as we might suspect). 
# 
# Frailty might directly affect outcomes to extent that it is usefully recorded.
# This was not the case for costs by assumption.
# 
# Perhaps suggests extra care is redundant in improving outcomes?
# 
# IVs reduce efficiency so bigger CIs.
# */
# 
# //question 5
# /*
# Potentially, yes. If we believe our IV estimates are preferable to our OLS
# estimates for costs then we'd select this. The decision regarding our best 
# estimates for outcomes would appear to be more straightforward: we both have
# less reason to think that an IV estimate is likely to be preferable and also
# our results do not differ much (compare CIs for different estimates).
# */
#   
#   //COST-EFFECTIVENESS
# *for discussion!