1  Simulation example

#Load packages
library(ggdag)
library(dagitty)
library(rstanarm)
library(tidybayes)
library(ggdist)
library(tidyverse)
library(gt)
library(rstan)

1.1 Data-generating process

The data-generating process is described via the directed acyclic graph (DAG) in Fig. 1.1. In this DAG, indoor air temperature \((\text{T})\) influences productivity \((\text{P})\) directly, while HVAC system \((\text{H})\) influences productivity \((\text{P})\) both directly (e.g., with air distribution and noise) and indirectly, passing through indoor air temperature \((\text{T})\).

dag_coords.ex1 <-
  data.frame(name = c('T', 'H', 'P'),
             x = c(1, 3.5, 6),
             y = c(1, 3, 1))

DAG.ex1 <-
  dagify(T ~ H,
         P ~ T + H,
         coords = dag_coords.ex1)

node_labels <- c(
  T = 'bold(T)', 
  H = 'bold(H)', 
  P = 'bold(P)'
)


ggplot(data = DAG.ex1, aes(x = x, y = y, xend = xend, yend = yend)) +
  geom_dag_text(aes(label = node_labels[name]),
                parse = TRUE,
                colour = 'black',
                size = 10,
                family = 'mono') +
  geom_dag_edges(arrow_directed = grid::arrow(length = grid::unit(10, 'pt'), type = 'open'),
                 edge_colour = 'black',
                 family = 'mono',
                 fontface = 'bold') +
  annotate('text', x = 1, y = 0.7, label = 'indoor air temperature', 
           size = 4, hjust = 0.4, colour = 'grey50') +
  annotate('text', x = 3.5, y = 3.3, label = 'HVAC system', 
           size = 4, hjust = 0.5, colour = 'grey50') +
  annotate('text', x = 6, y = 0.7, label = 'productivity', 
           size = 4, hjust = 0.5, colour = 'grey50') +
  coord_cartesian(xlim = c(0.5, 6.5), ylim = c(0.8, 3.2))  +
  theme_dag()

Fig. 1.1. Graphical representation via DAG of the data-generating process.

The DAG in Fig. 1.1 can be written as:

• \(T \sim f_{T}(H)\), read as ‘indoor air temperature \((\text{T})\) is some function of HVAC system \((\text{H})\)’.
• \(P \sim f_{P}(T, H)\), read as ‘productivity \((\text{P})\) is some function of indoor air temperature \((\text{T})\) and HVAC system \((\text{H})\)’.

1.2 Synthetic data set

To generate synthetic data, we defined the custom function data.sim_ex1(). This function takes as inputs the sample size n and generates synthetic data according to the DAG in Fig. 1.1.

data.sim_ex1 <- function(n) {
  b_hvac.prod = c(0, 5) #direct causal effect of H on P
  b_temp.prod = -1 #direct causal effect of T on P
  #Simulate HVAC
  H <- factor(sample(c('classic', 'alternative'), size = n, replace = TRUE))  
  #Simulate indoor air temperature
  T <- rnorm(n = n, mean = ifelse(H == 'classic', 21, 22), sd = 0.7)
  #Simulate productivity
  P <- rnorm(n = n, mean = ifelse(H == 'classic', 75 + b_hvac.prod[1] + b_temp.prod * T, 75 + b_hvac.prod[2] + b_temp.prod * T), sd = 2)  
  #Return tibble with simulated values
  return(tibble(H, T, P))
  }

From this data generation mechanism, we simulated the target population, which consists of one million observations.

set.seed(2025) #set random number for reproducibility
#Simulate the population
ex1_population <- data.sim_ex1(n = 1e6)
#Set HVAC reference category to 'classic'
ex1_population$H <- relevel(ex1_population$H, ref = 'classic')
#View the data frame
ex1_population

# A tibble: 1,000,000 × 3
   H               T     P
   <fct>       <dbl> <dbl>
 1 classic      20.6  51.1
 2 alternative  21.6  56.4
 3 alternative  22.3  58.2
 4 alternative  21.1  56.8
 5 classic      20.6  55.7
 6 classic      20.0  53.0
 7 alternative  21.8  56.4
 8 classic      21.7  53.9
 9 alternative  21.5  55.3
10 classic      20.5  57.4
# ℹ 999,990 more rows

From this population, we obtained one data set of five thousand observations using simple random sampling.

n_sims <- 1e3 #number of data sets to simulate

#Set random number for reproducibility
set.seed(2025)  
#Generate a vector of random numbers for reproducibility
ex1_random.seed <- sample(1:1e5, size = n_sims, replace = FALSE)

set.seed(ex1_random.seed[1])
#Sample one data set of 5,000 observations
ex1_sample.random <- 
  ex1_population %>% 
  slice_sample(n = 5e3) #take a simple random sample of size 5,000 

#View the data frame
ex1_sample.random

# A tibble: 5,000 × 3
   H               T     P
   <fct>       <dbl> <dbl>
 1 alternative  22.6  54.3
 2 classic      21.3  53.3
 3 classic      21.5  53.7
 4 alternative  22.5  58.2
 5 classic      22.5  55.6
 6 classic      21.3  52.1
 7 alternative  21.7  60.3
 8 classic      22.9  55.3
 9 alternative  22.7  58.6
10 classic      19.9  55.8
# ℹ 4,990 more rows

1.3 Data analysis

In this example, the target of our analysis is the total average causal effect, ACE (also known as total average treatment effect, ATE) of indoor air temperature \((\text{T})\) on productivity \((\text{P})\), which stands for the expected increase of \(\text{P}\) in response to a unit increase in \(\text{T}\) due to an intervention. The causal effect of interest is visualized in Fig. 1.2.

ggplot(data = DAG.ex1, aes(x = x, y = y, xend = xend, yend = yend)) +
  #visualize causal effect path
  geom_segment(x = 1, xend = 6, y = 1, yend = 1,
               linewidth = 14, lineend = 'round', colour = '#009E73', alpha = 0.05) +
  geom_dag_text(aes(label = node_labels[name]),
                parse = TRUE,
                colour = 'black',
                size = 10,
                family = 'mono') +
  geom_dag_edges(arrow_directed = grid::arrow(length = grid::unit(10, 'pt'), type = 'open'),
                 edge_colour = 'black',
                 family = 'mono',
                 fontface = 'bold') +
  annotate('text', x = 1, y = 0.7, label = 'indoor air temperature', 
           size = 4, hjust = 0.4, colour = 'grey50') +
  annotate('text', x = 3.5, y = 3.3, label = 'HVAC system', 
           size = 4, hjust = 0.5, colour = 'grey50') +
  annotate('text', x = 6, y = 0.7, label = 'productivity', 
           size = 4, hjust = 0.5, colour = 'grey50') +
  #causal effect number
  annotate('text', x = 3.5, y = 1.2, label = '-1', 
           size = 4.5, hjust = 0.5, colour = 'black', parse = TRUE) +

  coord_cartesian(xlim = c(0.5, 6.5), ylim = c(0.8, 3.2))  +
  theme_dag()

