This homework covers Modules 2 and 3.
Students may discuss homework exercises with one another. However, sharing, copying, or providing any part of a completed homework solution or code is an infraction of the University’s rules on Academic Integrity. Any violation will be punished as severely as possible. Final submissions must be uploaded to Gradescope using entry code XE3YJ6 (course ID 1370811). No email or hard copy will be accepted. For late submission policy and grading rubrics, please refer to the course website.
HW02_yourNetID.pdf. For example,
HW02_rqzhu.pdf. Please note that this must be a
.pdf file. .html format
cannot be accepted since it may not be readable in
Gradescope. All proofs must be typed in LaTeX format. Make all of your
code visible for grading.R or Python. The
supplied code uses R; equivalent Python
implementations are acceptable..Rmd file
as a template, be sure to remove this instruction
section.We learned a few estimators of the average treatment effect (ATE) in the lecture. In this question, we will compare their performance in a simulation study. Use the following data generating process. We simulate i.i.d. samples \((X,A,Y)\) with nonlinear effects. Let
\[ X=(X_1,X_2,X_3,X_4,X_5), \]
where the covariates are mutually independent, with \(X_1,X_2\sim\mathcal{N}(0,1)\), \(X_3\sim\mathrm{Unif}[-1,1]\), \(X_4\sim\mathrm{Bernoulli}(0.4)\), and \(X_5\sim\mathrm{Exp}(1)\). Generate \(A\mid X\sim\mathrm{Bernoulli}(e(X))\), with
\[ \operatorname{logit}\{e(X)\} = 0.6X_1+0.6X_2^2-0.8\mathbf{1}\{X_3>0\}. \]
The potential outcomes are generated with
\[ Y(0) = 1+0.8X_1+0.5X_4+0.3\log(1+X_5)+\varepsilon, \]
where \(\varepsilon\sim\mathcal{N}(0,1)\) is independent of \((X,A)\). Let \(Y(1)=Y(0)+\tau(X)\), with
\[ \tau(X) = 1-0.3X_2+0.4\mathbf{1}\{X_3>0\}+0.2\sqrt{X_5}. \]
\[ \tau=\mathbb{E}[\tau(X)] \]
using a Monte Carlo approach with 10,000 independent samples.
\[ Y=AY(1)+(1-A)Y(0) \]
instead of the potential outcomes. Simulate a dataset of size \(n=500\) from the above data generating process. Implement the IPW estimator
\[ \widehat{\tau}_{\text{IPW}} = \frac{1}{n} \sum_{i=1}^n \left( \frac{A_iY_i}{\widehat{e}(X_i)} - \frac{(1-A_i)Y_i}{1-\widehat{e}(X_i)} \right), \]
where \(\widehat{e}(X)\) is estimated using a logistic regression model with the main effects of \(X_1,\ldots,X_5\).
Repeat this simulation independently 500 times and report the mean, standard deviation, and bias of the IPW estimator. Compare your result with the ATE approximated in part (a) and briefly discuss your findings.
Randomly split each simulated dataset into two folds. For observations in each fold, fit both the propensity score and outcome models using only the other fold. Use these predictions in the two-fold cross-fitted doubly robust estimator
\[ \widehat{\tau}_{\text{DR}} = \frac{1}{n} \sum_{i=1}^n \left[ \widehat{m}_1(X_i)-\widehat{m}_0(X_i) + \frac{A_i\{Y_i-\widehat{m}_1(X_i)\}} {\widehat{e}(X_i)} - \frac{(1-A_i)\{Y_i-\widehat{m}_0(X_i)\}} {1-\widehat{e}(X_i)} \right]. \]
Report the mean, standard deviation, and bias of the DR estimator over the same 500 datasets used in part (b). Compare its performance with the IPW estimator and briefly explain your result, considering whether the conditions for double robustness hold in this setting.
Now consider the Hájek estimator, where \(\widehat e_i=\widehat e(X_i)\):
\[ \widehat{\tau}_{\text{Hajek}} = \frac{ \sum_{i=1}^n \frac{A_iY_i}{\widehat{e}_i} }{ \sum_{i=1}^n \frac{A_i}{\widehat{e}_i} } - \frac{ \sum_{i=1}^n \frac{(1-A_i)Y_i}{1-\widehat{e}_i} }{ \sum_{i=1}^n \frac{1-A_i}{1-\widehat{e}_i} }. \]
Show that the Hájek estimator is location invariant.
In practice, it may be difficult to know whether the propensity score model is correctly specified, and we may also have unobserved confounders. Consider the following simple regression model
\[ Y=\beta_0+X_1+A\cdot X_2+\varepsilon, \]
where \(X_1\) and \(X_2\) are jointly normally distributed with mean 0 but unknown covariance structure, and
\[ \varepsilon\sim\mathcal{N}(0,1). \]
Here \(A\in\{0,1\}\), and \(\varepsilon\) is independent of \((X_1,X_2,A)\). You may choose \(\beta_0\), but report its value.
In this question, suppose the true propensity score depends on both \(X_1\) and \(X_2\), but you only observe \(X_1\). When estimating the propensity score, you will therefore use a logistic regression model based only on \(X_1\).
You are asked to construct two models for the true propensity score
\[ P(A=1\mid X_1,X_2) \]
and the covariance structure of \(X_1\) and \(X_2\) such that
the IPW estimator for the ATE is systematically biased upward;
the IPW estimator for the ATE is systematically biased downward.
Explain your choices and why you expect the corresponding direction of bias. Use the IPW estimator defined in Question 1(b).
Demonstrate your idea using a simulation study with sample size \(n=200\) and 100 replications. Report the mean and standard deviation of the IPW estimator in both cases, and compare them with the true ATE, which can be obtained either analytically or using simulation.
The average treatment effect summarizes the treatment effect over the whole population. In many applications, however, the treatment effect may vary across individuals. The conditional average treatment effect (CATE) is defined as
\[ \tau(x) = \mathbb{E}[Y(1)-Y(0)\mid X=x]. \]
Consider the following data generating process. Let the covariates be mutually independent, with
\[ X_1,X_2,X_3\sim\mathcal{N}(0,1), \]
and generate \(A\mid X\sim\mathrm{Bernoulli}(e(X))\), where
\[ \operatorname{logit}\{e(X)\} = 0.5X_1-0.5X_2. \]
The potential outcomes are generated according to
\[ Y(0) = X_1+X_2^2+\varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,1), \]
and
\[ Y(1)=Y(0)+\tau(X), \]
where
\[ \tau(X)=1.5X_1-X_2. \]
The noise \(\varepsilon\) is independent of \((X,A)\), and the observed outcome is \(Y=AY(1)+(1-A)Y(0)\).
[5 pts] Write down the true CATE \(\tau(x)\). Assuming a larger outcome is preferred, what is the optimal individualized treatment rule?
[7 pts] Generate a training dataset with \(n=1000\) and an independent testing dataset of size 5000 for parts (b) and (c). Implement either the T-learner (separate outcome regressions for the two treatment groups) or the X-learner to estimate the CATE using the training data. You may use a regression method of your choice for the regression steps, but you should briefly justify your choice.
Based on your estimated CATE, define the individualized treatment rule
\[ \widehat{d}(X) = \mathbf{1}\{\widehat{\tau}(X)>0\}. \]
On the testing data, compare the estimated treatment rule with the true optimal treatment rule and report the proportion of observations for which the two decisions agree.
treat everyone;
treat no one;
use your estimated individualized treatment rule.
Since the potential outcomes are known in this simulation, calculate the average outcome under each treatment rule. Which rule performs the best?
Instead of first estimating the CATE and then constructing a treatment rule, outcome weighted learning directly estimates an individualized treatment rule by reformulating the problem as a weighted classification problem.
Suppose \(A\in\{-1,1\}\) is a binary treatment, \(X\) is a vector of covariates, and \(R\) is an outcome, where a larger value of \(R\) is preferred. Let \(d(X)\in\{-1,1\}\) denote a treatment rule.
\[ V(d) = \mathbb{E} \left[ \frac{ R\mathbf{1}\{A=d(X)\} }{ P(A\mid X) } \right]. \]
Explain briefly why maximizing the empirical version of this quantity can be viewed as a weighted classification problem when the outcome weights are nonnegative.
For the remaining parts, use the following simulation model.
set.seed(1)
n <- 800
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
A <- 2 * rbinom(n, 1, 0.5) - 1
side <- sign(x2 - sin(2 * pi * x1) / 3 - 0.5)
R <- rnorm(
n,
mean = ifelse(side == A, 1.5, 0.5),
sd = 0.5
)[3 pts] Based on the data generating mechanism, determine the true optimal treatment rule as a function of \(x_1\) and \(x_2\). Plot the observations and add the true decision boundary.
[4 pts] Estimate \(e(X)=P(A=1\mid X)\) using a logistic regression on \(x_1\) and \(x_2\). For binomial regression, code the response as \((A+1)/2\) or as a factor. Construct the outcome weights
\[ W_i = \frac{R_i+c} {\widehat{P}(A_i\mid X_i)}, \]
where \(\widehat P(A_i\mid X_i)\) is \(\widehat e(X_i)\) if \(A_i=1\), and \(1-\widehat e(X_i)\) if \(A_i=-1\). Choose and report a constant \(c\) that makes all weights nonnegative.
Use these weights to fit a weighted classifier for the treatment label \(A\). You may use a weighted logistic regression or another classification method discussed in class. Use the fitted classifier to obtain an estimated treatment rule.
A fitted random forest can be viewed as defining a similarity measure between observations. For a forest with \(B\) trees, define
\[ K(x,z) = \frac{1}{B} \sum_{b=1}^B \mathbf{1}\{x \text{ and } z \text{ belong to the same terminal node in tree } b\}. \]
\[ K_b(x,z) = \mathbf{1}\{x \text{ and } z \text{ belong to the same terminal node}\} \]
is a positive definite kernel.
\[ K(x,z) = \frac{1}{B}\sum_{b=1}^B K_b(x,z) \]
is also positive definite.
Briefly explain what a large value of \(K(x,z)\) means in terms of the similarity between observations \(x\) and \(z\).
A Mondrian tree generates a random recursive partition of the input space. Unlike a standard decision tree, the partition does not depend on the observed responses. In this question, you will implement a simple Mondrian tree and use it to construct a Mondrian forest.
\(\lambda\), the budget parameter controlling the complexity of the partition;
\(p\), the dimension of the input space.
Your function should return the fitted tree structure. You may store the tree in a matrix with the following columns:
node type: 1 for an internal node and
-1 for a leaf node;
splitting dimension;
cut location;
index of the left child;
index of the right child;
prediction value, which can be left as NA for
now.
You should also keep track of the lower and upper bounds of the bounding box associated with each node. To give you an idea, here is an example
node_type split_dim cut_loc left_child right_child pred_value
[1,] 1 2 0.1303658 2 3 NA
[2,] 1 1 0.2542666 4 5 NA
[3,] 1 1 0.1783499 6 7 NA
[4,] -1 NA NA NA NA NA
[5,] -1 NA NA NA NA NA## [,1] [,2]
## [1,] 0.0000000 0.0000000
## [2,] 0.0000000 0.0000000
## [3,] 0.0000000 0.1303658
## [4,] 0.0000000 0.0000000
## [5,] 0.2542666 0.0000000## [,1] [,2]
## [1,] 1.0000000 1.0000000
## [2,] 1.0000000 0.1303658
## [3,] 1.0000000 1.0000000
## [4,] 0.2542666 0.1303658
## [5,] 1.0000000 0.1303658Your code should contain a while or for
loop that follows the Mondrian splitting mechanism in the lecture.
Stop without splitting if the next event time would exceed the lifetime
\(\lambda\). You may use
rexp() in R or its equivalent in
Python to generate exponential random variables.
Run your function with \(p=2\) and \(\lambda=2\). Print up to the first five rows of the tree structure and the corresponding bounding boxes.
rect() function or an equivalent plotting function to draw
the bounding boxes of all terminal nodes.Generate 300 observations uniformly on \([0,1]^2\). Write a function that determines the terminal node containing each observation, and visualize the observations using a different color for each terminal node.
\[ Y = \sin(2\pi X_1)+2X_2^2+\varepsilon, \qquad \varepsilon\sim\mathcal{N}(0,1), \]
where \(X\) is uniformly distributed on \([0,1]^2\) and \(\varepsilon\) is independent of \(X\).
Generate 400 training observations and an independent sample of 1000 testing observations.
For each terminal node, calculate the average response of all training observations falling into that node and use this value as the prediction for the node. If a terminal node contains no training observations, use the mean response of the training sample used to fit that tree.
Fit a Mondrian tree with \(\lambda=4\) and report the mean squared error on the testing data. You may also experiment with other values of \(\lambda\).
Use the same training and testing data as in part (c). Report the testing MSE for the Mondrian forest. Experiment with different values of \(\lambda\), comparing the forest and single tree at the same values of \(\lambda\).
Is the Mondrian forest able to improve prediction performance? Briefly discuss your findings.