Assessing wellness programme

Executive Summary

This report investigates whether a workplace wellness program reduces employee medical costs. The analysis utilizes a Randomized Control Trial (RCT), considered the gold standard for causal inference.

The Business Question: Does offering a wellness program reduce medical spending?

Data Dictionary: Understanding the Variables

Key variables used in the analysis include:

• Random Assignment: 1 = Offered the program, 0 = Not offered.
• Participation Status: 1 = Enrolled, 0 = Did not enroll.
• Insurance Coverage Weight: Used to adjust for partial-year employees.
• Baseline Total Spending: Medical costs before the program (2015-2016).
• Outcome Total Spending: Medical costs during the program (2016-2017).
• Baseline Hospital Spending: Hospital costs (2015-2016).
• Baseline Office Visit Spending: Office visit costs (2015-2016).
• Outcome Hospital Spending: Hospital costs (2016-2017).
• Outcome Office Visit Spending: Office visit costs (2016-2017).
• Gender (Male): 1 = Male, 0 = Female.
• Age 37-49: 1 = 37-49 years old.
• Age 50+: 1 = 50+ years old.
• Race (White): 1 = White.

# Load required packages
library(tidyverse)
library(ggplot2)
library(modelsummary)
library(vtable)
library(haven)
library(estimatr)
library(ivreg)
library(sandwich)
library(skimr)

# Set global options
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, message = FALSE, warning = FALSE)

# Load data
my_data <- read_dta("https://reifjulian.github.io/illinois-wellness-data/data/stata/claims.dta")

# Rename columns using Base R based on data dictionary
names(my_data)[names(my_data) == "treat"] <- "Random Assignment"
names(my_data)[names(my_data) == "hra_c_yr1"] <- "Participation Status"
names(my_data)[names(my_data) == "covg_0816_0717"] <- "Insurance Coverage Weight"
names(my_data)[names(my_data) == "spend_0715_0716"] <- "Baseline Total Spending"
names(my_data)[names(my_data) == "spend_0816_0717"] <- "Outcome Total Spending"
names(my_data)[names(my_data) == "spendHosp_0715_0716"] <- "Baseline Hospital Spending"
names(my_data)[names(my_data) == "spendOff_0715_0716"] <- "Baseline Office Visit Spending"
names(my_data)[names(my_data) == "spendHosp_0816_0717"] <- "Outcome Hospital Spending"
names(my_data)[names(my_data) == "spendOff_0816_0717"] <- "Outcome Office Visit Spending"
names(my_data)[names(my_data) == "male"] <- "Gender (Male)"
names(my_data)[names(my_data) == "age37_49"] <- "Age 37-49"
names(my_data)[names(my_data) == "age50"] <- "Age 50+"
names(my_data)[names(my_data) == "white"] <- "Race (White)"

# Create grouping variables for analysis
my_data$`Enrollment Status` <- ifelse(my_data$`Participation Status` == 1, "Enrolled", "Not Enrolled")
my_data$`Treatment Group` <- ifelse(my_data$`Random Assignment` == 1, "Treatment", "Control")

# Inspect data
skim(my_data)

── Data Summary ────────────────────────
                           Values 
Name                       my_data
Number of rows             4834   
Number of columns          60     
_______________________           
Column type frequency:            
  character                2      
  numeric                  58     
________________________          
Group variables            None   

Summary Statistics Among Employers With Workplace Wellness Programs

The initial analysis focuses on employees offered the program (Random Assignment == 1), comparing those who chose to enroll versus those who did not.

# Filter for employees offered the program
only_treat <- subset(my_data, `Random Assignment` == 1)

# Summary of employee characteristics by enrollment
sumtable(
  only_treat,
  vars = c(
    "Gender (Male)", "Age 37-49", "Age 50+", "Race (White)",
    "Baseline Hospital Spending", "Baseline Office Visit Spending", "Baseline Total Spending"
  ),
  group = "Enrollment Status",
  title = "Summary of Employee Characteristics by Wellness Program Enrollment",
  out = "kable"
) %>%
  kableExtra::kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))

Summary of Employee Characteristics by Wellness Program Enrollment

Variable	N	Mean	SD	N	Mean	SD
Enrollment Status	Enrolled			Not Enrolled
Gender (Male)	1848	0.4	0.49	1452	0.46	0.5
Age 37-49	1848	0.34	0.47	1452	0.33	0.47
Age 50+	1848	0.31	0.46	1452	0.34	0.47
Race (White)	1848	0.84	0.37	1452	0.84	0.37
Baseline Hospital Spending	1318	222	433	870	316	1078
Baseline Office Visit Spending	1318	59	153	870	56	168
Baseline Total Spending	1318	424	779	870	527	1331

# Summary of health outcomes by enrollment
sumtable(
  only_treat,
  vars = c("Outcome Hospital Spending", "Outcome Office Visit Spending", "Outcome Total Spending"),
  group = "Enrollment Status",
  title = "Summary of Health Outcomes by Wellness Program Enrollment",
  out = "kable"
) %>%
  kableExtra::kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))