Fig. 1.2. Graphical representation via DAG of the data-generating process. The green line indicates the causal question of interest, and the number on the path indicates the total average causal effect.

1.3.1 Identification

The first step to answer the causal question of interest is identification. Identification answers a ‘theoretical’ question by determining whether a causal effect can, in principle, be estimated from observed data. The backdoor criterion and its generalization, the adjustment criterion, allow us to understand whether our causal effect of interest can be identified and, if so, which variables we should (or should not) statistically adjust for (i.e., the adjustment set) to estimate the causal effect from the data.

Given its simplicity, we will first apply the backdoor criterion to identify valid adjustment sets to estimate the causal effect of interest. If the backdoor criterion is not applicable, we will apply its generalization, the adjustment criterion.

Backdoor criterion

Applying the backdoor criterion revealed the existence of a backdoor path (i.e., a non-causal path) from indoor air temperature \((\text{T})\) to productivity \((\text{P})\). Specifically, the backdoor path is \(\text{T} \leftarrow \text{H} \rightarrow \text{P}\). Since HVAC system \((\text{H})\) is a common cause of \(\text{T}\) and \(\text{P}\), association can flow from \(\text{T}\) to \(\text{P}\) through \(\text{H}\). As a result, there is confounding. To close this backdoor path we need to adjust for \(\text{H}\).

Given the DAG in Fig. 1.2, we can use the adjustmentSets() function to identify the adjustment set algorithmically. It is essential to note that this function applies the adjustment criterion and not the backdoor criterion. As such, the adjustmentSets() function can find adjustment sets even when the backdoor criterion is not applicable.

adjustmentSets(DAG.ex1,
               exposure = 'T', #indoor air temperature
               outcome = 'P', #productivity
               type = 'all', 
               effect = 'total', 
               max.results = Inf)

{ H }

As expected, the resulting adjustment set includes HVAC system \((\text{H})\). Therefore, to get the correct estimate of the total average causal effect, we need to adjust for \(\text{H}\); failing to do so will lead to bias.

1.3.2 Estimation

Following the identification step is the estimation step. This step addresses a statistical question by determining how the causal effect identified in the previous step can be estimated. To perform this step, we used a parametric (model-based) estimator, specifically, linear regression. This was possible because we designed the illustrative examples to be simple and with a linear relationship between the variables. This way, we limited the complexity of the examples themselves and shifted the focus to the application of the backdoor criterion to define ‘correct’ adjustment sets.

For transparency and understanding, all (implicit) assumptions used for this illustrative example are (explicitly) provided in Table 1.1.

Table 1.1. Summary description of the simulation example

Research question	Total average causal effect (ACE) of indoor air temperature (T) on productivity (P).
Assumptions	Random sample (simple random sampling): everyone in the population has an equal chance of being selected into the sample.
	Limited random variability: large sample size.
	Independence of observations: each observation represents independent bits of information.
	No confounding: the DAG includes all shared causes among the variables.
	No model error: perfect functional form specification.
	No measurement error: all variables are measured perfectly.
Variables	Indoor air temperature (T): continuous variable [unit: °C]
	HVAC (H): categorical variable ['classic'; 'alternative']
	Productivity (P): continuous variable [unit: -]

To carry out the estimation step, we utilized linear regression within both the frequentist and Bayesian frameworks. Specifically, we will run two regression models:

• Model.1 will include only indoor air temperature \((\text{T})\) as predictor;
• Model.2 will include both indoor air temperature \((\text{T})\) and HVAC system \((\text{H})\) as predictors.

The results of the fitted statistical models (i.e., Model.1 and Model.2) are presented here.

Frequentist framework

#Fit the linear regression model with T and P (Model 1)
ex1_Model.1 <-
  lm(formula = P ~ T,
     data = ex1_sample.random)
#View of the model summary
summary(ex1_Model.1)

Call:
lm(formula = P ~ T, data = ex1_sample.random)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.7246 -2.0653  0.0144  1.9700 10.0478 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 40.56735    1.00976   40.17   <2e-16 ***
T            0.71541    0.04693   15.24   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.869 on 4998 degrees of freedom
Multiple R-squared:  0.04443,   Adjusted R-squared:  0.04423 
F-statistic: 232.4 on 1 and 4998 DF,  p-value: < 2.2e-16

#Fit the linear regression model with T , H and P (Model 2)
ex1_Model.2 <-
  lm(formula = P ~ T + H,
     data = ex1_sample.random)
#View of the model summary
summary(ex1_Model.2)

Call:
lm(formula = P ~ T + H, data = ex1_sample.random)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.6805 -1.3390 -0.0401  1.3199  8.8070 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  75.15431    0.84857   88.56   <2e-16 ***
T            -1.00927    0.04036  -25.00   <2e-16 ***
Halternative  5.08393    0.06979   72.84   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.998 on 4997 degrees of freedom
Multiple R-squared:  0.5365,    Adjusted R-squared:  0.5364 
F-statistic:  2892 on 2 and 4997 DF,  p-value: < 2.2e-16

The estimated coefficients for the two models are then plotted in Fig. 1.3.

b_temp.prod = -1 

data.frame(model = c('Model.1', 'Model.2'),
           estimate = c(coef(ex1_Model.1)['T'], coef(ex1_Model.2)['T']),
           lower.95.CI = c(confint(ex1_Model.1, level = 0.95, type = 'Wald')['T', 1], confint(ex1_Model.2, level = 0.95, type = 'Wald')['T', 1]),
           upper.95.CI = c(confint(ex1_Model.1, level = 0.95, type = 'Wald')['T', 2], confint(ex1_Model.2, level = 0.95, type = 'Wald')['T', 2])) %>%
  
ggplot(aes(x = estimate, y = model, xmin = lower.95.CI, xmax = upper.95.CI)) + 
  geom_vline(xintercept = b_temp.prod, alpha = 0.8, linetype = 'dashed', colour = 'blue') + 
  geom_linerange() +
  geom_point(shape = 21, size = 2, fill = 'white', stroke = 1) +
  scale_x_continuous('Estimate', 
                     breaks = seq(from = -5, to = 5, by = 0.5),
                     limits = c(-1.25, 1.25)) +
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA),
        axis.title.y = element_blank())

Fig. 1.3. Estimates of the temperature coefficient for the two models. The white dots represent the point estimate, and the black lines represent the 95% confidence intervals. The blue dashed line represents the parameter used to generate the data.

