Lab01-Rubins potential outcomes framework

Author

Yiling Zhou

The perfect doctor

Table 1 below displays a small data set that is used for practice purposes only. Two unobserved (and imagined) potential outcomes are recorded for each patient, denoting years of post-treatment survival under each of two treatments. Suppose the “perfect doctor” knows each patient’s potential outcomes and as a result chooses the best treatment for each patient. She admisters bedrest to patients who would benefit from it most (D=0𝐷=0) and ventilators to those who would benefit from it more (D=1𝐷=1).

Exercise1

Question: Recreate table 1 shown below, inserting values appropriately for the three empty colums: (i) The column labeled δi𝛿𝑖: please enter each patient’s treatment effect (ii) The column labeled D𝐷: the optimal treatment for this patient (iii) The column labeled Y𝑌: the observed outcomes. Calculate the average treatment effect (ATE) and the average treatment effect for the treated (ATT) when comparing the outcome of the ventilators treatment with that of the bedrest treatment and comment as to which type of intervention is more effective on average. Finally, explain under which conditions might SUTVA be violated for treatments of covid-19 in the scenario described above.

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.3     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.3     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

data <- tribble(
  ~patient, ~`Y(0)`, ~`Y(1)`, ~Age, ~`delta`, ~`D`, ~`Y`,
  1, 10, 1,  29, -9, 0, 10,
  2,  5, 1,  35, -4, 0, 5,
  3,  4, 1,  19, -3, 0, 4,
  4,  6, 5,  45, -1, 0, 6,
  5,  1, 5,  65, 4, 1, 5,
  6,  7, 6,  50, -1, 0, 7, 
  7,  8, 7,  77, -1, 0, 8, 
  8, 10, 7,  18, -3, 0, 10, 
  9,  2, 8,  85, 6, 1, 8,
  10, 6, 9,  96, 3, 1, 9,
  11, 7, 10, 77, 3, 1, 10) 

print(data)

# A tibble: 11 × 7
   patient `Y(0)` `Y(1)`   Age delta     D     Y
     <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
 1       1     10      1    29    -9     0    10
 2       2      5      1    35    -4     0     5
 3       3      4      1    19    -3     0     4
 4       4      6      5    45    -1     0     6
 5       5      1      5    65     4     1     5
 6       6      7      6    50    -1     0     7
 7       7      8      7    77    -1     0     8
 8       8     10      7    18    -3     0    10
 9       9      2      8    85     6     1     8
10      10      6      9    96     3     1     9
11      11      7     10    77     3     1    10

data %>% 
  summarise(ATE = mean(delta))

# A tibble: 1 × 1
     ATE
   <dbl>
1 -0.545

data %>% 
  filter(D == 1) %>% 
  summarise(ATT = mean(delta))

# A tibble: 1 × 1
    ATT
  <dbl>
1     4

Answer: In total population, bedrest is more effective, ATE is -0.545. For treated population, ventilators is more effective, ATT >0.

If patients are in the same hospital room and their interactions or encouragement of one another affect their outcomes, the SUTVA would also be violated, as an individual’s potential outcome is influenced by the presence or characteristics of other individuals in the same setting. Or if the individual potential outcome varies along with the time or the type of the drugs.

Exercise2

Question: Calculate the simple difference in outcomes (SDO), showing the details of your calculation. Is the SDO a good estimation for the ATE? Finally, check whether the SDO is equal to the sum of the ATT and the selection bias, E[Y(0)|T=1]−E[Y(0)|T=0]𝐸[𝑌(0)|𝑇=1]−𝐸[𝑌(0)|𝑇=0].

ATE <- data %>% 
  summarise(ATE = mean(delta))
ATE$ATE

[1] -0.5454545

ATT <- data %>% 
  filter(D == 1) %>% 
  summarise(ATT = mean(delta))
ATT$ATT

[1] 4

ObservedY_E <- data %>% 
  group_by(D) %>% 
  summarise(EY = mean(Y))
SDO <-  ObservedY_E$EY[ObservedY_E$D==1] - ObservedY_E$EY[ObservedY_E$D==0]
SDO