Summary of Health Outcomes by Wellness Program Enrollment

Variable	N	Mean	SD	N	Mean	SD
Enrollment Status	Enrolled			Not Enrolled
Outcome Hospital Spending	1338	268	687	870	371	1446
Outcome Office Visit Spending	1338	73	175	870	58	159
Outcome Total Spending	1338	517	1064	870	635	1715

Now, let’s visualize the spending differences.

# Plot distribution of healthcare spending by enrollment
ggplot(
  data = only_treat,
  aes(
    x = log(`Outcome Total Spending` + 1),
    fill = as.factor(`Enrollment Status`),
    color = as.factor(`Enrollment Status`)
  )
) +
  geom_density(alpha = 0.4, na.rm = TRUE) +
  labs(
    title = "Distribution of Total Healthcare Spending by Enrollment",
    x = "Total Spending (2016–2017)",
    y = "Density",
    fill = "Enrollment Status",
    color = "Enrollment Status"
  )

# Box plot of log outcome spending to visualize outliers
ggplot(
  data = only_treat,
  aes(x = `Enrollment Status`, y = log(`Outcome Total Spending` + 1), fill = `Enrollment Status`)
) +
  geom_boxplot() +
  labs(
    title = "Box Plot of Healthcare Spending by Enrollment",
    x = "Enrollment Status",
    y = "Total Spending (2016–2017)",
    fill = "Enrollment Status"
  )

Interpretation: The “Selection Bias” Trap

The table and plots indicate that employees who enrolled in the program spent significantly less money before the program even started (Baseline Spending).

This illustrates Selection Bias. Participants in wellness programs are typically those who are already healthy and health-conscious. Consequently, observed lower costs may reflect pre-existing health status rather than program effectiveness. Comparing enrollees to non-enrollees without adjustment leads to biased conclusions.

Regression Analysis (The “Biased View”)

Standard statistical analysis (OLS Regression) provides estimates if selection bias is ignored.

# Baseline model
baseline <- lm(`Outcome Total Spending` ~ `Participation Status`, data = only_treat, weights = `Insurance Coverage Weight`)

# Add additional models with controls
with_controls <- lm(`Outcome Total Spending` ~ `Participation Status` + `Gender (Male)` + `Age 37-49` + `Age 50+` + `Race (White)`, data = only_treat, weights = `Insurance Coverage Weight`)

# Present in a regression table
models <- list(baseline, with_controls)
modelsummary(models,
  coef_omit = "(Intercept)",
  statistic = c("std.error", "statistic", "conf.int"), fmt = 4,
  gof_map = c("nobs", "r.squared"),
  vcov = "HC1",
  title = "OLS Estimates of the Effect of Wellness Program on Total Spending"
)

OLS Estimates of the Effect of Wellness Program on Total Spending

	(1)	(2)
Participation Status	-137.3006	-152.4476
	(68.5928)	(66.3829)
	(-2.0017)	(-2.2965)
	[-271.8138, -2.7873]	[-282.6273, -22.2680]
Gender (Male)		-276.5082
		(60.4891)
		(-4.5712)
		[-395.1299, -157.8865]
Age 37-49		69.0354
		(55.3430)
		(1.2474)
		[-39.4946, 177.5653]
Age 50+		390.6435
		(84.3748)
		(4.6299)
		[225.1810, 556.1059]
Race (White)		171.1125
		(54.8499)
		(3.1196)
		[63.5494, 278.6755]
Num.Obs.	2208	2208
R2	0.002	0.032

Interpretation: The “Healthy User” Effect

The regression results show a negative coefficient for the program effect, suggesting that enrolling reduces spending by about $137 to $152.

However, this result is misleading as it captures the “Healthy User” effect. The analysis credits the program for the participants’ underlying good health. Since “healthiness” is not perfectly measured by age and gender alone, the regression cannot fully separate the program’s effect from the participants’ pre-existing health. Relying on these results would lead to incorrect business decisions.

Comparison of Treatment and Control Groups (The “Reality Check”)

The Randomized Control Trial (RCT) design addresses this issue. In an RCT, individuals are randomly assigned to be offered the program or not. This randomization creates two statistically identical groups.

# Create a table summarizing the characteristics for treated and untreated enrollees
sumtable(my_data,
  vars = c("Gender (Male)", "Age 37-49", "Age 50+", "Race (White)", "Baseline Hospital Spending", "Baseline Office Visit Spending", "Baseline Total Spending"),
  group = "Treatment Group", title = "Comparison of Baseline Characteristics by Treatment Group", out = "kable"
) %>%
  kableExtra::kable_styling(full_width = TRUE, bootstrap_options = c("striped", "hover", "condensed"))

Comparison of Baseline Characteristics by Treatment Group

Variable	N	Mean	SD	N	Mean	SD
Treatment Group	Control			Treatment
Gender (Male)	1534	0.43	0.49	3300	0.43	0.49
Age 37-49	1534	0.34	0.47	3300	0.33	0.47
Age 50+	1534	0.32	0.47	3300	0.33	0.47
Race (White)	1534	0.84	0.37	3300	0.84	0.37
Baseline Hospital Spending	1035	283	683	2188	259	759
Baseline Office Visit Spending	1035	67	400	2188	58	159
Baseline Total Spending	1035	506	1099	2188	465	1035