Fig. 1.3 shows the estimates (point estimate and 95% confidence interval) of the coefficient for temperature for Model.1 and Model.2.

For Model.1 we found a positive coefficient between indoor air temperature \((\text{T})\) and productivity \((\text{P})\) equal to 0.715 with 95% confidence interval (CI) [0.623, 0.807]. Since the 95% CI does not include zero, the regression coefficient is statistically significantly different from zero at the 0.05 level (p-value = 2.56e-51). Additionally, since the 95% CI does not include the data-generating parameter for indoor air temperature (i.e., b_temp.prod = -1), we can deduce that the estimated coefficient for indoor air temperature is statistically significantly different from -1 at the 0.05 level (although this will be the case for all numbers within the 95% confidence interval). We can test this more formally by using the linearHypothesis() function. Specifically,

#Test for b_temp.prod = -1
car::linearHypothesis(ex1_Model.1, 'T = -1') # car:: access the function from the car package without loading it in the environment 

Linear hypothesis test:
T = - 1

Model 1: restricted model
Model 2: P ~ T

  Res.Df   RSS Df Sum of Sq    F    Pr(>F)    
1   4999 52123                                
2   4998 41130  1     10994 1336 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The resulting p-value is 1.97e-259, which indicates that we can reject the null hypothesis (i.e., T = -1) at the 0.05 level. This result suggests that the regression coefficient for indoor air temperature (i.e., 0.715) is statistically significantly different from -1. Since the estimated causal effect from Model.1 is biased, we would erroneously conclude that an increase of 1°C in indoor air temperature causes an increase in productivity by 0.715 units (95% CI [0.623, 0.807]).

For Model.2 we found a negative coefficient between indoor air temperature \((\text{T})\) and productivity \((\text{P})\) equal to -1.009 with 95% CI [-1.088, -0.930]. Since the 95% CI excludes zero, the regression coefficient is statistically significantly different from zero at the 0.05 level (p-value = 3.99e-130). Additionally, since the 95% CI includes the data-generating parameter for indoor air temperature (i.e., b_temp.prod = -1), we can deduce that the estimated coefficient for indoor air temperature is not statistically significantly different from -1 at the 0.05 level (although this will be the case for all numbers within the 95% confidence interval). We can test this more formally by using the linearHypothesis() function. Specifically,

#Test for T = -1
car::linearHypothesis(ex1_Model.2, 'T = -1') # car:: access the function from the car package without loading it in the environment 

Linear hypothesis test:
T = - 1

Model 1: restricted model
Model 2: P ~ T + H

  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1   4998 19948                           
2   4997 19948  1   0.21046 0.0527 0.8184

The resulting p-value is 0.818, which indicates that we fail to reject the null hypothesis (i.e., T = -1) at the 0.05 level. This result suggests that the regression coefficient for indoor air temperature (i.e., -1.009) is not statistically significantly different from -1. Since the estimated causal effect from Model.2 is unbiased, we would correctly conclude that an increase of 1°C in indoor air temperature causes a decrease in productivity by -1.009 units (95% CI [-1.088, -0.930]).

Importantly, for both Model.1 and Model.2, the 95% confidence interval means that if we were to repeat the sampling process and calculate the interval many times, 95% of those calculated intervals would contain the true population parameter. To highlight this, we can repeat the analysis by fitting the two models to one thousand data sets randomly selected from our population.

The for-loop shown in the code below performs the following operations. First, sample (using simple random sampling) a data set of 5,000 observations from the target population. Subsequently, perform linear regression using Model.1 and Model.2 and store the estimated coefficients for temperature, its standard error and 95% confidence interval in the data frame coefs_ex1. This operation is repeated a thousand times, resulting in the data frame coefs_ex1 containing the estimates (point estimate, standard error and confidence interval) of a thousand random samples of size 5,000 using both Model.1 and Model.2.

n_model <- c('mod.1', 'mod.2')     

n_row <- n_sims*length(n_model)

#Create an empty data frame
empty.df <- data.frame(matrix(NA, nrow = n_row, ncol = 7))
#Rename the data frame columns
colnames(empty.df) <- c('sim.id', 'estimate', 'se', 'CI_2.5', 'CI_97.5',
                        'model', 'coverage')

#Sample a thousand data sets of 5,000 observations and perform linear regression 
coefs_ex1 <- empty.df #assign the empty data frame
k = 1
for (i in 1:n_sims){
  set.seed(ex1_random.seed[i]) #set unique seed for each simulation 
#Sample data set from population   
  sample.random <- 
    ex1_population %>% 
    slice_sample(n = 5e3) #take a simple random sample of size 5,000
#Fit models
  for (j in 1:length(n_model)){
    if (n_model[j] == 'mod.1'){
      fit <- lm(formula = P ~ T,
                data = sample.random)
    } else {
      fit <- lm(formula = P ~ T + H,
                data = sample.random)
    }  
#Compile matrix  
  coefs_ex1[k, 1] <- i #simulation ID
  coefs_ex1[k, 2] <- coef(fit)['T'] #point estimate
  coefs_ex1[k, 3] <- summary(fit)$coef['T','Std. Error'] #standard error
  coefs_ex1[k, 4:5] <- confint(fit, level = 0.95, type = 'Wald')['T', ] #confidence interval (Wald)
  coefs_ex1[k, 6] <- n_model[j] #sample size
  k = k + 1
  }
}
coefs_ex1 <- as_tibble(coefs_ex1)
#View the data frame
coefs_ex1

# A tibble: 2,000 × 7
   sim.id estimate     se CI_2.5 CI_97.5 model coverage
    <int>    <dbl>  <dbl>  <dbl>   <dbl> <chr> <lgl>   
 1      1    0.715 0.0469  0.623   0.807 mod.1 NA      
 2      1   -1.01  0.0404 -1.09   -0.930 mod.2 NA      
 3      2    0.706 0.0477  0.613   0.800 mod.1 NA      
 4      2   -0.959 0.0415 -1.04   -0.878 mod.2 NA      
 5      3    0.673 0.0472  0.581   0.766 mod.1 NA      
 6      3   -1.02  0.0402 -1.10   -0.946 mod.2 NA      
 7      4    0.668 0.0468  0.576   0.759 mod.1 NA      
 8      4   -1.02  0.0405 -1.10   -0.937 mod.2 NA      
 9      5    0.717 0.0475  0.624   0.811 mod.1 NA      
10      5   -1.03  0.0418 -1.11   -0.949 mod.2 NA      
# ℹ 1,990 more rows

The coverage is defined by setting its value to 1 if the confidence interval overlaps the data-generating parameter for temperature (i.e., b_temp.prod = -1) and 0 otherwise.