[1] 0.8571429

selection <- data %>% 
  group_by(D) %>% 
  summarise(EY0 = mean(`Y(0)`))
selection_bias <- selection$EY0[2] - selection$EY0[1]
selection_bias

[1] -3.142857

difference_ATE_SDO <- SDO-ATE$ATE
difference_ATE_SDO

[1] 1.402597

difference_ATT_SDO <- SDO-ATT$ATT
difference_ATT_SDO

[1] -3.142857

Answer: SDO is not a good measure for ATE, difference between ATE (-0.545) and SDO (0.857) is 1.402597; SDO is equal to ATT plus selection bias: selection_bias = difference_ATT_SDO.

Exercise 3

Question: Compare the treatment effect for both groups: for those treated with a ventilator and for those treated with bedrest. What explains the difference in the average effect? Now compare all four measures of effects. What are the advantages and disadvantages of each? Is the ATE equal to the mean of the ATU and the ATT? Why or why not?

ATE <- data %>% 
  summarise(ATE = mean(delta))
ATE$ATE

[1] -0.5454545

ATT <- data %>% 
  filter(D == 1) %>% 
  summarise(ATT = mean(delta))
ATT$ATT

[1] 4

ATU <- data %>% 
  filter(D == 0) %>% 
  summarise(ATU = mean(delta))
ATU$ATU

[1] -3.142857

ObservedY_E <- data %>% 
  group_by(D) %>% 
  summarise(EY = mean(Y))
SDO <-  ObservedY_E$EY[ObservedY_E$D==1] - ObservedY_E$EY[ObservedY_E$D==0]
SDO

[1] 0.8571429

data %>% 
  group_by(D) %>% 
  summarise(mean(Age))

# A tibble: 2 × 2
      D `mean(Age)`
  <dbl>       <dbl>
1     0        39  
2     1        80.8

###whether ATE = average(ATT, ATU)
P_D_0 <- sum(data$D==0)/nrow(data)
P_D_1 <- sum(data$D==1)/nrow(data)
average <- ATT$ATT *P_D_1 + ATU$ATU*P_D_0
average

[1] -0.5454545

###ATE =  mean (ATT, ATU)
(ATT$ATT  + ATU$ATU)/2

[1] 0.4285714

Answer: ATT (in only treated group) is not equal to ATE (in all population, average of ATT and ATU), because untreated group (D=0) and treated group (D=1) are different, which means treated group is different from the total population, and this is supported by the selection bias (3.142857), i.e., the potential outcome if untreated is different between untreated group and treated group. Age in the treated group (80.75) is higher than that in untreated group (39). Age is more likely to be confounder.

If there is a perfect doctor, it is better to know ATT, ATU, rather than ATE. ATT is average treatment effect in group (D=1), which means this subpopulation will benefit more from ventilators. ATU is average treatment effect in group (D=0), which means this subpopulation will benefit more from bedrest. However, ATE is better if there is no perfect doctor.

If we calculate the mean of the ATU and the ATT using (ATU+ATT)/2, then it is not equal to ATE, because the treated group and untreated group have different portion in the total population. Thus, the weighted mean of ATU and ATT, average=ATU*P(D=0) + ATT*P(D=1) will be equal to ATE.

Using regression to estimate effects

The following exercises demonstrate that regression is a useful tool to estimate average outcomes and treatment effects in the different groups. Notice that in contrast to the role that regressions play in traditional statistics, here standard errors and significance are not of primary concern. Instead, we are interested in using regression to calculate average effects proper.

Exercise 4

Question: Calculate the outcome, conditional on getting the bedrest treatment 𝔼[Y|D=0]𝐸[𝑌|𝐷=0]. Now estimate the following regression, comparing the coefficients α𝛼 and δ𝛿 to the statistics you’ve previously calculated. What did you find? How would you explain these finding?

fit <- lm(Y ~ D, data = data)
fit %>% summary()


Call:
lm(formula = Y ~ D, data = data)

Residuals:
   Min     1Q Median     3Q    Max 
