Han

Causal Inference

The role of propensity scores in observational study

Han

十月 2022

library(arm)
library(ggplot2)
library(rlang)
library(RColorBrewer)
library(stats)
library(Hmisc)
library(caret)
library(randomForest)
library(nnet)
library(dplyr)
library(MatchIt)
library(reshape2)
library(knitr)
library(kableExtra)

# Set global ggplot theme
theme_set(theme_minimal())
# Define global color palettes
col.RdBl.2 <- brewer.pal(3, 'RdBu')[-2]

Objective

This assignment will give you the opportunity to practice several different propensity score approaches to causal inference. In addition you will be asked to interpret the resulting output and discuss the assumptions necessary for causal inference.

Problem Statement

In this assignment will use data from a constructed observational study. The data and an associated data dictionary are available in this folder. The treatment group for the study that the data are drawn from is the group of children who participated in the IHDP intervention discussed in class. The research question of interest focuses on the effect of the IHDP intervention on age 3 IQ scores for the children that participated in it. The data for the comparison sample of children was pulled from the National Longitudinal Study of Youth during a similar period of time that the data were collected for the IHDP study.

Question 1: Load the data and choose confounders (Step 1)

A = read.csv('/Users/zhanghan/R_projects/ihdp-causality-master/ihdp.csv')


corr_mat <- cor(A)
cat(dim(A))

## 747 27

corr_mat = rbind(corr_mat['treatment',], corr_mat['outcome',])
confounders = names(which((abs(corr_mat[1,]) > 0.03) * (abs(corr_mat[2,])>0.02) >0))
dt = A[confounders]
confounders = confounders[!(confounders %in% c('treatment', 'outcome'))]
cat('筛选后的长度',length(confounders))

## 筛选后的长度 15

# Subset data for analysis:
#  (a) Reduce data frame to include only observations for children whose b.w. < 3000 g
#  (b) Include outcome, treatment indicator, and selected covariates

# New data frame with outcome in 1st column, treatment indicator in 2nd column
# and covariates in the remaining columns
head(dt)

# List of the variable names for the confounder variables chosen
colnames(dt)[3:length(dt)]

##  [1] "feature0"  "feature1"  "feature3"  "feature4"  "feature5"  "feature8"  "feature9"  "feature13"
##  [9] "feature15" "feature16" "feature17" "feature21" "feature22" "feature23" "feature24"

---

Question 2: Estimate the propensity score (Step 2)

ps.m1 <- glm(treatment ~ ., data=dt[, c('treatment',confounders)], family=binomial(link='logit'))
dt$psc <- ps.m1$fitted.values

---

Question 3: Restructure your data through matching. [Or at least create the weights variable that will let you to do so in the following steps] (Step 3)

(a) The first thing you need to be clear on before restructuring your data is the estimand. Given the description above about the research question, what is the estimand of interest?

The research question of interest focuses on the effect of the intervention on the children that participated in it. Since we are interested on the effect of those who received treatment, the estimand of interest is the (average treatment effect for the treated).

(b) First please perform one-to-one nearest neighbor matching with replacement using your estimated propensity score from Question 2. Perform this matching using the matching command in the arm package. The “cnts” variable in the output reflects the number of times each control observation was used as a match (the length is equal to the number of control observations). Use the output of this function to create a weight variable that

• 1. equals one for treated observations and
     
• 2. equals the number of times used as a match for non- treated observations.

matches <- matching(z=dt$treatment, score=dt$psc, replace=FALSE)
dt[dt$treat==0, 'wt'] <- matches$cnts[dt$treatment==0]
dt[dt$treat==1, 'wt'] <- 1

---

Question 4: Check overlap and balance. (Step 4)

(a) Examining Overlap. Check overlap on the unmatched data using some diagnostic plots. Check overlap for the propensity scores as well as two other covariates.