#Calculate coverage
coefs_ex1 <-
  coefs_ex1 %>%
  mutate(coverage = case_when(CI_2.5 > b_temp.prod | CI_97.5 < b_temp.prod ~ 0,
                              CI_2.5 <= b_temp.prod & CI_97.5 >= b_temp.prod ~ 1,
                              .default = NA))
#View the data frame
coefs_ex1

# A tibble: 2,000 × 7
   sim.id estimate     se CI_2.5 CI_97.5 model coverage
    <int>    <dbl>  <dbl>  <dbl>   <dbl> <chr>    <dbl>
 1      1    0.715 0.0469  0.623   0.807 mod.1        0
 2      1   -1.01  0.0404 -1.09   -0.930 mod.2        1
 3      2    0.706 0.0477  0.613   0.800 mod.1        0
 4      2   -0.959 0.0415 -1.04   -0.878 mod.2        1
 5      3    0.673 0.0472  0.581   0.766 mod.1        0
 6      3   -1.02  0.0402 -1.10   -0.946 mod.2        1
 7      4    0.668 0.0468  0.576   0.759 mod.1        0
 8      4   -1.02  0.0405 -1.10   -0.937 mod.2        1
 9      5    0.717 0.0475  0.624   0.811 mod.1        0
10      5   -1.03  0.0418 -1.11   -0.949 mod.2        1
# ℹ 1,990 more rows

The results are then plotted in Fig. 1.4.

model_names <- c('mod.1' = 'Model.1', 
                 'mod.2' = 'Model.2')

ggplot(data = subset(coefs_ex1, sim.id <= 2e2), aes(x = sim.id, y = estimate, ymin = CI_2.5, ymax = CI_97.5, colour = as.factor(coverage))) + 
  geom_hline(yintercept = b_temp.prod, alpha = 0.8, linetype = 'dashed', colour = 'blue') + 
  geom_linerange() +
  geom_point(shape = 21, size = 1, fill = 'white', stroke = 0.5) +
  scale_colour_manual('', values = c('black', '#FF4040'),
                      breaks = c('1','0')) +
  scale_y_continuous('Estimate') +
  scale_x_continuous('Simulation ID') +
  facet_wrap(model~., 
             labeller = labeller(model = model_names),
             nrow = 2,
             scales = 'fixed') +
  theme(legend.position = 'none', 
        panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA),
        axis.title.y = element_blank())

Fig. 1.4. Estimates of the temperature coefficient (only the first hundreds of the thousand simulations are shown). The white dots represent the point estimate, and the black lines represent the 95% confidence intervals. In red are highlighted the confidence intervals that do not overlap the parameter used to generate the data (dashed blue line).

Fig. 1.4 shows the first 200 estimates (point estimate and confidence interval) of the coefficient for indoor air temperature for Model.1 and Model.2. If all the thousand simulations are considered, the frequency of the coverage of the calculated confidence intervals (i.e., how many times the confidence intervals overlap the data-generating parameter) is 0.0% and 94.9% for Model.1 and Model.2, respectively. Since the estimates of the causal effect for Model.1 are biased, the calculated confidence intervals for this model do not have the expected coverage. In fact, confidence intervals only quantify the uncertainty due to random error (i.e., sample variability), not systematic error (i.e., bias). Instead, the estimates from Model.2 are unbiased and the calculated 95% confidence intervals have the expected coverage (i.e., overlap the data-generating parameter 95% of the times).

We can also visualize all the 1,000 estimates of the coefficient for indoor air temperature for Model.1 and Model.2 (the estimates were shown as white dots in Fig. 1.4).

ggplot(data = coefs_ex1, aes(x = estimate, y = ifelse(after_stat(count) > 0, after_stat(count), NA))) +
  geom_histogram(binwidth = 0.05, fill = 'white', colour = 'grey50') +
  geom_vline(xintercept = -1, colour = 'blue', linetype = 'dashed') +
  scale_x_continuous('Estimate', 
                     breaks = seq(from = -5, to = 5, by = 0.5),
                     limits = c(-1.5, 1.5)) +
  facet_grid(model~., 
             labeller = labeller(model = model_names),
             scales = 'fixed') +
  theme(axis.title.y = element_blank(), 
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank(),
        panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA))

Fig. 1.5. Histogram of the 1,000 estimates of the temperature coefficient for the two models. The blue dashed line represents the parameter used to generate the data.

In Fig. 1.5, the estimated coefficients are visualized through a histogram. Here, we can see that the estimates for Model.1 are not clustered around the data-generating parameter (blue dashed line) having a mean estimate of 0.687 (with 0.043 standard deviation). In contrast, the estimates for Model.2 having a mean estimate of -1.005 (with 0.041 standard deviation) are centered around the data-generating parameter.

As expected, the inclusion of HVAC system \((\text{H})\) in the regression model (i.e., using Model.2) leads to the correct estimate of the total average causal effect, while its exclusion (i.e., using Model.1) leads to bias. Specifically, the estimated effect is reversed. This phenomenon is known as Simpson’s Paradox. Generally, Simpson’s Paradox happens when a trend that appears in different groups of data reverses (or disappears) when the groups are combined.

To better explain this apparent paradox, we can plot the data and the fitted regression lines for only the data with HVAC = classic, the data with HVAC = alternative, and the data combined (i.e., without considering the groups).

pred_ex1_Model.1_separete <-
  rbind(cbind(filter(ex1_sample.random, H == 'classic'), 
              predict(lm(formula = P ~ T, data = filter(ex1_sample.random, H == 'classic')), 
                      interval = 'confidence', level = 0.95), pred = 'classic'), 
        cbind(filter(ex1_sample.random, H == 'alternative'), 
              predict(lm(formula = P ~ T, data = filter(ex1_sample.random, H == 'alternative')),
                      interval = 'confidence', level = 0.95), pred = 'alternative'))

pred_ex1_Model.1 <- 
  cbind(ex1_sample.random, 
      predict(ex1_Model.1, interval = 'confidence', level = 0.95))

#Plot 
ex1_sample.random %>%
  ggplot(aes(x = T, y = P)) +
  geom_point(aes(colour = H), shape = 19, size = 2, alpha = 0.10, stroke = NA) +
  geom_ribbon(data = pred_ex1_Model.1_separete, aes(x = T , y = fit, ymin = lwr, ymax = upr, fill = H), alpha = 0.2, show.legend = NA) +
  geom_line(data = pred_ex1_Model.1_separete, aes(x = T , y = fit, colour = H)) + 
  geom_ribbon(data = pred_ex1_Model.1, aes(x = T , y = fit, ymin = lwr, ymax = upr, fill = 'all.data'), alpha = 0.2, show.legend = NA) +
  geom_line(data = pred_ex1_Model.1, aes(x = T , y = fit, colour = 'all.data')) + 
  scale_colour_manual('',
                      breaks = c('classic', 'alternative', 'all.data'),
                      labels = c('classic', 'alternative', 'all data'),
                      values = c('#D55E00','#0072B2', 'black')) +
  scale_fill_manual('',
                    breaks = c('classic', 'alternative', 'all.data'),
                    labels = c('classic', 'alternative', 'all data'),
                    values = c('#D55E00','#0072B2', 'black')) +
  scale_x_continuous('Indoor air temperature') + 
  scale_y_continuous('Productivity') +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) + 
  theme(legend.position = 'bottom',
        legend.direction = 'horizontal',
        panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA))