-3.143 -1.643  0.000  1.500  2.857 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.1429     0.8624   8.282 1.68e-05 ***
D             0.8571     1.4302   0.599    0.564    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.282 on 9 degrees of freedom
Multiple R-squared:  0.03838,   Adjusted R-squared:  -0.06847 
F-statistic: 0.3592 on 1 and 9 DF,  p-value: 0.5637

## Naive: the observed average ourtcome in control group
data %>% 
  filter(D==0) %>% 
  summarise(mean(Y))

# A tibble: 1 × 1
  `mean(Y)`
      <dbl>
1      7.14

Answer：The alpha/intercept is the mean observed outcome in the control group. The coefficient of D is the SDO. The intercept from the model and the observed average outcome in the control group from previous calculation (Naive) are exactly the same. The coefficient and the observed average treatment effect, SDO (Naive) are the same.

Exercise 5

Question: Now estimate the same regression, but this time, controlling for age, again comparing, the coefficient δ𝛿 to the statistics you’ve previously calculated. What did you find? How do you explain these results?

fit <- lm(Y ~ D + Age, data = data)

fit %>% summary()


Call:
lm(formula = Y ~ D + Age, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7390 -1.6630 -0.0858  1.3839  3.2811 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  6.35543    1.90736   3.332   0.0104 *
D            0.01419    2.33997   0.006   0.9953  
Age          0.02019    0.04309   0.469   0.6519  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.388 on 8 degrees of freedom
Multiple R-squared:  0.06407,   Adjusted R-squared:  -0.1699 
F-statistic: 0.2738 on 2 and 8 DF,  p-value: 0.7673

Answer: After adjusting age in the regression model, the coefficient of treatment is 0.014. If age is a confounder, the adjusted coefficient is more likely be a causal estimate.

Exercise 6

Estimate the following three regression models. The first model is the same as the one above. The second equation is the auxiliary regression of D𝐷 onto Xage𝑋𝑎𝑔𝑒. The third equation regresses Y𝑌 onto D̃ 𝐷~ which is the residual from the second equation. Compare the coefficient on D𝐷 from the first equation to the coefficient on D̃ 𝐷~ in the third equation. What does this tell you about how to interpret multivariate regressions?

# first model
fit1 <- lm(Y ~ D + Age, data = data)
# second model
fit2 <- lm(D ~ Age, data = data)
# extract residual
data$residuals <- fit2$residuals
# third model
fit3 <- lm(Y ~ residuals, data = data)
fit1 %>% summary()


Call:
lm(formula = Y ~ D + Age, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7390 -1.6630 -0.0858  1.3839  3.2811 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)  6.35543    1.90736   3.332   0.0104 *
D            0.01419    2.33997   0.006   0.9953  
Age          0.02019    0.04309   0.469   0.6519  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.388 on 8 degrees of freedom
Multiple R-squared:  0.06407,   Adjusted R-squared:  -0.1699 
F-statistic: 0.2738 on 2 and 8 DF,  p-value: 0.7673

fit2 %>% summary()


Call:
lm(formula = D ~ Age, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.68663 -0.16289  0.04442  0.17433  0.48323 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.403323   0.236116  -1.708  0.12178   
Age          0.014155   0.003925   3.606  0.00569 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3401 on 9 degrees of freedom
Multiple R-squared:  0.591, Adjusted R-squared:  0.5455 
F-statistic:    13 on 1 and 9 DF,  p-value: 0.005694

fit3 %>% summary()


Call:
lm(formula = Y ~ residuals, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4565 -1.9522  0.5426  2.0429  2.5456 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.45455    0.70156  10.626 2.15e-06 ***
residuals    0.01419    2.28040   0.006    0.995    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.327 on 9 degrees of freedom
Multiple R-squared:  4.302e-06, Adjusted R-squared:  -0.1111 
F-statistic: 3.872e-05 on 1 and 9 DF,  p-value: 0.9952

Answer: The coefficients from the first model and the third model are same. This is because the residuals from second regression is the part of D independent of the covariate age, and the residuals from the second regression represent the residual variation of D after accounting for age. Consequently, the third regression, which includes the residuals from the second regression, effectively conditions on age, behaving the same way as the first regression that explicitly conditions on age.