# Paired, inverted histograms: propensity scores
examine_overlap = function(dt, nbins=40, var='psc', var_name='Propensity Score', xlab='Propensity Score', ylim=c(-30,50), options=NULL){

    dt_plot = dt[, c('treatment', var)]
    colnames(dt_plot) = c('treatment', 'X')
    
    ggplot(dt_plot) +
        geom_histogram(data=dt_plot[dt_plot$treat==0,], bins=nbins, fill='grey', alpha=0.1, aes(X, y=..count.., color=col.RdBl.2[2])) +
        geom_histogram(data=dt_plot[dt_plot$treat==1,], bins=nbins, fill='grey', alpha=0.1, aes(X, y=-..count.., color=col.RdBl.2[1])) +
        coord_cartesian(ylim=ylim) +
        scale_color_manual(values=rev(col.RdBl.2), labels=c('Control', 'Treated'), name='') +
        labs(x=xlab, y='Frequency', title=paste('Overlap of ', var_name, ' between Groups', sep='')) +
        options 
}

# Overlap of propensity score
examine_overlap(dt, nbins=40, var='psc', var_name='Propensity Score', xlab='Propensity Score', ylim=c(-10,10))

# Overlap of preterm
examine_overlap(dt, nbins=20, var='feature0', var_name='Number Count', xlab='Feature0', ylim=c(-50,50))

(b) Interpreting Overlap. What do these plots reveal about the overlap required to estimate our estimand of interest.

These plots reveal we have enough overlap to estimate the average treatment effect on the treated. Between the treatment and control groups, we examined the overlap in propensity scores, preterm, and birth weight covariates. A few notes.
One, due to the imbalance in number of observations between treatment and control groups, comparing paired histograms can misleadingly imply lack of overlap. We addressed this issue by limiting the range of y-axis in the plots, effectively “zooming in”, so that all histogram bins are clearly visible.
Two, notice in the paired histograms for the propensity score that there are no control observations in the highest bin group of propensity scores (0.975-1.00). This may suggest there is not complete overlap, however, it is also an artifact of the histogram bin size. Doubling the bin size from 0.025 to 0.05 propensity score units eliminates this lack of overlap.
Lastly, there is not complete overlap in birth weight between the treatment and control groups. However, the complete overlap in the propensity score indicates that in the matched sample there should be overlap between groups across the multi-dimensional covariate space. Thus we are unconcerned with this lack of complete overlap in birth weight.

(c) Examining Balance. You will build your own function to check balance! This function should take as inputs the data frame created in Question 1, the vector with the covariate names chosen in Question 1, and the weights created in Question 2. It should output the following:

• 1. Mean in the unmatched treatment group
     
• 2. Mean in the unmatched control group
     
• 3. Mean in the matched treatment group
     
• 4. Mean in the matched control group
     
• 5. Unmatched mean difference (standardized for continuous variables, not standardized for binary variables)
     
• 6. Matched mean difference (standardized for continuous variables, not standardized for binary variables)
     
• 7. Ratio of standard deviations across unmatched groups (control/treated)
     
• 8. Ratio of standard deviations across matched groups (control/treated)