# Test whether characteristics are different
t.test(`Gender (Male)` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Gender (Male) by Random Assignment
t = -0.12372, df = 2990.4, p-value = 0.9015
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.03186374  0.02808120
sample estimates:
mean in group 0 mean in group 1 
      0.4256845       0.4275758 

t.test(`Age 37-49` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Age 37-49 by Random Assignment
t = 0.53784, df = 2973.3, p-value = 0.5907
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.02080130  0.03652648
sample estimates:
mean in group 0 mean in group 1 
      0.3402868       0.3324242 

t.test(`Age 50+` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Age 50+ by Random Assignment
t = -0.23005, df = 2996, p-value = 0.8181
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.03170200  0.02504403
sample estimates:
mean in group 0 mean in group 1 
      0.3233377       0.3266667 

t.test(`Race (White)` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Race (White) by Random Assignment
t = 0.45642, df = 3024.1, p-value = 0.6481
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -0.01707684  0.02743913
sample estimates:
mean in group 0 mean in group 1 
      0.8409387       0.8357576 

t.test(`Baseline Hospital Spending` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Baseline Hospital Spending by Random Assignment
t = 0.89886, df = 2235, p-value = 0.3688
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -28.39434  76.45223
sample estimates:
mean in group 0 mean in group 1 
       283.3604        259.3314 

t.test(`Baseline Office Visit Spending` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Baseline Office Visit Spending by Random Assignment
t = 0.67712, df = 1191.9, p-value = 0.4985
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -16.57693  34.04946
sample estimates:
mean in group 0 mean in group 1 
       66.71167        57.97540 

t.test(`Baseline Total Spending` ~ `Random Assignment`, data = my_data)

    Welch Two Sample t-test

data:  Baseline Total Spending by Random Assignment
t = 1.0014, df = 1923.8, p-value = 0.3168
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -39.07482 120.60487
sample estimates:
mean in group 0 mean in group 1 
       505.5786        464.8136 

Interpretation: Why “No Difference” is Good

T-tests were conducted to compare the “Treatment” group (offered the program) and the “Control” group (not offered). The results show no statistically significant difference between them in terms of age, gender, or past spending.

This confirms that the randomization was successful. The two groups are effectively identical. Any difference in future spending can be confidently attributed to the program, not to pre-existing differences.

IV Estimates (The “Causal Truth”)

The Instrumental Variables (IV) technique was utilized to isolate the true effect. The random assignment (Random Assignment) serves as an instrument for enrollment (Participation Status), filtering out selection bias.

# Replace workplace wellness variable with 0 for observations in the control group
my_data$`Participation Status`[is.na(my_data$`Participation Status`)] <- 0

# Estimate the first stage
first_stage <- lm(`Participation Status` ~ `Random Assignment`, data = my_data, weights = `Insurance Coverage Weight`)

# Estimate the reduced form
rf <- lm(`Outcome Total Spending` ~ `Random Assignment`, data = my_data, weights = `Insurance Coverage Weight`)

# Estimate 2SLS
iv <- ivreg(`Outcome Total Spending` ~ `Participation Status` | `Random Assignment`, data = my_data, weights = `Insurance Coverage Weight`)

# Report results in a regression table
iv_models <- list(first_stage, rf, iv)
modelsummary(iv_models,
  coef_omit = "(Intercept)",
  vcov = "robust",
  statistic = c("std.error", "statistic", "conf.int"), fmt = 4,
  gof_map = c("nobs", "r.squared"),
  title = "IV Estimates of the Effect of Wellness Program on Total Spending"
)

IV Estimates of the Effect of Wellness Program on Total Spending

	(1)	(2)	(3)
Random Assignment	0.6135	10.8341
	(0.0071)	(48.5229)
	(86.8366)	(0.2233)
	[0.5997, 0.6274]	[-84.3046, 105.9728]
Participation Status			17.6582
			(79.0955)
			(0.2233)
			[-137.4242, 172.7405]
Num.Obs.	3239	3239	3239
R2	0.334	0.000	-0.000

Interpretation: The Final Verdict

• OLS (Biased): Suggested the program saves ~$150.
• IV (True Causal Effect): The coefficient is not statistically significant (effectively zero).

The difference is stark. Adjusting for selection bias reveals that the savings disappear. The program does not cause employees to spend less. The OLS result was an artifact of healthy people self-selecting into the program.

Conclusion

Recommendation to the Employer:

Do not count on this wellness program to cut medical costs.

Rigorous analysis shows that the program has zero impact on reducing healthcare spending. The apparent savings in simple reports are misleading—reflecting only that healthy employees are more likely to join wellness programs.

If the program is maintained for morale, recruiting, or employee happiness, that is a valid choice. However, if the primary goal is financial ROI from medical cost reduction, this program is not the solution. Alternative strategies or a redesign targeting high-risk employees should be considered.