Fig. 1.6. Visualization of the reversed regression line direction in the total sample compared to the subgroups.

In Fig. 1.6 the three colored lines represent the regression lines for only the data with HVAC = classic (in orange), the data with HVAC = alternative (in blue), and the combined data (in black). The regression lines in orange and blue are obtained by fitting the following regression models:

• subset of data with H == 'classic'

#Filter data with H == 'classic'
ex1_sample.random_classic <- filter(ex1_sample.random, H == 'classic')
#Fit the linear regression model for only data with H == 'classic'
ex1_Model.1_classic <- lm(formula = P ~ T, data = ex1_sample.random_classic)
#View of the model summary
ex1_Model.1_classic

Call:
lm(formula = P ~ T, data = ex1_sample.random_classic)

Coefficients:
(Intercept)            T  
    74.1501      -0.9614  

• subset of data with H == 'alternative'

#Filter data with H == 'alternative'
ex1_sample.random_alternative <- filter(ex1_sample.random, H == 'alternative')
#Fit the linear regression model for only data with H == 'alternative'
ex1_Model.1_alternative <- lm(formula = P ~ T, data = ex1_sample.random_alternative)      
#View of the model summary
ex1_Model.1_alternative        

Call:
lm(formula = P ~ T, data = ex1_sample.random_alternative)

Coefficients:
(Intercept)            T  
      81.36        -1.06  

We can visualize the estimated coefficients for these new regression models and compare them with the ones estimated from Model.1 and Model.2.

#Plot
data.frame(model = c('Model.1', 'Model.1_classic', 'Model.1_alternative', 'Model.2'),
           estimate = c(coef(ex1_Model.1)['T'], 
                        coef(ex1_Model.1_classic)['T'], 
                        coef(ex1_Model.1_alternative)['T'], 
                        coef(ex1_Model.2)['T']),
           lower.95.CI = c(confint(ex1_Model.1, level = 0.95, type = 'Wald')['T', 1],
                           confint(ex1_Model.1_classic, level = 0.95, type = 'Wald')['T', 1],
                           confint(ex1_Model.1_alternative, level = 0.95, type = 'Wald')['T', 1],
                           confint(ex1_Model.2, level = 0.95, type = 'Wald')['T', 1]),
           upper.95.CI = c(confint(ex1_Model.1, level = 0.95, type = 'Wald')['T', 2], 
                           confint(ex1_Model.1_classic, level = 0.95, type = 'Wald')['T', 2], 
                           confint(ex1_Model.1_alternative, level = 0.95, type = 'Wald')['T', 2], 
                           confint(ex1_Model.2, level = 0.95, type = 'Wald')['T', 2])) %>%
  
ggplot(aes(x = estimate, y = model, xmin = lower.95.CI, xmax = upper.95.CI)) + 
  geom_vline(xintercept = b_temp.prod, alpha = 0.8, linetype = 'dashed', colour = 'blue') + 
  geom_linerange() +
  geom_point(shape = 21, size = 2, fill = 'white', stroke = 1) +
  scale_x_continuous('Estimate', 
                     breaks = seq(from = -5, to = 5, by = 0.5),
                     limits = c(-1.25, 1.25)) +
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA),
        axis.title.y = element_blank())

Fig. 1.7. Estimates of the temperature coefficient for the four models. The white dots represent the point estimate, and the black lines represent the 95% confidence intervals. The blue dashed line represents the parameter used to generate the data.

Fig. 1.7, we can see that only the estimated coefficient for Model.1 is biased. As explained in the identification step, to identify the causal effect of indoor air temperature \((\text{T})\) on productivity \((\text{P})\), we need to block the backdoor path (i.e., the non-causal path) \(\text{T} \leftarrow \text{H} \rightarrow \text{P}\). Fitting a separate regression model to the subset of data with H == 'classic' and one to H == 'alternative' is equivalent to conditioning for \(\text{H}\), closing the backdoor path and allowing unbiased estimation of the causal effect of indoor air temperature \((\text{T})\) on productivity \((\text{P})\) (although this comes at the price of a larger standard error and consequently wider confidence intervals compared to Model.2, where \(\text{H}\) is added as a predictor).

Bayesian framework

Click to expand

#Fit the linear regression model with T and P (Model 1)
ex1_Model.1 <- stan_glm(formula = P ~ T,
                        family = gaussian(),
                        data = ex1_sample.random,
                        #Prior coefficients
                        prior = normal(location = 0, scale = 2),
                        #Prior intercept
                        prior_intercept = normal(location = 56, scale = 10),
                        #Prior sigma
                        prior_aux = exponential(rate = 0.5),
                        iter = 4000, warmup = 1000, 
                        save_warmup = TRUE,
                        chains = 4, cores = 4,
                        seed = 2025) #for reproducibility   
#View of the model
ex1_Model.1

stan_glm
 family:       gaussian [identity]
 formula:      P ~ T
 observations: 5000
 predictors:   2
------
            Median MAD_SD
(Intercept) 40.6    1.0  
T            0.7    0.0  

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.9    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

#Fit the linear regression model with T , H and P (Model 2)
ex1_Model.2 <- stan_glm(formula = P ~ T + H,
                        family = gaussian(),
                        data = ex1_sample.random,
                        #Prior coefficients
                        prior = normal(location = c(0, 0), scale = c(2, 5)),
                        #Prior intercept
                        prior_intercept = normal(location = 56, scale = 10),
                        #Prior sigma
                        prior_aux = exponential(rate = 0.5),
                        iter = 4000, warmup = 1000, 
                        save_warmup = TRUE,
                        chains = 4, cores = 4,
                        seed = 2025) #for reproducibility
#View of the model
ex1_Model.2

stan_glm
 family:       gaussian [identity]
 formula:      P ~ T + H
 observations: 5000
 predictors:   3
------
             Median MAD_SD
(Intercept)  75.1    0.9  
T            -1.0    0.0  
Halternative  5.1    0.1  

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.0    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

In Bayesian analysis, there are important diagnostics that have to be carried out in order to assess the convergence and efficiency of the Markov Chains. This is done by using the monitor() function which computes summaries of MCMC (Markov Chain Monte Carlo) draws and monitor convergence. Specifically, we will look at Rhat, Bulk_ESS and Tail_ESS metrics.