checkBalance = function(df, confounders, wt_var) {

    # Subset data to confounders only
    df2 <- df[, confounders]

    binary = apply(df2, 2, function(x) all(x %in% 0:1)) # Binary Variable Indicator
    trt = df$treatment == 1                   # Treatment Indicator
    ctr = df$treatment == 0                   # Control Indicator
    n_trt = sum(trt)                      # Treatment Sample Size
    n_ctrl = sum(ctr)                     # Control Sample Size
    wts.trt = df[trt, wt_var]              # Set treatment weights
    wts.ctr = df[ctr, wt_var]  # Set control weights
  
    # Means
    trt.means = apply(df2[trt,], 2, mean)
    trt.means_w = apply(df2[trt,], 2, wtd.mean, wts.trt)
    ctr.means = apply(df2[ctr,], 2, mean)
    ctr.means_w = apply(df2[ctr,], 2, wtd.mean, wts.ctr)
  
    # Variances
    trt.var = apply(df2[trt,], 2, var)
    trt.var_w = apply(df2[trt,], 2, function(x) sum(wts.trt * (x - wtd.mean(x, wts.trt))^2) / (sum(wts.trt) - 1))
    ctr.var = apply(df2[ctr,], 2, var)
    ctr.var_w = apply(df2[ctr,], 2, function(x) sum(wts.ctr * (x - wtd.mean(x, wts.ctr))^2) / (sum(wts.ctr) - 1))
  
    # Standardized Mean Differences
    mean.diff = (trt.means-ctr.means) / sqrt(trt.var)
    # 处理离散变量
    mean.diff.bin = (trt.means-ctr.means)
    diff = mean.diff * (1-binary) + mean.diff.bin * binary
  
    mean.diff_w = (trt.means_w-ctr.means_w) / sqrt(trt.var_w)
    mean.diff.bin_w = (trt.means_w-ctr.means_w)
    diff.m = mean.diff_w*(1-binary) + mean.diff.bin_w*binary
  
    # Ratios of standard deviations
    ratio = sqrt(ctr.var)/sqrt(trt.var)
    ratio[binary] = NA

    ratio.m = sqrt(ctr.var_w)/sqrt(trt.var_w)
    ratio.m[binary] = NA
  
    # Return results
    result = cbind(trt.means, ctr.means, trt.means_w, ctr.means_w, diff, diff.m, ratio, ratio.m)
    colnames(result) = c('mean_t', 'mean_c', 'mean_t_match', 'mean_c_match', 'diff', 'diff_m', 'ratio', 'ratio_m')
    return(result)
}

(d) How do you interpret the resulting balance? In particular what are your concerns with regard to covariates that are not well balanced (write about 5 or 6 sentences).

Given we are using propensity scores to match controls observations as empirical counterfactuals to treatment observations, we wish to see better balance in the matched sample compared to the unmatched sample. Specifically, we are looking to reduce mean differences and achieve ratios of standard deviations close to 1 in the matched sample.
First, we examine the difference in means (standardized for continuous covariates) in the matched sample, and notice similar group means (mean difference < 0.10) in nearly half of the covariates (7 of 15). The fact that the means of several mother-related covariates (momage, b.marr, first, black, hispanic) are dissimilar between the treated and control groups in the matched sample is of concern. This imbalance in means between groups could increase bias in treatment effect estimation. Additionally, the ratio of standard deviations in the matched sample is only close to 1.0 for 1 of 3 continuous covariates. The greater variance in mother’s age at birth in the treatment group compared to the control group is also concerning for the same reasons of possibly increasing bias in treatment effect estimation.

ps.m1.bal <- checkBalance(df=dt, confounders=confounders, wt_var='wt')

# Rounded to two decimal places
round(ps.m1.bal,2)

##           mean_t mean_c mean_t_match mean_c_match  diff diff_m ratio ratio_m
## feature0    0.21  -0.05         0.21         0.15  0.28   0.06  1.10    1.02
## feature1    0.18  -0.04         0.18         0.21  0.25  -0.04  1.15    1.08
## feature3   -0.22   0.05        -0.22        -0.21 -0.30  -0.02  1.12    1.06
## feature4   -0.14   0.03        -0.14        -0.23 -0.17   0.08  0.95    1.02
## feature5    0.21  -0.05         0.21         0.34  0.28  -0.13  1.09    1.15
## feature8    0.68   0.49         0.68         0.68  0.19   0.00    NA      NA
## feature9    0.29   0.38         0.29         0.32 -0.09  -0.04    NA      NA
## feature13   1.58   1.44         1.58         1.60  0.28  -0.06  1.00    0.99
## feature15   0.94   0.97         0.94         0.96 -0.03  -0.03    NA      NA
## feature16   0.69   0.57         0.69         0.70  0.12  -0.01    NA      NA
## feature17   0.99   0.96         0.99         0.99  0.04   0.00    NA      NA
## feature21   0.04   0.09         0.04         0.04 -0.06   0.00    NA      NA
## feature22   0.01   0.09         0.01         0.00 -0.07   0.01    NA      NA
## feature23   0.06   0.14         0.06         0.07 -0.09  -0.01    NA      NA
## feature24   0.27   0.13         0.27         0.22  0.14   0.05    NA      NA

