The perfect doctor

Lab specs can be found here

Exercise 1

Question: Recreate table 1 shown below, inserting values appropriately for the three empty colums: (i) The column labeled \(\delta_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.

# Write your code here

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) %>% 
  gt()
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
 
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) -> data

Answer: 1. ATE:

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

  1. The treatment ventilators is more effective than the bedrest in the ventilators-treated group, while the bedrest is more effective in the bedrest-treated (control) group.

  2. If there are multiple versions of so-called ventilators or bedrest and each version could lead to a different potential outcome, the SUTVA would be violated. Or if the patients are in a same hospital room and encourage each other, then the potential outcome would be affected by other patients.

Exercise 2

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]\).

# Write your code here
## SDO:
data %>% 
  group_by(D) %>% 
  reframe(mean = mean(Y)) %>% 
  summarise(SDO = mean[2] - mean[1])
## selection bias:
data %>% 
  group_by(D) %>% 
  reframe(mean = mean(`Y(0)`)) %>% 
  summarise(bias = mean[2] - mean[1])

Answer: The SDO equals 0.86 while the ATT equals 4, which are different. The selection bias between two groups is -0.314. i.e., SDO = ATT + bias.

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?

# Write your code here
## ATE:
data %>% 
  summarise(ATE = mean(delta))
## ATT:
data %>% 
  filter(D == 1) %>% 
  summarise(ATT = mean(delta))
## ATU:
data %>% 
  filter(D == 0) %>% 
  summarise(ATU = mean(delta))
## ATO:
data %>% 
  filter(Age>=50 & Age <=77) %>% 
  summarise(ATO = mean(delta))

Answer: Mostly, age modifies the treatment effect in each group. Patients in the ventilators-treated group aged older than those in the bedrest-treated group. The overall ATE is than closed to zero because it is (weighted) average across positive and negative effect.

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 \(\mathbb{E}[Y|D=0]\). Now estimate the following regression, comparing the coefficients \(\alpha\) and \(\delta\) to the statistics you’ve previously calculated. What did you find? How would you explain these finding?

# Write your code here
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
# mean outcome in the control group
data %>% 
  filter(D == 0) %>% 
  summarise(mean(Y))
# mean difference between two group
data %>% 
  group_by(D) %>% 
  summarise(mean = mean(Y)) %>% 
  summarise(naive = mean[2] - mean[1])

Answer: The coefficient alpha/intercept is the mean observed outcome in the control group. The coefficient \(\delta\) is the average treatment effect difference on the observed outcome. The model-based mean outcome and the unadjusted average treatment effect in the control group is exactly the same as those from the naive calculation. Because the model is saturated with only one binary variable, it is mathematically equivalence to the naive calculation which can be considered as non-parametric MLE.

Exercise 5

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

# Write your code here
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: When controlling for age in the regression model, the coefficient of treatment is more closed to zero. If age is indeed a confounder, then the adjusted coefficient is more likely be a causal estimate.

Exercise 6

Question: 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 \(X_{age}\). The third equation regresses \(Y\) onto \(\tilde{D}\) which is the residual from the second equation. Compare the coefficient on \(D\) from the first equation to the coefficient on \(\tilde{D}\) in the third equation. What does this tell you about how to interpret multivariate regressions?

# Write your code here
# 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
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 on D from first model and third model are exactly same. Because the second regression build a independence between treatment D and covariate age, and the residual of second regression is the residual variation of treatment D conditional on age, the third regression behaves same as the first regression which explicitly conditional on age in the regression formula. I think this kind of technique is so-called orthogonalized regression.