\documentclass[notitlepage]{article} \usepackage{../../Math598} \usepackage{listings} \usepackage{numprint} \usepackage{bbm} \usepackage{amsfonts,amsmath} \usepackage{fancyhdr} \usepackage{lastpage} \usepackage{enumitem} \usepackage{bbold} \usepackage{hyperref} \hypersetup{ colorlinks=true, linkcolor=blue, filecolor=magenta, urlcolor=blue, pdftitle={P3}, bookmarks=false } \usepackage{tikz} \usetikzlibrary{shapes,decorations,arrows,calc,arrows.meta,fit,positioning} \def\E{\Expect} \lstloadlanguages{R} \definecolor{keywordcolor}{rgb}{0,0.6,0.6} \definecolor{delimcolor}{rgb}{0.461,0.039,0.102} \definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432} \lstdefinestyle{Rsettings}{ basicstyle=\ttfamily, breaklines=true, showstringspaces=false, keywords={if, else, function, theFunction, tmp}, % Write as many keywords otherkeywords={}, commentstyle=\itshape\color{Rcommentcolor}, keywordstyle=\color{keywordcolor}, moredelim=[s][\color{delimcolor}]{"}{"}, } \lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}} \lfoot{MATH 598/782 Project \ExNumber~Solutions} \cfoot{ } \rfoot{{\small{\em Page \thepage \ of \pageref{LastPage}}}} \pagestyle{empty} \def\ExNumber{3} \parindent 0in \renewcommand{\headrulewidth}{0.0pt} \renewcommand{\footrulewidth}{0.0pt} \def\bU{\mathbf{U}} \def\bR{\mathbf{R}} \begin{document} <>= library(knitr) # global chunk options opts_chunk$set(cache=TRUE, autodep=TRUE) options(scipen=6) options(repos=c(CRAN="https://cloud.r-project.org/")) knit_hooks$set(document = function(x) {sub('\\usepackage[]{color}', '\\usepackage{xcolor}', x, fixed = TRUE)}) @ \pagestyle{fancy} \begin{center} {\textsclarge{MATH 598/782 -- Project \ExNumber~Solutions }} \end{center} \bigskip \tikzset{ -Latex,auto,node distance =1 cm and 1 cm,semithick, state/.style ={circle, draw, minimum width = 0.7 cm}, box/.style ={rectangle, draw, minimum width = 0.7 cm, fill=lightgray}, point/.style = {circle, draw, inner sep=0.08cm,fill,node contents={}}, bidirected/.style={Latex-Latex,dashed}, el/.style = {inner sep=3pt, align=left, sloped} } \textit{For the $p \times 1$ system of estimating equations \[ \sum_{i=1}^n \bU_i(\theta) = \zerovec \] with $\E[\bU(\theta_0)] = \zerovec$, we have that \[ \sqrt{n}(\widehat \theta_n - \theta_0) \CiD Normal_p(\zerovec,\bV) \] where \[ \bV = \mathcal{J}^{-1} \mathcal{I} \mathcal{J}^{-\top} \] with \[ \mathcal{I} = \E[ \bU(\theta_0) \bU(\theta_0)^\top] \qquad \qquad \mathcal{J} = \E[\dot{\bU}(\theta_0)] \] both $(p \times p)$ matrices, and \[ \dot{\bU}(\theta_0) = -\dfrac{\partial \bU(\theta)}{\partial \theta^\top} \bigg |_{\theta = \theta_0} \] } \textit{We proceed by computing estimates, $\widehat I_n$ and $\widehat J_n$, of the two matrices based on the observed data, and then computing \[ \widehat \bV = \widehat J_n^{-1} \widehat I_n \widehat J_n^{-\top} \] to estimate the asymptotic variance; this approach is known as \textit{sandwich} or \textit{robust} variance estimation. } \bigskip \textit{Based on the data set from Project 2 \begin{center} \url{https://www.math.mcgill.ca/dstephens/598-Causal-2024/Exercises/Project2/P2data.csv} \end{center} which comprises four predictor variables $X_1,X_2,X_3,X_4$, a binary treatment variable $Z$, and outcome variable $Y$, compute an estimate of the asymptotic variance of the parameter of interest $\psi$ based on the sandwich estimation procedure and knitr sheets for the following models and estimating procedures.} \begin{enumerate}[label=(\alph*)] \item \textit{An outcome regression model \[ Y_i = \beta_0 + \beta_1 X_{1i} + \psi Z_i + \epsilon_i \] estimated using ordinary least squares.} \medskip \textbf{Solution:} In this case, we can use the standard results for OLS. There are no additional parameters to estimate, and we have that \[ \bU(\theta) = \bx^\top (Y - \bx \theta) \qquad \dot{\bU}(\theta_0) = \bx^\top \bx \] where $\theta = (\beta_0,\beta_1,\psi)^\top$ and $\bx = (1,X_1,Z)$. Therefore the asymptotic variance is \[ V = \{\E[\bx^\top \bx]\}^{-1}\E[(Y-\bx \theta_0)^2 \bx^\top \bx] \{\E[\bx^\top \bx]\}^{-1} \] that we can estimate using the quantities \[ \widehat I_n = \frac{1}{n} \sum_{i=1}^n (y_i - \bx_i \widehat \theta_n)^2 \bx_i^\top \bx_i \qquad \qquad \widehat J_n = \frac{1}{n} \sum_{i=1}^n \bx_i^\top \bx_i \] before extracting the component for $\psi$. <>= fn<-"https://www.math.mcgill.ca/dstephens/598-Causal-2024/Exercises/Project2/P2data.csv" P3.data<-read.csv(fn) n<-nrow(P3.data) Xm<-cbind(1,P3.data$X1,P3.data$Z) Y<-P3.data$Y Z<-P3.data$Z X1<-P3.data$X1 X2<-P3.data$X2 theta.hat<-solve(t(Xm) %*% Xm) %*% (t(Xm) %*% Y) theta.hat All.ests<-All.se<-rep(0,6) All.ests[1]<-theta.hat[3] eps.y<-Y - Xm %*% theta.hat I.n<-crossprod(Xm*as.vector(eps.y))/n; round(I.n,5) J.n<-crossprod(Xm)/n; round(J.n,5) #Asymptotic variance V.n<-solve(J.n) %*% (I.n %*% solve(J.n));round(V.n,5) #Approximate estimate variance round(V.n/n,5) se.vec<-sqrt(diag(V.n)/n); round(se.vec,5) All.se[1]<-se.vec[3] @ Therefore the approximate standard errors for the three estimates are the square roots of the diagonal of this matrix divided by $n$, that is \begin{table}[ht] \centering \begin{tabular}{|c|c|c|} \hline & Est. & s.e. \\[3pt] \hline & & \\[-6pt] $\widehat \beta_0$ & \Sexpr{round(theta.hat[1],4)} & \Sexpr{round(se.vec[1],4)} \\[3pt] $\widehat \beta_1$ & \Sexpr{round(theta.hat[2],4)} & \Sexpr{round(se.vec[2],4)} \\[3pt] $\widehat \psi$ & \Sexpr{round(theta.hat[3],4)} & \Sexpr{round(se.vec[3],4)} \\[3pt] \hline \end{tabular} \end{table} \hfill{4 Marks} For this problem, we can compare with the ordinary \texttt{lm} results (that assume correct specification, and estimates the homoscedastic error variance using the standard approach) <>= fit.a<-lm(Y~X1+Z,data=P3.data) round(summary(fit.a)$coef,4) se.standard<-summary(fit.a)$coef[,2] @ These standard errors are similar to, but slightly different from, the sandwich versions. \begin{table}[ht] \centering \begin{tabular}{|c|c|c|} \hline & Sandwich & \texttt{lm}\\[3pt] \hline & & \\[-6pt] $\widehat \beta_0$ & \Sexpr{round(se.vec[1],4)} & \Sexpr{round(se.standard[1],4)}\\[3pt] $\widehat \beta_1$ & \Sexpr{round(se.vec[2],4)} & \Sexpr{round(se.standard[2],4)} \\[3pt] $\widehat \psi$ & \Sexpr{round(se.vec[3],4)} & \Sexpr{round(se.standard[3],4)}\\[3pt] \hline \end{tabular} \end{table} We can also use the \texttt{sandwich} package to obtain the sandwich variance estimates. <>= library(sandwich) sandwich(fit.a) @ This matrix is the $\widehat{V}_n/n$ matrix, and matches the above calculation. \medskip \item \textit{A propensity score regression model \[ Y_i = \beta_0 + \psi Z_i + \phi e(X_{1i},X_{2i};\widehat \alpha) + \epsilon_i \] where $\alpha$ is estimated using maximum likelihood in a logistic regression of $Z$ on $(X_1,X_2)$, including an intercept, that is $Z|X_1=x_1,X_2 = x_2 \sim Bernoulli(e(x_1,x_2;\alpha))$ with $\alpha = (\alpha_0,\alpha_1,\alpha_2)^\top$, and \[ e(x_1,x_2;\alpha) = \frac{\exp\{\alpha_0 + \alpha_1 x_1 + \alpha_2 x_2 \}}{1+\exp\{\alpha_0 + \alpha_1 x_1 + \alpha_2 x_2 \}}. \]} \bigskip \textbf{Solution:} For this propensity score regression model with estimated propensity score, we have \[ \bU(\theta) = \begin{pmatrix} (Y - \beta_0 - \psi Z - \phi e(X;\alpha))\\[12pt] Z (Y - \beta_0 - \psi Z - \phi e(X;\alpha))\\[12pt] e(X;\alpha) (Y - \beta_0 - \psi Z - \phi e(X;\alpha))\\[12pt] \bx^\top (Z - e(X; \alpha)) \end{pmatrix} = \begin{pmatrix} \varepsilon_Y \\[6pt] Z \varepsilon_Y \\[6pt] e(X;\alpha) \varepsilon_Y \\[6pt] \bx^\top \epsilon_Z \end{pmatrix}. \] where $\bx=(1,x_1,x_2)$. We have therefore that \[ \mathcal{I} = \E \begin{bmatrix} \varepsilon_Y^2 & Z \varepsilon_Y^2 & e \varepsilon_Y^2 & \bx \epsilon_Z \varepsilon_Y \\[6pt] Z \varepsilon_Y^2 & Z^2 \varepsilon_Y^2 & e Z \varepsilon_Y^2 & \bx Z \epsilon_Z \varepsilon_Y \\[6pt] e \varepsilon_Y^2 & e Z \varepsilon_Y^2 & e^2 \varepsilon_Y^2 & \bx e \epsilon_Z \varepsilon_Y\\[6pt] \bx^\top \epsilon_Z \varepsilon_Y & \bx^\top Z \epsilon_Z \varepsilon_Y & \bx^\top e \epsilon_Z \varepsilon_Y & \bx^\top \bx \epsilon_Z^2 \end{bmatrix} \qquad \qquad \mathcal{J} = \E \begin{bmatrix} 1 & Z & e & \phi \bx e(1-e)\\[6pt] Z & Z^2 & Z e & \phi Z \bx e(1-e) \\[6pt] e & Z e & e^2 & (\phi e - \varepsilon_Y) \bx e(1-e)\\[6pt] \zerovec & \zerovec & \zerovec & \bx^\top \bx e (1-e) \end{bmatrix}. \] This follows in the third line as by the chain rule \begin{align*} \frac{\partial}{\partial \alpha} \left\{e(X;\alpha) (Y - \beta_0 - \psi Z - \phi e(X;\alpha)) \right\} & = - \phi e(X;\alpha) \frac{\partial e(X;\alpha)}{\partial \alpha} + \frac{\partial e(X;\alpha)}{\partial \alpha} (Y - \beta_0 - \psi Z - \phi e(X;\alpha))\\[6pt] & = \frac{\partial e(X;\alpha)}{\partial \alpha}( - \phi e(X;\alpha) + \varepsilon_Y)\\[6pt] & = e(X;\alpha) (1-e(X;\alpha)) \bx^\top (\varepsilon_Y - \phi e(X;\alpha)). \end{align*} Note that \[ \frac{\partial e(X;\alpha)}{\partial \alpha} = e(X;\alpha) (1-e(X;\alpha)) \bx^\top \] as under the logistic regression propensity model \[ \frac{\partial e(x)}{\partial x} = \frac{\partial}{\partial x} \left\{\frac{e^x}{1+e^x} \right\} = \frac{e^x}{(1+e^x)^2} = \frac{e^x}{1+e^x} \frac{1}{1+e^x} = \frac{e^x}{1+e^x} \left(1-\frac{e^x}{1+e^x} \right) = e(x) (1-e(x)) \] <>= eX<-fitted(glm(Z~X1+X2,family=binomial,data=P3.data)) vx<-eX*(1-eX) Ptmp<-P3.data Ptmp$eX<-eX fit.psr<-lm(Y~Z+eX,data=Ptmp) psr.ests<-coef(fit.psr) All.ests[2]<-psr.ests[2] Xa<-cbind(1,X1,X2) eps.y<-residuals(fit.psr) eps.z<-Z-eX I.n<-crossprod(cbind(eps.y,Z*eps.y,eX*eps.y,Xa*eps.z))/n; round(I.n,5);I.n.psr<-I.n J.n<-matrix(0,6,6) Xp<-cbind(1,Z,eX) J.n[1:3,1:3]<-crossprod(Xp)/n J.n[4:6,4:6]<-crossprod(Xa*sqrt(vx))/n J.n[1,4]<-mean(psr.ests[3]*vx) J.n[1,5]<-mean(psr.ests[3]*X1*vx) J.n[1,6]<-mean(psr.ests[3]*X2*vx) J.n[2,4]<-mean(psr.ests[3]*Z*vx) J.n[2,5]<-mean(psr.ests[3]*Z*X1*vx) J.n[2,6]<-mean(psr.ests[3]*Z*X2*vx) J.n[3,4]<-mean((psr.ests[3]*eX-eps.y)*vx) J.n[3,5]<-mean((psr.ests[3]*eX-eps.y)*X1*vx) J.n[3,6]<-mean((psr.ests[3]*eX-eps.y)*X2*vx) J.n.psr<-J.n round(J.n,5) #Asymptotic variance V.n<-solve(J.n) %*% (I.n %*% t(solve(J.n))); round(V.n,5); V.n.psr<-V.n #Approximate estimate variance round(V.n/n,5) se.vec<-sqrt(diag(V.n)/n); se.vec All.se[2]<-se.vec[2] @ Thus the $\psi$ parameter has estimate \Sexpr{round(psr.ests[2],5)} with standard error \Sexpr{round(se.vec[2],5)}. We can use the \texttt{sandwich} package to compute the sandwich variance in the outcome model that assumes that the plug-in propensity model is not in fact estimated. <>= library(sandwich) sandwich::sandwich(fit.psr) @ Therefore this estimated standard error is $\sqrt{\Sexpr{round(sandwich(fit.psr)[2,2],5)}}$, that is \Sexpr{round(sqrt(sandwich(fit.psr)[2,2]),5)}. However, this variance does not accurately reflect the sampling variability. \hfill{6 Marks} \pagebreak \item \textit{An inverse probability weighting (IPW) procedure using the propensity score model from (b).} \medskip \textbf{Solution:} For this, we can simply deploy the code from knitr sheet 12, and either of the IPW estimators. In this case, we have to estimate the parameters using the logistic regression model. \begin{itemize} \item \textbf{IPW version 1:} For the elementary application of IPW, there is a single unknown parameter of interest, and the estimator is \[ \widehat \psi = \frac{1}{n} \sum_{i=1}^n \left(\frac{Z_i}{e(X)} - \frac{1-Z_i}{1-e(X_i)} \right)Y_i \] For the unknown propensity score case, the score is \[ \bU(\theta) = \begin{pmatrix} \left(\dfrac{Z}{e(X;\alpha)} - \dfrac{1-Z}{1-e(X; \alpha)} \right)Y - \psi \\[12pt] \bx^\top (Z - e(X; \alpha)) \end{pmatrix} = \begin{pmatrix} \zeta_Y \\[6pt] \bx^\top \epsilon_Z \end{pmatrix}. \] say. Therefore \[ \mathcal{I} = \E \begin{bmatrix} \zeta_Y^2 & \bx \epsilon_Z \zeta_Y \\[6pt] \bx^\top \epsilon_Z \zeta_Y & \bx^\top \bx \epsilon_Z^2 \end{bmatrix} \qquad \qquad \mathcal{J} = \E \begin{bmatrix} 1 & \left(\dfrac{Z}{e^2} + \dfrac{1-Z}{(1-e)^2} \right) Y \bx e(1-e) \\[12pt] \zerovec & \bx^\top \bx e (1-e) \end{bmatrix} \] as \[ \frac{\partial }{\partial \alpha} \left\{ \left(\dfrac{Z}{e(X;\alpha)} - \dfrac{1-Z}{1-e(X; \alpha)} \right)Y - \psi \right\} = \left( - \dfrac{Z}{e(X;\alpha)^2} - \dfrac{1-Z}{(1-e(X; \alpha))^2} \right)Y \frac{\partial e(X;\alpha) }{\partial \alpha} \] with the final derivative as given above. <>= eX<-fitted(glm(Z~X1+X2,family=binomial,data=P3.data)) wv<-Z/eX-(1-Z)/(1-eX) psi.hat.ipw1<-mean(wv*Y) psi.hat.ipw1 All.ests[3]<-psi.hat.ipw1 Xa<-cbind(1,X1,X2) eps.y<-wv*Y-psi.hat.ipw1 eps.z<-Z-eX I.n<-crossprod(cbind(eps.y,Xa*eps.z))/n; round(I.n,5) vx<-eX*(1-eX) J.n<-matrix(0,4,4) J.n[1,1]<-1 J.n[1,2]<-mean(vx*Y*(Z/eX^2+(1-Z)/(1-eX)^2)) J.n[1,3]<-mean(X1*vx*Y*(Z/eX^2+(1-Z)/(1-eX)^2)) J.n[1,4]<-mean(X2*vx*Y*(Z/eX^2+(1-Z)/(1-eX)^2)) J.n[2:4,2:4]<-crossprod(Xa*sqrt(vx))/n round(J.n,5) #Asymptotic variance V.n<-solve(J.n) %*% (I.n %*% t(solve(J.n))); round(V.n,5) #Approximate estimate variance round(V.n/n,5) se.vec<-sqrt(diag(V.n)/n); se.vec All.se[3]<-se.vec[1] @ Therefore the estimate is \Sexpr{round(psi.hat.ipw1,5)} with standard error \Sexpr{round(se.vec[1],5)}. \medskip \item \textbf{IPW version 2:} For the second version of the IPW estimator, we consider weighted least squares and the objective function \[ \sum_{i=1}^n R_i (Y_i - \beta_0 - Z_i \psi)^2 \] where \[ R_i = \left(\dfrac{1-Z_i}{1-e(X_i)}+ \dfrac{Z_i}{e(X_i)} \right). \] This yields the estimator \[ \widehat \psi = \dfrac{\displaystyle\sum_{i=1}^n \dfrac{Z_i Y_i}{e(X_i)}}{\displaystyle \sum_{i=1}^n \dfrac{Z_i}{e(X_i)}} - \dfrac{\displaystyle \sum_{i=1}^n \dfrac{(1-Z_i) Y_i}{1-e(X_i)}}{\displaystyle\sum_{i=1}^n \dfrac{(1-Z_i)}{1-e(X_i)}}. \] For the unknown propensity score case, the score vector is \[ \bU(\theta) = \begin{pmatrix} R (Y - \beta_0 - Z \psi) \\[6pt] R Z (Y - \beta_0 - Z \psi) \\[6pt] \bx^\top (Z - e(X; \alpha)) \end{pmatrix} = \begin{pmatrix} R \varepsilon_Y \\[6pt] R Z \varepsilon_Y \\[6pt] \bx^\top \epsilon_Z \end{pmatrix} \] Therefore \[ \mathcal{I} = \E \begin{bmatrix} R^2 \varepsilon_Y^2 & R^2 Z \varepsilon_Y^2 & \bx R \epsilon_Z \varepsilon_Y \\[6pt] R^2 Z \varepsilon_Y^2 & R^2 Z^2 \varepsilon_Y^2 & \bx R Z \epsilon_Z \varepsilon_Y \\[6pt] \bx^\top R \epsilon_Z \varepsilon_Y & \bx^\top R Z \epsilon_Z \varepsilon_Y & \bx^\top \bx \epsilon_Z^2 \end{bmatrix} \qquad \qquad \mathcal{J} = \E \begin{bmatrix} R & R Z & -\dot{\bR}^\top \varepsilon_Y \\[6pt] R Z & RZ^2 & -\dot{\bR}^\top Z \varepsilon_Y \\[6pt] \zerovec & \zerovec & \bx^\top \bx e (1-e) \end{bmatrix}. \] where \[ \dot{\bR} = \frac{d R}{d \alpha} = \frac{1-Z}{1-e(X)} e(X) \bx^\top - \frac{Z}{e(X)} (1-e(X)) \bx^\top. \] \end{itemize} <>= w0<-(1-Z)/(1-eX) w1<-Z/eX w<-w0+w1 R<-w fit.ipw<-lm(Y~Z,weights=w) summary(fit.ipw)$coef ipw.ests<-coef(fit.ipw) All.ests[4]<-ipw.ests[2] eps.y<-Y-ipw.ests[1]-ipw.ests[2]*Z eps.z<-Z-eX I.n<-crossprod(cbind(R*eps.y,R*Z*eps.y,Xa*eps.z))/n; round(I.n,5);I.n.ipw2<-I.n Rdot1<-((1-Z)/(1-eX))*eX-(Z/eX)*(1-eX) Rdot2<-X1*Rdot1 Rdot3<-X2*Rdot1 J.n<-matrix(0,5,5) J.n[1,1]<-mean(R) J.n[1,2]<-mean(R*Z) J.n[1,3]<-mean(-Rdot1*eps.y) J.n[1,4]<-mean(-Rdot2*eps.y) J.n[1,5]<-mean(-Rdot3*eps.y) J.n[2,1]<-mean(R*Z) J.n[2,2]<-mean(R*Z^2) J.n[2,3]<-mean(-Rdot1*Z*eps.y) J.n[2,4]<-mean(-Rdot2*Z*eps.y) J.n[2,5]<-mean(-Rdot3*Z*eps.y) J.n[3:5,3:5]<-crossprod(Xa*sqrt(vx))/n; round(J.n,5);J.n.ipw2<-J.n #Asymptotic variance V.n<-solve(J.n) %*% (I.n %*% t(solve(J.n))); round(V.n,5); V.n.ipw2<-V.n #Approximate estimate variance round(V.n/n,5) se.vec<-sqrt(diag(V.n)/n); se.vec All.se[4]<-se.vec[2] @ Therefore the estimate is \Sexpr{round(ipw.ests[2],5)} with standard error \Sexpr{round(se.vec[2],5)}. \hfill{4 Marks} \bigskip \item \textit{An augmented outcome regression model \[ Y_i = \beta_0 + \psi Z_i + \phi_0 \frac{1-Z_i}{1-e(X_{i1},X_{i2};\widehat \alpha)} + \phi_1 \frac{Z_i}{e(X_{i1},X_{i2};\widehat \alpha)} \] with $\alpha$ estimated using the propensity score model from (b).} \medskip \textbf{Solution:} For augmented outcome regression model with estimated propensity score, if \[ R_0 = \frac{1-Z}{1-e(X;\alpha)} \qquad R_1 = \frac{Z}{e(X;\alpha)} \] we have \[ \bU(\theta) = \begin{pmatrix} (Y - \beta_0 - \psi Z - \phi_0 R_0 - \phi_1 R_1)\\[12pt] Z (Y - \beta_0 - \psi Z - \phi_0 R_0 - \phi_1 R_1)\\[12pt] R_0 (Y - \beta_0 - \psi Z - \phi_0 R_0 - \phi_1 R_1)\\[12pt] R_1 (Y - \beta_0 - \psi Z - \phi_0 R_0 - \phi_1 R_1)\\[12pt] \bx^\top (Z - e(X; \alpha)) \end{pmatrix} = \begin{pmatrix} \varepsilon_Y \\[6pt] Z \varepsilon_Y \\[6pt] R_0 \varepsilon_Y \\[6pt] R_1 \varepsilon_Y \\[6pt] \bx^\top \epsilon_Z \end{pmatrix}. \] where $\bx=(1,x_1,x_2)$. We have therefore that \[ \mathcal{I} = \E \begin{bmatrix} \varepsilon_Y^2 & Z \varepsilon_Y^2 & R_0 \varepsilon_Y^2 & R_1 \varepsilon_Y^2 & \bx \epsilon_Z \varepsilon_Y \\[6pt] Z \varepsilon_Y^2 & Z^2 \varepsilon_Y^2 & Z R_0 \varepsilon_Y^2 & Z R_1 \varepsilon_Y^2 & \bx Z \epsilon_Z \varepsilon_Y \\[6pt] R_0 \varepsilon_Y^2 & Z R_0 \varepsilon_Y^2 & R_0^2 \varepsilon_Y^2 & R_0 R_1 \varepsilon_Y^2 & \bx R_0 \epsilon_Z \varepsilon_Y\\[6pt] R_1 \varepsilon_Y^2 & Z R_1 \varepsilon_Y^2 & R_0 R_1 \varepsilon_Y^2 & R_1^2 \varepsilon_Y^2 & \bx R_1 \epsilon_Z \varepsilon_Y\\[6pt] \bx^\top \epsilon_Z \varepsilon_Y & \bx^\top Z \epsilon_Z \varepsilon_Y & \bx^\top R_0 \epsilon_Z \varepsilon_Y & \bx^\top R_0 \epsilon_Z \varepsilon_Y & \bx^\top \bx \epsilon_Z^2 \end{bmatrix} \] and \[ \mathcal{J} = \E \begin{bmatrix} 1 & Z & R_0 & R_1 & \phi_0 \dot{\bR}_0 + \phi_1 \dot{\bR}_1 \\[6pt] Z & Z^2 & Z R_0 & Z R_1 & Z (\phi_0 \dot{\bR}_0 + \phi_1 \dot{\bR}_1) \\[6pt] R_0 & Z R_0 & R_0^2 & R_0 R_1 & R_0 (\phi_0 \dot{\bR}_0 + \phi_1 \dot{\bR}_1) - \dot{\bR}_0 \varepsilon_Y\\[6pt] R_1 & Z R_1 & R_0 R_1 & R_1^2 & R_1 (\phi_0 \dot{\bR}_0 + \phi_1 \dot{\bR}_1) - \dot{\bR}_1 \varepsilon_Y \\[6pt] \zerovec & \zerovec & \zerovec & \zerovec & \bx^\top \bx e (1-e) \end{bmatrix}. \] where \begin{align*} \dot{\bR}_0 & = \frac{\partial }{\partial \alpha } \left\{ \frac{1-Z}{1-e(X;\alpha)} \right\} = \frac{1-Z}{1-e(X;\alpha)} e(X;\alpha) \bx^\top = R_0 e(X;\alpha) \bx^\top \\[6pt] \dot{\bR}_1 & = \frac{\partial }{\partial \alpha} \left\{ \frac{Z}{e(X;\alpha)} \right\} = - \frac{Z}{e(X;\alpha)} (1-e(X;\alpha)) \bx^\top = - R_1 (1-e(X;\alpha)) \bx^\top \end{align*} <>= eX<-fitted(glm(Z~X1+X2,family=binomial,data=P3.data)) vx<-eX*(1-eX) Ptmp<-P3.data Ptmp$eX<-eX R0<-(1-Z)/(1-eX) R1<-Z/eX Ptmp$R0<-R0 Ptmp$R1<-R1 Ptmp$R<-R0+R1 fit.aor<-lm(Y~Z+R0+R1,data=Ptmp) aor.ests<-coef(fit.aor) All.ests[5]<-aor.ests[2] Xa<-cbind(1,X1,X2) eps.y<-residuals(fit.aor) eps.z<-Z-eX I.n<-crossprod(cbind(eps.y,Z*eps.y,R0*eps.y,R1*eps.y,Xa*eps.z))/n; I.n.aor<-I.n round(I.n,5) J.n<-matrix(0,7,7) Xp<-cbind(1,Z,R0,R1) dR0<-R0*eX dR1<--R1*(1-eX) dterm<-(aor.ests[3]*dR0+aor.ests[4]*dR1) J.n[1:4,1:4]<-crossprod(Xp)/n J.n[5:7,5:7]<-crossprod(Xa*sqrt(vx))/n J.n[1,5]<-mean(dterm) J.n[1,6]<-mean(dterm*X1) J.n[1,7]<-mean(dterm*X2) J.n[2,5]<-mean(Z*dterm) J.n[2,6]<-mean(Z*dterm*X1) J.n[2,7]<-mean(Z*dterm*X2) J.n[3,5]<-mean((R0*dterm-dR0*eps.y)) J.n[3,6]<-mean((R0*dterm-dR0*eps.y)*X1) J.n[3,7]<-mean((R0*dterm-dR0*eps.y)*X2) J.n[4,5]<-mean((R1*dterm-dR1*eps.y)) J.n[4,6]<-mean((R1*dterm-dR1*eps.y)*X1) J.n[4,7]<-mean((R1*dterm-dR1*eps.y)*X2) round(J.n,5);J.n.aor<-J.n #Asymptotic variance V.n<-solve(J.n) %*% (I.n %*% t(solve(J.n))); round(V.n,5); V.n.aor<-V.n #Approximate estimate variance round(V.n/n,5) se.vec<-sqrt(diag(V.n)/n); se.vec All.se[5]<-se.vec[2] @ Therefore the estimate is \Sexpr{round(aor.ests[2],5)} with standard error \Sexpr{round(se.vec[2],5)}. \hfill{6 Marks} \medskip \textbf{Note:} In this model, the parameter $\psi$ does not quantify the causal effect of $Z$ on $Y$ unless the outcome mean model component is correctly specified. If this were true, we would expect to see the estimates of $\phi_0$ and $\phi_1$ close to zero: <>= round(summary(fit.aor)$coef,5) @ The coefficients of $R_0$ and $R_1$ are seemingly different from zero on this basis of the OLS analysis; using the sandwich form of the asymptotic variance, we have the following equivalent calculation: <>= aor.summary<-matrix(0,4,4) aor.summary[,1]<-aor.ests aor.summary[,2]<-se.vec[1:4] aor.summary[,3]<-aor.ests/se.vec[1:4] aor.summary[,4]<-2*pnorm(-abs(aor.ests/se.vec[1:4])) rownames(aor.summary)<-c('(Intercept)','Z','R0','R1') colnames(aor.summary)<-c('Estimate','Std. Error','t value','Pr(>|t|)') round(aor.summary,5) @ <>= newd0<-Ptmp;newd0$Z<-0;newd0$R0<-1/(1-eX);newd0$R1<-0 newd1<-Ptmp;newd1$Z<-1;newd1$R1<-1/eX;newd1$R0<-0 mu0.est<-mean(predict(fit.aor,newdata=newd0)) mu1.est<-mean(predict(fit.aor,newdata=newd1)) c(mu0.est,mu1.est,mu1.est-mu0.est) All.ests[6]<-mu1.est-mu0.est @ Therefore, using the prediction method based on the fitted AOR, we have that the estimated average treatment effect is \Sexpr{round(mu1.est-mu0.est,5)}. <>= library(sandwich) round(sandwich(fit.aor),5) @ \end{enumerate} \pagebreak \textbf{Additional Results:} We can check the sense of these results using the bootstrap <>= nboot<-10000; ests.mat<-matrix(0,nrow=nboot,ncol=6) for(ib in 1:nboot){ ind<-sample(1:n,size=n,rep=T) Pboot<-P3.data[ind,] ps1<-fitted(glm(Z~X1+X2,family=binomial,data=Pboot)) Pboot$eX<-ps1; Pboot$R0<-(1-Pboot$Z)/(1-ps1); Pboot$R1<-Pboot$Z/ps1 Pboot$R<-Pboot$R0+Pboot$R1 wv<-Pboot$R1-Pboot$R0 ests.mat[ib,1]<-coef(lm(Y~X1+Z,data=Pboot))[3] ests.mat[ib,2]<-coef(lm(Y~Z+eX,data=Pboot))[2] ests.mat[ib,3]<-mean(wv*Pboot$Y) ests.mat[ib,4]<-coef(lm(Y~Z,weights=R,data=Pboot))[2] fit.aor<-lm(Y~Z+R0+R1,data=Pboot) ests.mat[ib,5]<-coef(fit.aor)[2] newd0<-Pboot;newd0$Z<-0;newd0$R0<-1/(1-eX);newd0$R1<-0 newd1<-Pboot;newd1$Z<-1;newd1$R1<-1/eX;newd1$R0<-0 mu0.est<-mean(predict(fit.aor,newdata=newd0)) mu1.est<-mean(predict(fit.aor,newdata=newd1)) ests.mat[ib,6]<-mu1.est-mu0.est } @ <>= colnames(ests.mat)<-c('OR','PSR','IPW1','IPW2','AOR','AOR-ATE') par(mar=c(4,3,1,1)) boxplot(ests.mat,pch=19,cex=0.5) cmat<-rbind(apply(ests.mat,2,mean),All.ests);rownames(cmat)<-c('Boot','Est.');cmat smat<-rbind(apply(ests.mat,2,sd),All.se);rownames(smat)<-c('Boot','Est.');smat @ \textbf{Computing the sandwich estimate numerically:} We can use the \texttt{nleqslv} package to return estimates that solve the estimating equations, and also to compute the $\dot{\bU}$ estimated values in the \texttt{jacobian} argument. \bigskip \textbf{Outcome regression:} <>= library(nleqslv) or.eqn<-function(xp,dfx){ be<-xp[1:3] Xp<-cbind(1,dfx$X1,dfx$Z) r1<-t(Xp) %*% (dfx$Y - Xp %*% be) return(r1) } or.hess<-function(xp,dfx){ be<-xp[1:3] Xp<-cbind(1,dfx$X1,dfx$Z) epsy<-(dfx$Y - Xp %*% be) hmat<-crossprod(cbind(epsy,dfx$X1*epsy,dfx$Z*epsy)) return(hmat) } xs<-rep(0,3) fit.nl<-nleqslv(xs,fn=or.eqn,dfx=Ptmp, jacobian=TRUE,method='Newton',control=list(ftol=1e-10,maxit=500,allowSingular=TRUE)) est<-fit.nl$x J.hat<-fit.nl$jac/n I.hat<-or.hess(est,Ptmp)/n V.hat.nl.or<-solve(J.hat) %*% I.hat %*% t(solve(J.hat))/n @ \bigskip \textbf{Propensity score regression:} <>= library(nleqslv) expit<-function(x){return(1/(1+exp(-x)))} psr.eqn<-function(xp,dfx){ be<-xp[1:3] al<-xp[4:6] Xa<-cbind(1,dfx$X1,dfx$X2) psv<- expit(Xa %*% al) Xp<-cbind(1,dfx$Z,psv) r1<-t(Xp) %*% (dfx$Y - Xp %*% be) r2<-t(Xa) %*% (dfx$Z - psv) return(c(r1,r2)) } psr.hess<-function(xp,dfx){ be<-xp[1:3] al<-xp[4:6] Xa<-cbind(1,dfx$X1,dfx$X2) psv<- expit(Xa %*% al) Xp<-cbind(1,dfx$Z,psv) epsy<-(dfx$Y - Xp %*% be) epsz<-as.vector(dfx$Z - psv) hmat<-crossprod(cbind(epsy,dfx$Z*epsy,psv*epsy,Xa*epsz)) return(hmat) } xs<-rep(0,6) fit.nl<-nleqslv(xs,fn=psr.eqn,dfx=Ptmp, jacobian=TRUE,method='Newton',control=list(ftol=1e-10,maxit=500,allowSingular=TRUE)) est<-fit.nl$x J.hat<-fit.nl$jac/n I.hat<-psr.hess(est,Ptmp)/n V.hat.nl.psr<-solve(J.hat) %*% I.hat %*% t(solve(J.hat))/n @ \bigskip \textbf{IPW version 1:} <>= library(nleqslv) expit<-function(x){return(1/(1+exp(-x)))} ipw1.eqn<-function(xp,dfx){ be<-xp[1] al<-xp[2:4] Xa<-cbind(1,dfx$X1,dfx$X2) psv<- expit(Xa %*% al) Rv<-as.vector(dfx$Z/psv-(1-dfx$Z)/(1-psv)) r1<-rep(1,n)%*%(Rv*dfx$Y-be) r2<-t(Xa) %*% (dfx$Z - psv) return(c(r1,r2)) } ipw1.hess<-function(xp,dfx){ be<-xp[1] al<-xp[2:4] Xa<-cbind(1,dfx$X1,dfx$X2) psv<- expit(Xa %*% al) Rv<-dfx$Z/psv-(1-dfx$Z)/(1-psv) epsy<-Rv*dfx$Y-be epsz<-as.vector(dfx$Z - psv) hmat<-crossprod(cbind(epsy,Xa*epsz)) return(hmat) } xs<-rep(0,4) fit.nl<-nleqslv(xs,fn=ipw1.eqn,dfx=Ptmp, jacobian=TRUE,method='Newton',control=list(ftol=1e-10,maxit=500)) est<-fit.nl$x J.hat<-fit.nl$jac/n I.hat<-ipw1.hess(est,Ptmp)/n V.hat.nl.ipw1<-solve(J.hat) %*% I.hat %*% t(solve(J.hat))/n @ \bigskip \textbf{IPW version 2:} <>= library(nleqslv) expit<-function(x){return(1/(1+exp(-x)))} ipw2.eqn<-function(xp,dfx){ be<-xp[1:2] al<-xp[3:5] Xa<-cbind(1,dfx$X1,dfx$X2) Xp<-cbind(1,dfx$Z) psv<- expit(Xa %*% al) Rv<-(1-dfx$Z)/(1-psv)+dfx$Z/psv r1<-t(Xp) %*% (Rv*(dfx$Y - Xp %*% be)) r2<-t(Xa) %*% (dfx$Z - psv) return(c(r1,r2)) } ipw2.hess<-function(xp,dfx){ be<-xp[1:2] al<-xp[3:5] Xa<-cbind(1,dfx$X1,dfx$X2) Xp<-cbind(1,dfx$Z) psv<- expit(Xa %*% al) Rv<-(1-dfx$Z)/(1-psv)+dfx$Z/psv epsy<-(dfx$Y - Xp %*% be) epsz<-as.vector(dfx$Z - psv) hmat<-crossprod(cbind(Rv*epsy,Rv*dfx$Z*epsy,Xa*epsz)) return(hmat) } xs<-rep(0,5) fit.nl<-nleqslv(xs,fn=ipw2.eqn,dfx=Ptmp, jacobian=TRUE,method='Newton',control=list(ftol=1e-10,maxit=500)) est<-fit.nl$x J.hat<-fit.nl$jac/n I.hat<-ipw2.hess(est,Ptmp)/n V.hat.nl.ipw2<-solve(J.hat) %*% I.hat %*% t(solve(J.hat))/n @ \bigskip \textbf{Augmented Outcome Regression:} <>= library(nleqslv) expit<-function(x){return(1/(1+exp(-x)))} aor.eqn<-function(xp,dfx){ be<-xp[1:4] al<-xp[5:7] Xa<-cbind(1,dfx$X1,dfx$X2) psv<- expit(Xa %*% al) R0v<-(1-dfx$Z)/(1-psv) R1v<-dfx$Z/psv Xp<-cbind(1,dfx$Z,R0v,R1v) r1<-t(Xp) %*% (dfx$Y - Xp %*% be) r2<-t(Xa) %*% (dfx$Z - psv) return(c(r1,r2)) } aor.hess<-function(xp,dfx){ be<-xp[1:4] al<-xp[5:7] Xa<-cbind(1,dfx$X1,dfx$X2) psv<- expit(Xa %*% al) R0v<-(1-dfx$Z)/(1-psv) R1v<-dfx$Z/psv Xp<-cbind(1,dfx$Z,R0v,R1v) epsy<-(dfx$Y - Xp %*% be) epsz<-as.vector(dfx$Z - psv) hmat<-crossprod(cbind(epsy,dfx$Z*epsy,R0v*epsy,R1v*epsy,Xa*epsz)) return(hmat) } xs<-rep(0,7) fit.nl<-nleqslv(xs,fn=aor.eqn,dfx=Ptmp, jacobian=TRUE,method='Newton',control=list(ftol=1e-10,maxit=500,allowSingular=TRUE)) xs<-fit.nl$x fit.nl<-nleqslv(xs,fn=aor.eqn,dfx=Ptmp, jacobian=TRUE,method='Newton',control=list(ftol=1e-10,maxit=500)) est<-fit.nl$x J.hat<-fit.nl$jac/n I.hat<-aor.hess(est,Ptmp)/n V.hat.nl.aor<-solve(J.hat) %*% (I.hat %*% t(solve(J.hat)))/n @ <>= smat<-rbind(apply(ests.mat,2,sd),All.se, c(sqrt(V.hat.nl.or[3,3]),sqrt(V.hat.nl.psr[2,2]),sqrt(V.hat.nl.ipw1[1,1]), sqrt(V.hat.nl.ipw2[2,2]),sqrt(V.hat.nl.aor[2,2]),0)) rownames(smat)<-c('Boot','Est.','NL');smat @ Note that it is sometimes useful to compute some of the expectations analytically if it is possible, assuming correct specification. For example, we have in this binary treatment case that \[ \E[\epsilon_Z^2] = \E_X[ \E_{Z|X}[(Z-e(X))^2|X]] \equiv \E[e(X)(1-e(X))] \] meaning that there are two possible estimators based on the data, namely \[ \frac{1}{n} \sum_{i=1}^n (Z_i - e(X_i))^2 \qquad \text{and} \qquad \frac{1}{n} \sum_{i=1}^n e(X_i) (1-e(X_i)). \] The second estimator assumes that $\E[Z|X] = e(X)$, that is, that we have correct specification of the PS, that is required for IPW to give a consistent result in estimation. \end{document}