能够观察到在配对后样本更加均衡。均值差异和方差比例更加合理，说明倾向性得分确实能够对于样本进行匹配。

# Plot standardized mean Differences
plotMeanDiff = function(balance_object, name){
  df = data.frame(balance_object[,c(5,6)])
  df$var = as.factor(row.names(df))
  colnames(df) = c('unmatched','matched','var')
  # knitr::kable(df)
  df = melt(df)
  # knitr::kable(df)
  ggplot(df, aes(x = value, y = var)) +
    geom_point(aes(shape = variable, color=variable)) +
    geom_vline(xintercept=0, color='blue') +
    geom_vline(xintercept=-.05, color='orange', linetype='dotted') +
    geom_vline(xintercept=.05, color='orange', linetype='dotted') +
    geom_vline(xintercept=-.1, color='red', linetype='dotted') +
    geom_vline(xintercept=.1, color='red', linetype='dotted') +
    coord_cartesian(xlim=c(-0.5,0.5)) +
    labs(x = 'Standardized Difference in Means', y = 'predictor', title = name,
         caption='Orange dotted line indicates +/-0.05 bounds and red lines indicate +/-0.1 bounds.\n Values outside of range x=[-0.5,0.5] not displayed.')
}

# Plot Ratios of Standard Deviations
plotSDRatios = function(balance_object,name){
  df = data.frame(balance_object[,c(7,8)])
  df = df[!is.na(df)[,1],]
  df$var = as.factor(row.names(df))
  colnames(df) = c('unmatched','matched','var')
  df = melt(df)

  ggplot(df, aes(x = value, y = var, shape = variable, color = variable)) + 
    geom_point() + 
    geom_vline(xintercept=1, color='blue') + 
    labs(x = 'Ratios of Standard Deviations', y = 'predictor', title=name)
}

plotMeanDiff(ps.m1.bal, 'Propensity Score Model: Logistic Regression')

## Using var as id variables

plotSDRatios(ps.m1.bal, 'Propensity Score Model: Logistic Regression')

## Using var as id variables

(e) Unit test. 子集回归

# Subset data
temp <- A[c('treatment', 'feature4', 'feature5')]

# Fit propensity score model (logistic)
ps.temp <- glm(treatment ~ feature4 + feature5, data=temp, family=binomial(link='logit'))

# Generate propensities scores
temp$psc <- ps.temp$fitted.values

# 1-1 nearest neighbor matching with replacement
ps.temp.matches <- matching(z=temp$treatment, score=temp$psc, replace=TRUE)
temp[temp$treatment==0, 'wt'] <- ps.temp.matches$cnts[temp$treatment==0]
temp[temp$treatment==1, 'wt'] <- 1

# Rounded to 3 decimal places
round(checkBalance(temp, c('feature4', 'feature5'), 'wt'),3)

##          mean_t mean_c mean_t_match mean_c_match   diff diff_m ratio ratio_m
## feature4 -0.142  0.032       -0.142       -0.212 -0.168  0.067 0.952   0.868
## feature5  0.212 -0.049        0.212        0.172  0.280  0.044 1.086   0.967

---

Question 5: Repeat steps 2-4 within the matching framework.

We use the following propensity score models to try to achieve better balance:
- Probit regression

We use the following matching methods to try to achieve better balance:
- Mahalanobis Distance Matching
- Matching without Replacement (Using propensity scores from original logistic regression)
- Nearest Neighbor Matching with Replacement

Brief descriptions of each of the methods:

(1) Probit Regression

A probit regression is equivalent to the logistic regression except that it replaces the logistic with the normal distribution:

\[Pr(y_{i} = 1) = \Phi(X_{i}B)\]

Where \(\Phi\) is the normal cumulative distribution function. There is strong overlap between the control and treatment distributions over the propensity score so no observations were dropped. We use matching with replacement.

(2) Mahalanobis Matching

Mahalanobis matching uses a distance measure in multivariate space that takes into account variances of variables as well as covariances between them. We match, with replacement, based on the calculated Mahalanobis distance.

(3) Matching without Replacement (Using propensity scores from original logistic regression)

We use the original propensity scores from Question 2. This time, we match without replacement.

(4) Nearest Neighbor Matching with Replacement

We use K-nearest neighbor matching with k = 2. This matching method selects the 2 control subjects with the closest distance to the treated subject (using default=logit). Matches are chosen for treatment subjects one at a time using the default order (largest to smallest). We allow for replacement.

Fit new propensity score models

# Model 2: Probit
ps.m2 <- glm(treatment ~ ., data=dt[, c('treatment', confounders)], family=binomial(link='probit'))
dt$psc_probit <- ps.m2$fitted.values
confusionMatrix(factor(ifelse(dt$psc_probit>0.5,1,0)), factor(dt$treatment))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   0   1
##          0 605 136
##          1   3   3
##                                           
##                Accuracy : 0.8139          
##                  95% CI : (0.7841, 0.8412)
##     No Information Rate : 0.8139          
##     P-Value [Acc > NIR] : 0.5227          
##                                           
##                   Kappa : 0.0264          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.99507         
##             Specificity : 0.02158         
##          Pos Pred Value : 0.81646         
##          Neg Pred Value : 0.50000         
##              Prevalence : 0.81392         
##          Detection Rate : 0.80991         
##    Detection Prevalence : 0.99197         
##       Balanced Accuracy : 0.50832         
##                                           
##        'Positive' Class : 0               
## 

Matching with new models

# Rematching

# Probit
matches <- matching(z=dt$treatment, score=dt$psc_probit, replace=TRUE)
dt[dt$treatment==1, 'w_probit'] <- 1
dt[dt$treatment==0, 'w_probit'] <- matches$cnts[dt$treatment==0]

# Matching without Replacement using Original Logit Model
matches_nr <- matching(z=dt$treatment, score=dt$psc, replace=FALSE)
# matches_nr_idx <- c(which(dt$treatment==1))
# dt[matches_nr_idx, 'w_logit_no_rep'] = 1
dt[dt$treatment==1, 'w_logit_no_rep'] = 1
dt[dt$treatment==0, 'w_logit_no_rep'] = matches_nr$cnts[dt$treatment==0]

# dt$w_logit_no_rep[is.na(dt$w_logit_no_rep)] = 0

Calculate the Mahalanobis distance for matching: \[d(x,y) = \sqrt{(x - y)^T S^{-1} (x-y)}\] We constructed our own Mahalanobis Distance function. As shown below, it gets the same results as the MatchIt package.

# Mahalanobis Distance function, d(x,y)
myMH = function(trt, ctr, inv.cov, data){
  x = as.matrix(data[trt,])
  y = as.matrix(data[ctr,])
  diff = x - y
  sqrt(rowSums((diff %*% inv.cov) * diff))
}

inv_cov_mat = solve(cov(dt[,confounders]))
trt = row.names(dt[dt$treat==1,])
ctr = row.names(dt[dt$treat==0,])
cat('trt长度', length(trt))

## trt长度 139

# Compute distances for all pairwise distances between treated and control
mahalo_dist = outer(trt, ctr, FUN = myMH, inv.cov = inv_cov_mat, data = dt[,confounders])

# Find Matches (Smallest Distance) 
matches = apply(mahalo_dist,1, function(x) min(which(x == min(x))))
cat('matches dim',length(matches))

## matches dim 139

# knitr::kable(matches)
# Get Weights
m_freq = data.frame(table(matches))
# cat('长度',dim(m_freq))