#Diagnostics for model 1
monitor(ex1_Model.1$stanfit)  

Inference for the input samples (4 chains: each with iter = 4000; warmup = 0):

                    Q5      Q50      Q95     Mean  SD  Rhat Bulk_ESS Tail_ESS
(Intercept)       39.0     40.6     42.2     40.6 1.0     1    11517     8203
T                  0.6      0.7      0.8      0.7 0.0     1    11445     8323
sigma              2.8      2.9      2.9      2.9 0.0     1    12774     8565
mean_PPD          55.9     55.9     56.0     55.9 0.1     1    12512    11016
log-posterior -12372.1 -12369.4 -12368.4 -12369.7 1.2     1     5351     7267

For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
effective sample size for bulk and tail quantities respectively (an ESS > 100 
per chain is considered good), and Rhat is the potential scale reduction 
factor on rank normalized split chains (at convergence, Rhat <= 1.05).

#Diagnostics for model 2
monitor(ex1_Model.2$stanfit) 

Inference for the input samples (4 chains: each with iter = 4000; warmup = 0):

                    Q5      Q50      Q95     Mean  SD  Rhat Bulk_ESS Tail_ESS
(Intercept)       73.7     75.1     76.5     75.1 0.9     1     8648     8845
T                 -1.1     -1.0     -0.9     -1.0 0.0     1     8495     8786
Halternative       5.0      5.1      5.2      5.1 0.1     1     8517     8685
sigma              2.0      2.0      2.0      2.0 0.0     1    12597     8666
mean_PPD          55.9     55.9     56.0     55.9 0.0     1    11453    11001
log-posterior -10565.3 -10562.3 -10561.0 -10562.6 1.4     1     5905     7324

For each parameter, Bulk_ESS and Tail_ESS are crude measures of 
effective sample size for bulk and tail quantities respectively (an ESS > 100 
per chain is considered good), and Rhat is the potential scale reduction 
factor on rank normalized split chains (at convergence, Rhat <= 1.05).

Rhat is a metric used to assess the convergence of Markov Chain Monte Carlo (MCMC) simulations. It helps determine if the MCMC chains have adequately explored the target posterior distribution. Specifically, it compares the between- and within-chain estimates for model parameters: If chains have not mixed well (i.e., the between- and within-chain estimates do not agree), R-hat is larger than 1. A general rule of thumb is to use the sample only if R-hat is less than 1.05; a larger value suggests that the chains have not mixed well, and the results might not be reliable. In our two models, all Rhat are equal to 1 indicating that the chains have mixed well and have adequately explored the target posterior distribution.

Bulk_ESS and Tail_ESS stand for ‘Bulk Effective Sample Size’ and ‘Tail Effective Sample Size,’ respectively. Since MCMC samples are not truly independent (they are correlated), these metrics assess the sampling efficiencies, that is, they help evaluate how efficiently the MCMC sampler is exploring the parameter space.

• Bulk_ESS is a useful measure for sampling efficiency in the bulk (center) of the distribution (e.g., efficiency of mean and median estimates);
• Tail_ESS is a useful measure for sampling efficiency in the tails of the distribution (e.g., efficiency of variance and tail quantile estimates). A general rule of thumb is that both Bulk-ESS and Tail-ESS should be at least 100 (approximately) per Markov Chain in order to be reliable and indicate that estimates of respective posterior quantiles are reliable. In our two models, all Bulk-ESS and Tail-ESS are well above 400 (i.e., 100 multiplied by 4, the number of chains we used) indicating that estimates of posterior quantiles are reliable.

Since we have established that the posteriors are reliable, we can now explore the model estimates.

The estimated coefficients for the two models are then plotted in Fig. 1.8.

#Extract draws from model 1 
post_ex1_Model.1 <-
  ex1_Model.1 %>% 
  spread_draws(T) %>% #extract draws from the fitted model
  mutate(model = 'Model.1') #add a new column to specify that the model

#Extract draws from model 2
post_ex1_Model.2 <-
  ex1_Model.2 %>% 
  spread_draws(T) %>% #extract draws from the fitted model
  mutate(model = 'Model.2') #add a new column to specify that the model

#Combine draws
plot.post <- rbind(post_ex1_Model.1, post_ex1_Model.2)

# Plot
plot.post  %>%
  ggplot(aes(y = model, x = T)) +
  stat_slabinterval(point_interval = 'mean_hdi',
                    .width = c(.95)) +
  geom_vline(xintercept = b_temp.prod, alpha = 0.8, linetype = 'dashed', colour = 'blue') + 
  scale_x_continuous('Estimate', 
                     breaks = seq(from = -5, to = 5, by = 0.5),
                     limits = c(-1.25, 1.25)) +
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA),
        axis.title.y = element_blank())

Fig. 1.8. Posterior distribution of the temperature coefficient for the two models. The black line and dot at the bottom of each distribution represent the highest density interval (HDI) and the mean, respectively.

Fig. 1.8 shows the estimates (mean and 95% HDI) of the coefficient for temperature for Model.1 and Model.2.

For Model.1 we found a positive coefficient between indoor air temperature \((\text{T})\) and productivity \((\text{P})\) (mean = 0.714, 95% HDI [0.624, 0.803]). The estimated causal effect is therefore biased, leading to the wrong conclusion that an increase of 1°C in indoor air temperature causes an increase in productivity by 0.714 units, on average. In contrast, for Model.2 we found a negative coefficient between indoor air temperature \((\text{T})\) and productivity \((\text{P})\) (mean = -1.008, 95% HDI [-1.087, -0.929]). The estimated causal effect is therefore unbiased, leading to the correct conclusion that an increase of 1°C in indoor air temperature causes a decrease in productivity by -1.008 units, on average.

Importantly, for both Model.1 and Model.2, the 95% HDI is the range of parameter values within which the most credible 95% of the posterior distribution falls. Unlike a frequentist confidence interval, the Bayesian 95% HDI has a direct probabilistic meaning: every point inside the HDI has a higher probability density than any point outside the interval. Therefore, given the model, the prior and the data, we can say that there is a 95% probability that the data-generating parameter (i.e., b_temp.prod = -1) lies within the HDI. However, since Model.1 leads to a biased estimate, we will reach the wrong conclusion by stating that there is a 95% probability that the data-generating parameter lies within the [0.624, 0.803] interval. This probability is 0%.

As expected, the inclusion of HVAC system \((\text{H})\) in the regression model (i.e., using Model.2) leads to the correct estimate of the total average causal effect, while its exclusion (i.e., using Model.1) leads to bias. Specifically, the estimated effect is reversed. This phenomenon is known as Simpson’s Paradox. Generally, Simpson’s Paradox happens when a trend that appears in different groups of data reverses (or disappears) when the groups are combined.