ctr.idx = data.frame(idx=seq(1:sum(dt$treat==0)))
matches = merge(ctr.idx,m_freq,by.x='idx',by.y='matches',all.x=T)

# Set Weights
dt[dt$treat==0, 'wt_mh'] = matches$Freq*0.8129496
dt[is.na(dt$wt_mh), 'wt_mh'] = 0
dt[dt$treat==1, 'wt_mh'] = 1

# Check against MatchIt function  
df = dt[,c('treatment',confounders)]
zz <- matchit(treatment ~ ., data=df, method="nearest",
              distance="mahalanobis", replace=TRUE)

# Check weights derived from our Mahalanobis matches to the MatchIt matches (up to 5 decimal places)
all(round(dt$wt_mh, 5) == round(zz$weights, 5))

## [1] TRUE

KNN (K=2) Matching with Replacement

# KNN (K=2) 
df = dt[,c('treatment',confounders)]
k2 <- matchit(treatment ~ ., data=df, method="nearest",
              ratio=2, replace=T)
dt$wt_k2 = k2$weights

Check overlap of matched sample

title <- ggtitle('Propensity Score Model: Probit Regression')
examine_overlap(dt[dt$w_probit!=0,], nbins=40, var='psc_probit', var_name='Propensity Score', xlab='Propensity Score', ylim=c(-10,10), title)

title <- ggtitle('Matching Method: Logistic w/o Replacement')
examine_overlap(dt[dt$w_logit_no_rep!=0,], nbins=40, var='psc', var_name='Propensity Score', xlab='Propensity Score', ylim=c(-10,10), title)

title <- ggtitle('Matching Method: Mahalanobis')
examine_overlap(dt[dt$wt_mh!=0,], nbins=20, var='feature4', var_name='Number Weeks Baby Born Preterm', xlab='Weeks Preterm', ylim=c(-50,50), title)

title <- ggtitle('Matching Method: Nearest Neighbor (k=2)')
examine_overlap(dt[dt$wt_k2!=0,], nbins=20, var='feature4', var_name='Number Weeks Baby Born Preterm', xlab='Weeks Preterm', ylim=c(-50,20), title)

Check balance

logit.bal <- checkBalance(dt, confounders, 'wt')
probit.bal <- checkBalance(dt, confounders, 'w_probit')
logit_nr.bal <- checkBalance(dt, confounders, 'w_logit_no_rep')
mh.bal <- checkBalance(dt, confounders, 'wt_mh')
mk2.bal <- checkBalance(dt, confounders, 'wt_k2')

plotMeanDiff(logit.bal, 'Propensity Score Model: Logistic Regression')

## Using var as id variables

plotMeanDiff(probit.bal, 'Propensity Score Model: Probit Regression')

## Using var as id variables

plotMeanDiff(logit_nr.bal, 'Matching Method: Logistic w/o Replacement')

## Using var as id variables

plotMeanDiff(mh.bal, 'Matching Method: Mahalanobis')

## Using var as id variables

plotMeanDiff(mk2.bal, 'Matching Method: Nearest Neighbor (k=2)')

## Using var as id variables

# Plot SD Ratios
temp = rbind(logit.bal, probit.bal, logit_nr.bal, mh.bal, mk2.bal)

temp = temp[order(rownames(temp)),]
row.names(temp) = paste(row.names(temp),rep(c('logit','probit','logit_no_replace','mahalanobis','neighbor'),times=15),sep='_')
plotSDRatios(temp, 'Comparing SD Ratios of All Propensity Score Models\n and Matching Methods Used')

## Using var as id variables

---

Question 6: Repeat steps 2-4, but this time using IPTW.

ps.IPTW <- glm(treatment ~ ., data=dt[, c('treatment',confounders)], family=binomial(link='logit'))
dt$psc_iptw <- ps.IPTW$fitted.values

For estimating the ATT, define the weight as follows:

\[\omega(Z,x) = Z + (1-Z) * \frac{\hat{\epsilon}(x)}{1-\hat{\epsilon}(x)}\]

dt$r_att = dt$psc + (1-dt$treat) * dt$psc / (1-dt$psc)

iptw.balance <- checkBalance(dt, confounders, 'r_att')

---

Question 7: Comparative balance table. Create a table with columns 6 and 8 from your function for each of the matching and weighting methods performed above. Which approach would you choose and why? (1-2 paragraphs at most)

compareBalance <- function(model_bals, model_names, col) {
    # model_bals: needs to be a list() of outputs from checkBalance
    # model_names: vector of model names in same order as bals
    # col: col to extract for output: [1,8]

    # number of models
    n_models <- length(model_bals)
    # subset out columns 6 & 8
    temp.ls <- lapply(model_bals, function(bal) bal[,col])
    # cbind list of df's into one df
    temp.df <- data.frame(do.call('cbind', temp.ls))
    # kable output
    names(temp.df) <- model_names
    if (interactive()) {
        return( round(temp.df, 2) )
    } else {
        out <- knitr::kable(temp.df, digits=2, align='c', format='latex')
        if (col==6) {
            out <- out %>% add_header_above(c(' '=1, 'Mean Difference (matched)'=n_models))
        } else {
            out <- out %>% add_header_above(c(' '=1, 'Ratio sd (matched)'=n_models))
        }
        return( out )
    }
}

model_bals <- list(logit.bal, probit.bal, logit_nr.bal, mh.bal, mk2.bal, iptw.balance)
model_names <- c('Logistic', 'Probit', 'Logistic (no replacement)', 'Mahalanobis Matching', 'K=2 Nearest Neighbors', 'IPTW')

compareBalance(model_bals, model_names, 6)

compareRatios <- compareBalance(model_bals, model_names, 8)

我们主要关注期望的差异

---

Question 8: Estimate the treatment effect for the restructured datasets implied by Questions 4-6 (Step 5)

Estimate the effect of the treatment on the treated for each of your five approaches by fitting a regression with weights equal to the number of times each observation appears in the matched sample (that is, use your weights variable from above) or using IPTW weights. Report the treatment effect and standard error for each approach.

temp <- dt[,c('outcome', 'treatment', confounders)]

mod1 = summary(lm(outcome ~ ., data=temp, weights=dt$wt))$coefficients
mod1.out = c(mod1['treatment', 'Estimate'], mod1['treatment', 'Std. Error'])

mod2 = summary(lm(outcome ~ ., data=temp, weights=dt$w_probit))$coefficients
mod2.out = c(mod2['treatment', 'Estimate'], mod2['treatment', 'Std. Error'])

mod3 = summary(lm(outcome ~ ., data=temp, weights=dt$w_logit_no_rep))$coefficients
mod3.out = c(mod3['treatment', 'Estimate'], mod3['treatment', 'Std. Error'])

mod4 = summary(lm(outcome ~ ., data=temp, weights=dt$wt_mh))$coefficients
mod4.out = c(mod4['treatment', 'Estimate'], mod4['treatment', 'Std. Error'])

mod5 = summary(lm(outcome ~ ., data=temp, weights=dt$wt_k2))$coefficients
mod5.out = c(mod5['treatment', 'Estimate'], mod5['treatment', 'Std. Error'])

mod6 = summary(lm(outcome ~ ., data=temp, weights=dt$r_att))$coefficients
mod6.out = c(mod5['treatment', 'Estimate'], mod5['treatment', 'Std. Error'])

q8_results = rbind(mod1.out, mod2.out, mod3.out, mod4.out, mod5.out, mod6.out)
row.names(q8_results) = c('logit (w/ replacement)', 'probit (w/ replacement)', 'logit (no replacement)', 'Mahalanobis Matching', 'K-2 Nearest Neighbors', 'IPTW')
colnames(q8_results) = c('Treatment Effect','Std. Error')
kable(q8_results, digits=2, align='c', format='markdown', caption='Treatment Effect Estimation by Propensity Score Model')