To better explain this apparent paradox, we can plot the data and the fitted regression lines for only the data with HVAC = classic, the data with HVAC = alternative, and the data combined (i.e., without considering the groups).

#Expected value of the posterior predictive distribution (epred) for model 1
epred_ex1_Model.1 <- 
  ex1_Model.1 %>%
  epred_draws(newdata = ex1_sample.random, 
              seed = 2025) %>%
  group_by(T) %>%
  mean_hdi(.epred, .width = .95) %>%
  ungroup()

#Expected value of the posterior predictive distribution (epred) for model 2
epred_ex1_Model.2 <- 
  ex1_Model.2 %>%
  epred_draws(newdata = ex1_sample.random, 
              seed = 2025) %>%
  group_by(T, H) %>%
  mean_hdi(.epred, .width = .95) %>%
  ungroup()

#Plot 
ex1_sample.random %>%
  ggplot(aes(x = T, y = P)) +
  geom_point(aes(colour = H), shape = 19, size = 2, alpha = 0.15, stroke = NA) +
  
  geom_ribbon(data = epred_ex1_Model.2, aes(x = T , y = .epred, ymin = .lower, ymax = .upper, fill = H), alpha = 0.2, show.legend = NA) +
  geom_line(data = epred_ex1_Model.2, aes(x = T , y = .epred, colour = H)) + 
  
  geom_ribbon(data = epred_ex1_Model.1, aes(x = T , y = .epred, ymin = .lower, ymax = .upper, fill = 'all.data'), alpha = 0.2, show.legend = NA) +
  geom_line(data = epred_ex1_Model.1, aes(x = T , y = .epred, colour = 'all.data')) + 
  scale_colour_manual('',
                      breaks = c('classic', 'alternative', 'all.data'),
                      labels = c('classic', 'alternative', 'all data'),
                      values = c('#D55E00','#0072B2', 'black')) +
  scale_fill_manual('',
                    breaks = c('classic', 'alternative', 'all.data'),
                    labels = c('classic', 'alternative', 'all data'),
                    values = c('#D55E00','#0072B2', 'black')) +
  scale_x_continuous('Indoor air temperature') + 
  scale_y_continuous('Productivity') +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) + 
  theme(legend.position = 'bottom',
        legend.direction = 'horizontal',
        panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA))

Fig. 1.9. Visualization of the reversed regression line direction in the total sample compared to the subgroups.

In Fig. 1.9 the three colored lines represent the regression lines for only the data with HVAC = classic (in orange), the data with HVAC = alternative (in blue), and the combined data (in black). The regression lines in orange and blue are obtained by fitting the following regression models:

• subset of data with H == 'classic'

#Filter data with H == 'classic'
ex1_sample.random_classic <- filter(ex1_sample.random, H == 'classic')
#Fit the linear regression model for only data with H == 'classic'
ex1_Model.1_classic <- stan_glm(formula = P ~ T,
                                    family = gaussian(),
                                    data = ex1_sample.random_classic,
                                    #Prior coefficients
                                    prior = normal(location = 0, scale = 5),
                                    #Prior intercept
                                    prior_intercept = normal(location = 56, scale = 10),
                                    #Prior sigma
                                    prior_aux = exponential(rate = 0.5),
                                    iter = 4000, warmup = 1000, 
                                    save_warmup = TRUE,
                                    chains = 4, cores = 4,
                                    seed = 2025) #for reproducibility
#View of the model
ex1_Model.1_classic

stan_glm
 family:       gaussian [identity]
 formula:      P ~ T
 observations: 2551
 predictors:   2
------
            Median MAD_SD
(Intercept) 74.1    1.2  
T           -1.0    0.1  

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.0    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

• subset of data with H == 'alternative'

#Filter data with H == 'alternative'
ex1_sample.random_alternative <- filter(ex1_sample.random, H == 'alternative')
#Fit the linear regression model for only data with H == 'alternative'
ex1_Model.1_alternative <- stan_glm(formula = P ~ T,
                                    family = gaussian(),
                                    data = ex1_sample.random_alternative,
                                    #Prior coefficients
                                    prior = normal(location = 0, scale = 5),
                                    #Prior intercept
                                    prior_intercept = normal(location = 56, scale = 10),
                                    #Prior sigma
                                    prior_aux = exponential(rate = 0.5),
                                    iter = 4000, warmup = 1000, 
                                    save_warmup = TRUE,
                                    chains = 4, cores = 4,
                                    seed = 2025) #for reproducibility
#View of the model
ex1_Model.1_alternative

stan_glm
 family:       gaussian [identity]
 formula:      P ~ T
 observations: 2449
 predictors:   2
------
            Median MAD_SD
(Intercept) 81.4    1.3  
T           -1.1    0.1  

Auxiliary parameter(s):
      Median MAD_SD
sigma 2.0    0.0   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

We can visualize the estimated coefficients for these new regression models and compared them with the one estimated from Model.1 and Model.2.

#Extract draws from model 1 
post_ex1_Model.1_classic <-
  ex1_Model.1_classic %>% 
  spread_draws(T) %>% #extract draws from the fitted model
  mutate(model = 'Model.1_classic') #add a new column to specify that the model

#Extract draws from model 2
post_ex1_Model.1_alternative <-
  ex1_Model.1_alternative %>% 
  spread_draws(T) %>% #extract draws from the fitted model
  mutate(model = 'Model.1_alternative') #add a new column to specify that the model

#Combine draws
plot.post1 <- rbind(post_ex1_Model.1, post_ex1_Model.1_classic, post_ex1_Model.1_alternative, post_ex1_Model.2)

#Plot
plot.post1  %>%
  ggplot(aes(y = model, x = T)) +
  stat_slabinterval(point_interval = 'mean_hdi',
                    .width = c(.95)) +
  geom_vline(xintercept = b_temp.prod, alpha = 0.8, linetype = 'dashed', colour = 'blue') + 
  scale_x_continuous('Estimate', 
                     breaks = seq(from = -5, to = 5, by = 0.5),
                     limits = c(-1.5, 1.5)) +
  theme(panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(colour = 'black', fill = NA),
        axis.title.y = element_blank())

Fig. 1.10. Posterior distribution of the temperature coefficient for the four models. The black line and dot at the bottom of each distribution represent the highest density interval (HDI) and the mean, respectively.