Treatment Effect Estimation by Propensity Score Model

	Treatment Effect	Std. Error
logit (w/ replacement)	3.95	0.15
probit (w/ replacement)	3.86	0.14
logit (no replacement)	3.95	0.15
Mahalanobis Matching	3.98	0.15
K-2 Nearest Neighbors	3.85	0.14
IPTW	3.85	0.14

---

Question 9: Assumptions What assumptions are necessary to interpret the estimates from the propensity score approaches causally? List and describe briefly.

(a) Ignorability

Ignorability holds if the covariates in the propensity score model are the only confounding covariates and we match on the propensity score. So we need to assume that we have controlled for all the potential confounders. Then we can unbiasedly estimate \(E[Y(1)| Z=1, e(X)]\) with \(\bar{Y}_{Z=1,e(X)}\) and \(E[Y(0)| Z=1, e(X)]\) with \(\bar{Y}_{Z=0,e(X)}\).

(b) Sufficient Overlap

If we are interested in the effect of the treatment on the treated, we want to make sure that for each treatment group member there is a control group member that is sufficiently similar. Then, we can use this control group member as an empirical counterfactual.

Similarly, if we are interested in the effect of the treatment on the control, we want to make sure that for each control group member there is a treatment group member that is sufficiently similar. Then, we can use this treatment group member as an empirical counterfactual.

If we are interested in the average treatment effect, we want to make sure that for each subject in either treatment or control, there is an empirical counterfactual (i.e. a sufficiently similar control subject for a treatment subject and vise versa).

(c) Appropriate specification of the propensity score model / balance achieved

An appropriate specification of the propensity score model ensures that we have good overlap and balance. Overlap is important with regard to drawing causal inference due to reasons listed above. Balance is important because imbalance forces us to rely more on the correct functional form of our model. Incorrect functional forms would lead to biased estimates of the treatment effect.

(d) SUTVA

Stable unit treatment value assumption states that the treatment effect does not depend on the particular configuration of treatment assignment. That is, there is no dilution or concentration of the treatment effect.

---

Question 10: Causal Interpretation Provide a causal interpretation of one of your estimates above. Remember to specify the counterfactual and to be clear about whom you are making inferences. Also make sure to use causal (counterfactual) language.

For children who participated in the IDHP intervention, their IQ scores at age 3 were 4 points higher than had they not participated in the intervention.

---

Question 11: Comparison to linear regression Fit a regression of your outcomes to the treatment indicator and covariates.

(a) Report your estimate and standard error.

lm1 <- lm(outcome ~ ., data=dt[, c(confounders, 'outcome', 'treatment')])
round(summary(lm1)$coefficients['treatment', c('Estimate', 'Std. Error'), drop=FALSE], 2)

##           Estimate Std. Error
## treatment     3.88       0.12

(b) Interpret your results non-causally

We expect that a group of children who participate in the IHDP intervention will have, on average, an 11.54 point increase in their age 3 IQ scores compared to a group of children who did not participate in the intervention, controlling for all other covariates.

(c) Why might we prefer the results from the propensity score approach to the linear regression results in terms of identifying a causal effect?

In terms of estimating causal effects, propensity score approaches are preferred over linear regression to account for the systematic differences between the treatment and control groups in key covariates that can affect the outcome. While linear regression can account for some group differences, the method is highly dependent on correct model specification in order avoid treatment effect estimation bias. In comparison, propensity score approaches are non-parametric so we avoid strong assumptions about correct model specification. Additionally, propensity score approaches decrease this bias by restructuring the sample of controls and/or treated observations so that there is overlap in the common support of the covariates between groups, and improved balance in distributions of the covariates between groups. Propensity score approaches decrease the treatment effect estimation bias caused by extrapolating into areas of the covariate space outside of the support of the data, and incorrect model specification of the treatment and covariates on the outcome.

However, propensity score approaches still rely on the strong assumption that there is no omitted variable bias; that is, all confounders have been controlled for.