Fig. 1.10, we can see that only the estimated coefficient for Model.1 is biased. As explained in the identification step, to identify the causal effect of indoor air temperature \((\text{T})\) on productivity \((\text{P})\), we need to block the backdoor path (i.e., the non-causal path) \(\text{T} \leftarrow \text{H} \rightarrow \text{P}\). Fitting a separate regression model to the subset of data with H == 'classic' and one to H == 'alternative' is equivalent to conditioning for \(\text{H}\), closing the backdoor path and allowing unbiased estimation of the causal effect of indoor air temperature \((\text{T})\) on productivity \((\text{P})\) (although this comes at the price of a larger standard error and consequently wider credible interval compared to Model.2, where \(\text{H}\) is added as a predictor).

Session info

Version information about R, the OS and attached or loaded packages.

sessionInfo()

R version 4.2.3 (2023-03-15 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 26200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.utf8 
[2] LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] rstan_2.32.7        StanHeaders_2.32.10 gt_1.1.0           
 [4] lubridate_1.9.4     forcats_1.0.1       stringr_1.6.0      
 [7] dplyr_1.1.4         purrr_1.2.0         readr_2.1.6        
[10] tidyr_1.3.1         tibble_3.3.0        ggplot2_4.0.1      
[13] tidyverse_2.0.0     ggdist_3.3.3        tidybayes_3.0.7    
[16] rstanarm_2.32.2     Rcpp_1.1.0          dagitty_0.3-4      
[19] ggdag_0.2.13       

loaded via a namespace (and not attached):
  [1] backports_1.5.0       plyr_1.8.9            igraph_2.2.1         
  [4] splines_4.2.3         svUnit_1.0.8          crosstalk_1.2.2      
  [7] rstantools_2.5.0      inline_0.3.21         digest_0.6.38        
 [10] htmltools_0.5.8.1     viridis_0.6.5         magrittr_2.0.4       
 [13] checkmate_2.3.3       memoise_2.0.1         tzdb_0.5.0           
 [16] graphlayouts_1.2.2    RcppParallel_5.1.11-1 matrixStats_1.5.0    
 [19] xts_0.14.1            timechange_0.3.0      colorspace_2.1-2     
 [22] ggrepel_0.9.6         rbibutils_2.4         xfun_0.54            
 [25] jsonlite_2.0.0        lme4_1.1-37           survival_3.8-3       
 [28] zoo_1.8-14            glue_1.8.0            reformulas_0.4.2     
 [31] polyclip_1.10-7       gtable_0.3.6          V8_8.0.1             
 [34] distributional_0.5.0  car_3.1-3             pkgbuild_1.4.8       
 [37] abind_1.4-8           scales_1.4.0          miniUI_0.1.2         
 [40] viridisLite_0.4.2     xtable_1.8-4          Formula_1.2-5        
 [43] stats4_4.2.3          DT_0.34.0             htmlwidgets_1.6.4    
 [46] threejs_0.3.4         arrayhelpers_1.1-0    RColorBrewer_1.1-3   
 [49] posterior_1.6.1       pkgconfig_2.0.3       loo_2.8.0            
 [52] farver_2.1.2          sass_0.4.10           utf8_1.2.6           
 [55] tidyselect_1.2.1      labeling_0.4.3        rlang_1.1.6          
 [58] reshape2_1.4.5        later_1.4.4           tools_4.2.3          
 [61] cachem_1.1.0          cli_3.6.5             generics_0.1.4       
 [64] evaluate_1.0.5        fastmap_1.2.0         yaml_2.3.10          
 [67] knitr_1.50            fs_1.6.6              tidygraph_1.3.1      
 [70] ggraph_2.2.2          nlme_3.1-168          mime_0.13            
 [73] xml2_1.5.0            compiler_4.2.3        bayesplot_1.14.0     
 [76] shinythemes_1.2.0     rstudioapi_0.17.1     curl_7.0.0           
 [79] tweenr_2.0.3          stringi_1.8.7         lattice_0.22-7       
 [82] Matrix_1.6-5          nloptr_2.2.1          markdown_2.0         
 [85] shinyjs_2.1.0         tensorA_0.36.2.1      vctrs_0.6.5          
 [88] pillar_1.11.1         lifecycle_1.0.4       Rdpack_2.6.4         
 [91] httpuv_1.6.16         QuickJSR_1.8.1        R6_2.6.1             
 [94] promises_1.5.0        gridExtra_2.3         codetools_0.2-20     
 [97] boot_1.3-32           colourpicker_1.3.0    MASS_7.3-58.2        
[100] gtools_3.9.5          withr_3.0.2           shinystan_2.6.0      
[103] parallel_4.2.3        hms_1.1.4             grid_4.2.3           
[106] coda_0.19-4.1         minqa_1.2.8           rmarkdown_2.30       
[109] S7_0.2.1              carData_3.0-5         otel_0.2.0           
[112] ggforce_0.5.0         shiny_1.11.1          base64enc_0.1-3      
[115] dygraphs_1.1.1.6     

Bååth, R., 2025. Bayesboot: An implementation of rubin’s (1981) bayesian bootstrap. https://doi.org/10.32614/CRAN.package.bayesboot

Barrett, M., 2024. Ggdag: Analyze and create elegant directed acyclic graphs.

Fox, J., Weisberg, S., 2019. An R companion to applied regression, Third. ed. Sage, Thousand Oaks CA.

Goodrich, B., Gabry, J., Ali, I., Brilleman, S., 2024. Rstanarm: Bayesian applied regression modeling via Stan.

Iannone, R., Cheng, J., Schloerke, B., Hughes, E., Lauer, A., Seo, J., Brevoort, K., Roy, O., 2025. Gt: Easily create presentation-ready display tables.

Kay, M., 2025. ggdist: Visualizations of distributions and uncertainty. https://doi.org/10.5281/zenodo.3879620

Kay, M., 2024b. tidybayes: Tidy data and geoms for Bayesian models. https://doi.org/10.5281/zenodo.1308151

Kay, M., 2024a. ggdist: Visualizations of distributions and uncertainty in the grammar of graphics. IEEE Transactions on Visualization and Computer Graphics 30, 414–424. https://doi.org/10.1109/TVCG.2023.3327195

Posit team, 2025. RStudio: Integrated development environment for r. Posit Software, PBC, Boston, MA.

R Core Team, 2023. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Stan Development Team, 2025. RStan: The R interface to Stan.

Textor, J., van der Zander, B., Gilthorpe, M.S., Liśkiewicz, M., Ellison, G.T., 2016. Robust causal inference using directed acyclic graphs: The R package ’dagitty’. International Journal of Epidemiology 45, 1887–1894. https://doi.org/10.1093/ije/dyw341

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L.D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T.L., Miller, E., Bache, S.M., Müller, K., Ooms, J., Robinson, D., Seidel, D.P., Spinu, V., Takahashi, K., Vaughan, D., Wilke, C., Woo, K., Yutani, H., 2019. Welcome to the tidyverse. Journal of Open Source Software 4, 1686. https://doi.org/10.21105/joss.01686