% Options for packages loaded elsewhere
\PassOptionsToPackage{unicode}{hyperref}
\PassOptionsToPackage{hyphens}{url}
%
\documentclass[
]{article}
\usepackage{amsmath,amssymb}
\usepackage{lmodern}
\usepackage{iftex}
\ifPDFTeX
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
  \usepackage{textcomp} % provide euro and other symbols
\else % if luatex or xetex
  \usepackage{unicode-math}
  \defaultfontfeatures{Scale=MatchLowercase}
  \defaultfontfeatures[\rmfamily]{Ligatures=TeX,Scale=1}
\fi
% Use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\IfFileExists{microtype.sty}{% use microtype if available
  \usepackage[]{microtype}
  \UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\makeatletter
\@ifundefined{KOMAClassName}{% if non-KOMA class
  \IfFileExists{parskip.sty}{%
    \usepackage{parskip}
  }{% else
    \setlength{\parindent}{0pt}
    \setlength{\parskip}{6pt plus 2pt minus 1pt}}
}{% if KOMA class
  \KOMAoptions{parskip=half}}
\makeatother
\usepackage{xcolor}
\IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available
\IfFileExists{bookmark.sty}{\usepackage{bookmark}}{\usepackage{hyperref}}
\hypersetup{
  pdftitle={Homework 08 Solutions},
  hidelinks,
  pdfcreator={LaTeX via pandoc}}
\urlstyle{same} % disable monospaced font for URLs
\usepackage[margin=1in]{geometry}
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\usepackage{longtable,booktabs,array}
\usepackage{calc} % for calculating minipage widths
% Correct order of tables after \paragraph or \subparagraph
\usepackage{etoolbox}
\makeatletter
\patchcmd\longtable{\par}{\if@noskipsec\mbox{}\fi\par}{}{}
\makeatother
% Allow footnotes in longtable head/foot
\IfFileExists{footnotehyper.sty}{\usepackage{footnotehyper}}{\usepackage{footnote}}
\makesavenoteenv{longtable}
\usepackage{graphicx}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}
% Set default figure placement to htbp
\makeatletter
\def\fps@figure{htbp}
\makeatother
\setlength{\emergencystretch}{3em} % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{-\maxdimen} % remove section numbering
\usepackage{amsmath}
\usepackage{amsthm}
\DeclareMathOperator{\Cov}{Cov}
\DeclareMathOperator{\Corr}{Corr}
\ifLuaTeX
  \usepackage{selnolig}  % disable illegal ligatures
\fi

\title{Homework 08 Solutions}
\author{}
\date{\vspace{-2.5em}11/10/2021}

\begin{document}
\maketitle

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(tidyverse)}
\FunctionTok{library}\NormalTok{(pander)}
\FunctionTok{library}\NormalTok{(nlme)}
\FunctionTok{library}\NormalTok{(janitor)}
\end{Highlighting}
\end{Shaded}

\hypertarget{problem-1}{%
\section{Problem 1}\label{problem-1}}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{data\_mat }\OtherTok{\textless{}{-}} \FunctionTok{matrix}\NormalTok{(}\FunctionTok{c}\NormalTok{(}
  \DecValTok{1}\NormalTok{, }\DecValTok{300}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{4}\NormalTok{, }\DecValTok{6}\NormalTok{,}
  \DecValTok{2}\NormalTok{, }\DecValTok{300}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{4}\NormalTok{, }\DecValTok{6}\NormalTok{,}
  \DecValTok{3}\NormalTok{, }\DecValTok{500}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{4}\NormalTok{, }\DecValTok{6}\NormalTok{,}
  \DecValTok{4}\NormalTok{, }\DecValTok{500}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{4}\NormalTok{, }\DecValTok{6}\NormalTok{,}
  \DecValTok{5}\NormalTok{, }\DecValTok{200}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{10}\NormalTok{, }\DecValTok{12}\NormalTok{,}
  \DecValTok{6}\NormalTok{, }\DecValTok{200}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{10}\NormalTok{, }\DecValTok{12}\NormalTok{,}
  \DecValTok{7}\NormalTok{, }\DecValTok{200}\NormalTok{, }\DecValTok{0}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{10}\NormalTok{, }\DecValTok{12}\NormalTok{,}
  \DecValTok{8}\NormalTok{, }\DecValTok{200}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{10}\NormalTok{, }\DecValTok{12}\NormalTok{), }\AttributeTok{byrow =} \ConstantTok{TRUE}\NormalTok{, }\AttributeTok{nrow =} \DecValTok{8}\NormalTok{)}

\FunctionTok{colnames}\NormalTok{(data\_mat) }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"Category"}\NormalTok{, }\StringTok{"num"}\NormalTok{, }\StringTok{"x"}\NormalTok{,}
                        \StringTok{"T"}\NormalTok{, }\StringTok{"y\_0"}\NormalTok{, }\StringTok{"y\_1"}\NormalTok{)}
\NormalTok{df }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(data\_mat)}
\end{Highlighting}
\end{Shaded}

\hypertarget{a}{%
\subsection{(a)}\label{a}}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{N }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num)}
\NormalTok{ace }\OtherTok{\textless{}{-}} \DecValTok{1}\SpecialCharTok{/}\NormalTok{N}\SpecialCharTok{*}\FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num }\SpecialCharTok{*}\NormalTok{ (df}\SpecialCharTok{$}\NormalTok{y\_1 }\SpecialCharTok{{-}}\NormalTok{ df}\SpecialCharTok{$}\NormalTok{y\_0))}
\NormalTok{ace}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 2
\end{verbatim}

We can actually directly see that for each individual, the difference
between \(y^1\) and \(y^0\) is 2 (which means the average is also this
value), but the above equation also shows us the the ACE is 2.

\hypertarget{b}{%
\subsection{(b)}\label{b}}

If we only view \(x\) as the only pre-treatment predictor (as suggested
in the beginning of the problem), then we can observe that there \(T\)
and \(x\) appear to be independent - aka \(P(T=t,X=x) = P(T=t)P(X=x)\),
this is explored specifically in part (d), but the below tables show
this through showing that \(P(T=t) = P(T=t|x = x)\) for
\(x \in \{0,1\}\).

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tab\_info }\OtherTok{\textless{}{-}}\NormalTok{ df }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{mutate}\NormalTok{(}\AttributeTok{T =} \FunctionTok{paste0}\NormalTok{(}\StringTok{"T"}\NormalTok{,T)) }\SpecialCharTok{\%\textgreater{}\%}
  \FunctionTok{group\_by}\NormalTok{(x,T) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{summarise}\NormalTok{(}\AttributeTok{count =} \FunctionTok{sum}\NormalTok{(num)) }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  tidyr}\SpecialCharTok{::}\FunctionTok{pivot\_wider}\NormalTok{(}\AttributeTok{id\_cols =}\NormalTok{ x, }\AttributeTok{names\_from =}\NormalTok{ T, }\AttributeTok{values\_from =}\NormalTok{ count) }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  janitor}\SpecialCharTok{::}\FunctionTok{adorn\_totals}\NormalTok{(}\AttributeTok{where =} \FunctionTok{c}\NormalTok{(}\StringTok{"row"}\NormalTok{, }\StringTok{"col"}\NormalTok{)) }
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## `summarise()` has grouped output by 'x'. You can override using the `.groups` argument.
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tab\_info}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##      x   T0   T1 Total
##      0  500  700  1200
##      1  500  700  1200
##  Total 1000 1400  2400
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tab\_info }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  janitor}\SpecialCharTok{::}\FunctionTok{adorn\_percentages}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##      x        T0        T1 Total
##      0 0.4166667 0.5833333     1
##      1 0.4166667 0.5833333     1
##  Total 0.4166667 0.5833333     1
\end{verbatim}

With our full fiew of the data we can see that \((y^0, y^1)\) pairing
actually differ between the first four categories and the last four
categories. Given that \(x\) doesn't capture this different it wouldn't
be too crazy to suggest that there are actually more pre-treatment
predictors of the outcome. If one just defines another grouping variable
\texttt{x2} to be \((1,1,1,1,0,0,0,0)\) to capture this we can see that
\texttt{x2} and \texttt{T} are not independent of each other.
None-the-less this goes beyond the problem description.

\hypertarget{c}{%
\subsection{(c)}\label{c}}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{N\_T }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num[df}\SpecialCharTok{$}\NormalTok{T }\SpecialCharTok{==} \DecValTok{1}\NormalTok{])}
\NormalTok{N\_not\_T }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num[df}\SpecialCharTok{$}\NormalTok{T }\SpecialCharTok{==} \DecValTok{0}\NormalTok{])}

\DecValTok{1}\SpecialCharTok{/}\NormalTok{N\_T }\SpecialCharTok{*} \FunctionTok{sum}\NormalTok{((df}\SpecialCharTok{$}\NormalTok{num }\SpecialCharTok{*}\NormalTok{ df}\SpecialCharTok{$}\NormalTok{y\_1)[df}\SpecialCharTok{$}\NormalTok{T }\SpecialCharTok{==} \DecValTok{1}\NormalTok{]) }\SpecialCharTok{{-}} 
  \DecValTok{1}\SpecialCharTok{/}\NormalTok{N\_not\_T }\SpecialCharTok{*} \FunctionTok{sum}\NormalTok{((df}\SpecialCharTok{$}\NormalTok{num }\SpecialCharTok{*}\NormalTok{ df}\SpecialCharTok{$}\NormalTok{y\_0)[df}\SpecialCharTok{$}\NormalTok{T }\SpecialCharTok{==} \DecValTok{0}\NormalTok{])}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 1.314286
\end{verbatim}

If there was no other confounders we'd expect this value to approximate
the ACE. Interestingly it does not, but that seems to be related to the
second commentary we had in part (b), which showed that the proportion
those in the treatment group where \texttt{x2\ =\ 1} was larger than the
proportion of those in the treatment group where \texttt{x2\ =\ 0}, and
this relateds to the smaller value we see for the treatment effect as
those in the \texttt{x2\ =\ 1} group had smaller outcomes under both
treatments than the \texttt{x2\ =\ 0} group.

\hypertarget{d}{%
\subsection{(d)}\label{d}}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prob\_T }\OtherTok{\textless{}{-}}\FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num }\SpecialCharTok{*}\NormalTok{ (df}\SpecialCharTok{$}\NormalTok{T }\SpecialCharTok{==} \DecValTok{1}\NormalTok{)) }\SpecialCharTok{/} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num)}

\NormalTok{prob\_T}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 0.5833333
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{prob\_T\_given\_x1 }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num }\SpecialCharTok{*}\NormalTok{ (df}\SpecialCharTok{$}\NormalTok{x }\SpecialCharTok{==} \DecValTok{1}\NormalTok{) }\SpecialCharTok{*}\NormalTok{ (df}\SpecialCharTok{$}\NormalTok{T }\SpecialCharTok{==} \DecValTok{1}\NormalTok{)) }\SpecialCharTok{/} 
  \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num }\SpecialCharTok{*}\NormalTok{ (df}\SpecialCharTok{$}\NormalTok{x }\SpecialCharTok{==} \DecValTok{1}\NormalTok{))}

\NormalTok{prob\_T\_given\_x1}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 0.5833333
\end{verbatim}

We find that, in the above population, \(P(T=1) =\) 0.5833, which is the
same as \(P(T=1|x=1)\).

\hypertarget{e}{%
\subsection{(e)}\label{e}}

Below are 2 approaches to estimate the desired coefficient and se. Note
that the coefficient values for the intercept and \(T\) directly related
to parts of the calculation in part (c). The intercept being the
\texttt{mean(y{[}T=0{]})} and the coefficient of \(T\) being the
difference between \texttt{mean(y{[}T=1{]})\ -\ mean(y{[}T=0{]})}.

\hypertarget{weighted-regression-approach}{%
\subsubsection{Weighted regression
approach}\label{weighted-regression-approach}}

A weighted regression solution is faster that expanding out the full
dataset, but requires a correct ion for the model's degrees of freedom.
See this
\href{https://www.r-bloggers.com/2015/09/linear-models-with-weighted-observations/}{blogpost}
for more comments.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df }\OtherTok{\textless{}{-}}\NormalTok{ df }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{mutate}\NormalTok{(}\AttributeTok{y =} \FunctionTok{ifelse}\NormalTok{(T, y\_1, y\_0))}

\NormalTok{wls }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(y }\SpecialCharTok{\textasciitilde{}}\NormalTok{ T }\SpecialCharTok{+}\NormalTok{ x, }\AttributeTok{data =}\NormalTok{ df, }\AttributeTok{weights =}\NormalTok{ num)}
\NormalTok{wls}\SpecialCharTok{$}\NormalTok{df.residual }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num) }\SpecialCharTok{{-}} \FunctionTok{length}\NormalTok{(wls}\SpecialCharTok{$}\NormalTok{coefficients) }
\CommentTok{\#\^{} correction for number of observations / degrees of freedom}

\NormalTok{summ\_w }\OtherTok{\textless{}{-}} \FunctionTok{summary}\NormalTok{(wls)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in summary.lm(wls): residual degrees of freedom in object suggest this
## is not an "lm" fit
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{summ\_w}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = y ~ T + x, data = df, weights = num)
## 
## Weighted Residuals:
##    Min     1Q Median     3Q    Max 
## -41.57 -39.14   6.29  53.34  60.61 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.4000     0.1058   60.51   <2e-16 ***
## T             1.3143     0.1163   11.30   <2e-16 ***
## x             0.0000     0.1147    0.00        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.81 on 2397 degrees of freedom
## Multiple R-squared:  0.05055,    Adjusted R-squared:  0.9972 
## F-statistic: 63.81 on 2 and 2397 DF,  p-value: < 2.2e-16
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{coef}\NormalTok{(summ\_w)[}\DecValTok{2}\NormalTok{,]}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##     Estimate   Std. Error      t value     Pr(>|t|) 
## 1.314286e+00 1.163413e-01 1.129682e+01 7.202413e-29
\end{verbatim}

\hypertarget{individualized-approach}{%
\subsubsection{Individualized approach}\label{individualized-approach}}

Another way would just be to create an ``individual'' data frame (where
each row represents one individual). Below is a messy way to create the
desired data.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df2 }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{x =} \FunctionTok{unlist}\NormalTok{(}\FunctionTok{sapply}\NormalTok{(}\DecValTok{1}\SpecialCharTok{:}\FunctionTok{nrow}\NormalTok{(df), }
                                    \ControlFlowTok{function}\NormalTok{(idx)\{}\FunctionTok{rep}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{x[idx], }
                                                      \AttributeTok{each =}\NormalTok{ df}\SpecialCharTok{$}\NormalTok{num[idx])\})),}
                  \AttributeTok{T =} \FunctionTok{unlist}\NormalTok{(}\FunctionTok{sapply}\NormalTok{(}\DecValTok{1}\SpecialCharTok{:}\FunctionTok{nrow}\NormalTok{(df), }
                                    \ControlFlowTok{function}\NormalTok{(idx)\{}\FunctionTok{rep}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{T[idx], }
                                                      \AttributeTok{each =}\NormalTok{ df}\SpecialCharTok{$}\NormalTok{num[idx])\})),}
                  \AttributeTok{y =} \FunctionTok{unlist}\NormalTok{(}\FunctionTok{sapply}\NormalTok{(}\DecValTok{1}\SpecialCharTok{:}\FunctionTok{nrow}\NormalTok{(df), }
                                    \ControlFlowTok{function}\NormalTok{(idx)\{}\FunctionTok{rep}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{y[idx], }
                                                      \AttributeTok{each =}\NormalTok{ df}\SpecialCharTok{$}\NormalTok{num[idx])\})))}

\NormalTok{ls\_all\_obs }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(y }\SpecialCharTok{\textasciitilde{}}\NormalTok{ T }\SpecialCharTok{+}\NormalTok{ x, }\AttributeTok{data =}\NormalTok{ df2)}
\NormalTok{summ }\OtherTok{\textless{}{-}} \FunctionTok{summary}\NormalTok{(ls\_all\_obs)}
\NormalTok{summ}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = y ~ T + x, data = df2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.400 -1.886 -1.714  3.600  4.286 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.400e+00  1.058e-01   60.51   <2e-16 ***
## T           1.314e+00  1.163e-01   11.30   <2e-16 ***
## x           2.901e-15  1.147e-01    0.00        1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.81 on 2397 degrees of freedom
## Multiple R-squared:  0.05055,    Adjusted R-squared:  0.04976 
## F-statistic: 63.81 on 2 and 2397 DF,  p-value: < 2.2e-16
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{coef}\NormalTok{(summ)[}\DecValTok{2}\NormalTok{,]}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##     Estimate   Std. Error      t value     Pr(>|t|) 
## 1.314286e+00 1.163413e-01 1.129682e+01 7.202413e-29
\end{verbatim}

\hypertarget{problem-2.}{%
\section{Problem 2.}\label{problem-2.}}

\hypertarget{a-1}{%
\subsection{(a)}\label{a-1}}

The treatment variable is no longer randomized in the way intended by
the researcher, because individuals have selected their own level of
treatment. This selection process is very likely confounded with the
outcome variable, meaning that there is some unobserved variable that is
correlated with both the treatment variable and the outcome.

For example, imagine that the therapy actually has no effect, but there
is an ``optimism'' variable that determines both how many sessions a
person attends and their emotional state at the end of the study. In
other words, the people who engage the most with the treatment are those
in the best state at the end of the study, but not because of the
therapy. Then the number of sessions attended will strongly predict the
outcome variable, but not because of an effect of the therapy.

An unobserved confounder could also obscure a real treatment effect.
Suppose that instead of an optimism variable, there is a variable
representing emotional wellbeing that is inversely correlated with the
number of sessions attended. That is, those who start out the most
depressed seek the most treatment, perhaps because they feel that they
have the most to gain. Suppose the treatment \emph{does} have a
beneficial effect that increases with the number of sessions attended,
and suppose that everyone seeks out just enough treatment to get them to
roughly the same level of emotional wellbeing at the end of the study.
Then a coefficient for number of treatment sessions could end up being 0
even though there is a treatment effect!

\emph{In sum:} whenever there is noncompliance in an experiment, there
is potential for confounding due to selection effects. Including a
variable for the amount of treatment received \emph{does not} address
this and can lead to highly misleading results.

\hypertarget{b-1}{%
\subsection{(b)}\label{b-1}}

This is perfectly set up for an Instrumental Variable analysis. The
assignment variable (assignment to treatment or control) is the
instrument, and the number of sessions attended is the treatment
variable. The assignment variable almost certainly affects the treatment
variable, but it should have no effect on the outcome through other
pathways. (That is, the mere fact of having been assigned to the
treatment or control group shouldn't affect participants' emotional
wellbeing at the end of the study.) An experimental design with
noncompliance is almost always a good candidate for an IV analysis.

An analysis based on propensity score matching won't work here unless we
have variables with which to construct the propensity score. The
propensity score is the probability of receiving treatment as a function
of some set of covariates \(X\). That is, the propensity score is a
function \(\pi(x) := P(A = 1|X = x)\), where \(A = 1\) means that
someone receives treatment. We can define it more broadly as the
probability of receiving a particular amount of treatment \(a\):
\(\pi_a(x) = P(A = a|X= x)\). We want as rich a set of covariates \(X\)
as possible in order to perform matching. If the investigator collected
these prior to the start of the study, then this analysis may be
possible.

\hypertarget{problem-3}{%
\section{Problem 3}\label{problem-3}}

\hypertarget{a-2}{%
\subsection{(a)}\label{a-2}}

A summary of the model is below. Diagnostic plots are in Figure
\ref{f:3a_diagnostic}, and the estimated autocorrelation between the
residuals is plotted in Figure \ref{f:3a_acf}.

The residuals are very roughly normally distributed, but they appear to
show an increasing trend relative to the fitted values. The
scale-location plot also appears to show that the variance is increasing
with fitted values, and the autocorrelation plot shows large
autocorrelations, all of which suggest that the linear model with
assumed independence between the residuals is not appropriate.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{carlsen }\OtherTok{\textless{}{-}} \FunctionTok{read.table}\NormalTok{(}\StringTok{"CarlsenQ.txt"}\NormalTok{, }\AttributeTok{header=}\ConstantTok{TRUE}\NormalTok{)}
\NormalTok{m1 }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(Sales }\SpecialCharTok{\textasciitilde{}}\NormalTok{ . }\SpecialCharTok{{-}}\NormalTok{Case }\SpecialCharTok{{-}}\NormalTok{ Time, }\AttributeTok{data=}\NormalTok{carlsen)}
\FunctionTok{pander}\NormalTok{(}\FunctionTok{summary}\NormalTok{(m1), }\AttributeTok{caption =} \StringTok{"Problem 3a: Summary of ordinary linear model."}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{longtable}[]{@{}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.25}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.15}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.18}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.14}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.15}}@{}}
\toprule
~ & Estimate & Std. Error & t value &
Pr(\textgreater\textbar t\textbar) \\
\midrule
\endhead
\textbf{(Intercept)} & 4421 & 1235 & 3.58 & 0.00336 \\
\textbf{Temp} & 78.18 & 33.89 & 2.307 & 0.03817 \\
\textbf{Sun} & 2.585 & 1.556 & 1.661 & 0.1206 \\
\textbf{Q2} & -879.7 & 1121 & -0.7845 & 0.4468 \\
\textbf{Q3} & -1552 & 1541 & -1.007 & 0.3323 \\
\textbf{Q4} & 715.4 & 474.3 & 1.508 & 0.1554 \\
\bottomrule
\end{longtable}

\begin{longtable}[]{@{}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.21}}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.31}}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.12}}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.24}}@{}}
\caption{Problem 3a: Summary of ordinary linear model.}\tabularnewline
\toprule
Observations & Residual Std. Error & \(R^2\) & Adjusted \(R^2\) \\
\midrule
\endfirsthead
\toprule
Observations & Residual Std. Error & \(R^2\) & Adjusted \(R^2\) \\
\midrule
\endhead
19 & 287.2 & 0.9537 & 0.9359 \\
\bottomrule
\end{longtable}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\# summary(m1)}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow =} \FunctionTok{c}\NormalTok{(}\DecValTok{2}\NormalTok{, }\DecValTok{2}\NormalTok{),  }\AttributeTok{mar =} \FunctionTok{c}\NormalTok{(}\DecValTok{4}\NormalTok{, }\DecValTok{5}\NormalTok{, }\DecValTok{2}\NormalTok{, }\DecValTok{1}\NormalTok{))}
\FunctionTok{plot}\NormalTok{(m1)}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw08_solutions_files/figure-latex/3a_diagnostic-1.pdf}
\caption{Problem 3a: Diagnostic plots for ordinary linear
model.\label{f:3a_diagnostic}}
\end{figure}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{acf}\NormalTok{(}\FunctionTok{resid}\NormalTok{(m1), }\AttributeTok{main =} \StringTok{""}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{figure}

{\centering \includegraphics{hw08_solutions_files/figure-latex/3a_acf-1} 

}

\caption{Problem 3a: acf plot for ordinary linear model residuals.\label{f:3a_acf}}\label{fig:3a_acf}
\end{figure}

\pagebreak

\hypertarget{b-2}{%
\subsection{(b)}\label{b-2}}

Note that the model can be estimated with either Maximum Likelihood (ML)
or Restricted Maximum Likelihood (REML). I have used the default, which
is REML, but either one is fine. They will produce different estimates.
The code should look like one of the following:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{model }\OtherTok{\textless{}{-}} \FunctionTok{gls}\NormalTok{(Sales }\SpecialCharTok{\textasciitilde{}}\NormalTok{ . }\SpecialCharTok{{-}}\NormalTok{Case }\SpecialCharTok{{-}}\NormalTok{ Time, }\AttributeTok{data=}\NormalTok{carlsen, }\AttributeTok{correlation=}\FunctionTok{corAR1}\NormalTok{())}
\NormalTok{model }\OtherTok{\textless{}{-}} \FunctionTok{gls}\NormalTok{(Sales }\SpecialCharTok{\textasciitilde{}}\NormalTok{ . }\SpecialCharTok{{-}}\NormalTok{Case }\SpecialCharTok{{-}}\NormalTok{ Time, }\AttributeTok{data=}\NormalTok{carlsen, }\AttributeTok{correlation=}\FunctionTok{corAR1}\NormalTok{(),}
             \AttributeTok{method =} \StringTok{"ML"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{enumerate}
\def\labelenumi{\roman{enumi}.}
\item
  A summary of this model is below. From this summary, the estimated
  value of \(\rho\), the lag-1 autocorrelation parameter, is 0.87. (If
  you use ML for estimation, the estimate is 0.79).
\item
  The coefficient values have changed in magnitude, though not in sign.
  The estimated standard errors in the AR1 model are much smaller than
  in the ordinary linear model, which is to be expected if the model is
  properly accounting for correlation among the residuals. As a result,
  both \texttt{Sun} and \texttt{Q4} are significant in the AR1 model but
  not in the ordinary linear model at conventional thresholds.
\end{enumerate}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{m2 }\OtherTok{\textless{}{-}} \FunctionTok{gls}\NormalTok{(Sales }\SpecialCharTok{\textasciitilde{}}\NormalTok{ . }\SpecialCharTok{{-}}\NormalTok{Case }\SpecialCharTok{{-}}\NormalTok{ Time, }\AttributeTok{data=}\NormalTok{carlsen, }\AttributeTok{correlation=}\FunctionTok{corAR1}\NormalTok{())}
\FunctionTok{cat}\NormalTok{(}\StringTok{"Problem 3b: Summary of AR1 model:"}\NormalTok{); }\FunctionTok{summary}\NormalTok{(m2)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Problem 3b: Summary of AR1 model:
\end{verbatim}

\begin{verbatim}
## Generalized least squares fit by REML
##   Model: Sales ~ . - Case - Time 
##   Data: carlsen 
##        AIC      BIC    logLik
##   212.9427 217.4623 -98.47134
## 
## Correlation Structure: AR(1)
##  Formula: ~1 
##  Parameter estimate(s):
##       Phi 
## 0.8688414 
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 5215.269  747.4869  6.977071  0.0000
## Temp          55.545   22.7958  2.436611  0.0300
## Sun            2.192    0.7802  2.809438  0.0148
## Q2           -36.255  662.9242 -0.054690  0.9572
## Q3          -371.191  941.7222 -0.394162  0.6999
## Q4           946.075  295.3927  3.202772  0.0069
## 
##  Correlation: 
##      (Intr) Temp   Sun    Q2     Q3    
## Temp -0.870                            
## Sun   0.005 -0.357                     
## Q2    0.914 -0.938  0.033              
## Q3    0.912 -0.960  0.100  0.992       
## Q4    0.747 -0.937  0.485  0.840  0.876
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.1843635 -0.2908069  0.4983337  0.5709198  1.1275844 
## 
## Residual standard error: 409.3627 
## Degrees of freedom: 19 total; 13 residual
\end{verbatim}

\hypertarget{c-1}{%
\subsection{(c)}\label{c-1}}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{T }\OtherTok{\textless{}{-}} \SpecialCharTok{{-}}\DecValTok{2}\SpecialCharTok{*}\NormalTok{(}\FunctionTok{as.numeric}\NormalTok{(}\FunctionTok{logLik}\NormalTok{(m1)) }\SpecialCharTok{{-}} \FunctionTok{as.numeric}\NormalTok{(}\FunctionTok{logLik}\NormalTok{(m2)))}
\NormalTok{pval }\OtherTok{\textless{}{-}} \FunctionTok{pchisq}\NormalTok{(T, }\AttributeTok{df =} \DecValTok{1}\NormalTok{, }\AttributeTok{lower.tail =} \ConstantTok{FALSE}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

The log likelihoods are -131 for the ordinary linear model and -98 for
the AR1 mode. The difference in degrees of freedom is 1, so we compare
the statistic \(-2(-131-(-98)) = 66\) to a \(\chi_{(1)}^2\)
distribution. The pvalue is \(\ensuremath{8.0841612\times 10^{-16}}\),
so we reject the null hypothesis and conclude that the autocorrelation
is not 0.

\hypertarget{d-1}{%
\subsection{(d)}\label{d-1}}

A summary of the model is below. Diagnostic plots are in Figure
\ref{f:3d_diagnostic}, and the autocorrelation plot of the residuals is
in Figure \ref{f:3d_acf}.

The diagnostic plots suggest problems with the fit. The residuals are
not very close to normally distributed, and the scale-location plot
appears to show a relationship between the fitted values and the
variance. The autocorrelation plot, however, shows much reduced
autocorrelation among the residuals, as hoped for.

The likelihood of the model is \texttt{-130}, as compared to
\texttt{-98} for the model in 3c. It appears to be more effective to
estimate the autocorrelation simultaneously with estimating the model,
as in 2c, rather than following a two-stage estimation procedure.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{r }\OtherTok{\textless{}{-}} \FunctionTok{resid}\NormalTok{(m1)}
\NormalTok{rho }\OtherTok{\textless{}{-}} \FunctionTok{cor}\NormalTok{(r[}\SpecialCharTok{{-}}\DecValTok{1}\NormalTok{], r[}\SpecialCharTok{{-}}\FunctionTok{length}\NormalTok{(r)])}
\NormalTok{Sigma }\OtherTok{\textless{}{-}} \FunctionTok{diag}\NormalTok{(}\FunctionTok{nrow}\NormalTok{(carlsen))}
\NormalTok{Sigma }\OtherTok{\textless{}{-}}\NormalTok{ rho}\SpecialCharTok{\^{}}\FunctionTok{abs}\NormalTok{(}\FunctionTok{row}\NormalTok{(Sigma) }\SpecialCharTok{{-}} \FunctionTok{col}\NormalTok{(Sigma))}
\NormalTok{S }\OtherTok{\textless{}{-}} \FunctionTok{chol}\NormalTok{(Sigma)}
\NormalTok{S\_inv }\OtherTok{\textless{}{-}} \FunctionTok{solve}\NormalTok{(}\FunctionTok{t}\NormalTok{(S))}

\NormalTok{Xstar }\OtherTok{\textless{}{-}}\NormalTok{ S\_inv }\SpecialCharTok{\%*\%} \FunctionTok{model.matrix}\NormalTok{(m1)}
\NormalTok{ystar }\OtherTok{\textless{}{-}}\NormalTok{ S\_inv }\SpecialCharTok{\%*\%}\NormalTok{ carlsen}\SpecialCharTok{$}\NormalTok{Sales}

\NormalTok{m3 }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(ystar }\SpecialCharTok{\textasciitilde{}}\NormalTok{ Xstar }\SpecialCharTok{{-}} \DecValTok{1}\NormalTok{)}
\FunctionTok{pander}\NormalTok{(}\FunctionTok{summary}\NormalTok{(m3), }
       \AttributeTok{caption =} \FunctionTok{paste}\NormalTok{(}\StringTok{"Problem 2d: Summary of ordinary linear model"}\NormalTok{,}
                       \StringTok{"with transformed variables."}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{longtable}[]{@{}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.32}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.15}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.18}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.14}}
  >{\centering\arraybackslash}p{(\columnwidth - 8\tabcolsep) * \real{0.15}}@{}}
\toprule
~ & Estimate & Std. Error & t value &
Pr(\textgreater\textbar t\textbar) \\
\midrule
\endhead
\textbf{Xstar(Intercept)} & 5090 & 793.5 & 6.415 & 2.29e-05 \\
\textbf{XstarTemp} & 61.01 & 24.64 & 2.476 & 0.02781 \\
\textbf{XstarSun} & 2.196 & 0.8827 & 2.488 & 0.02722 \\
\textbf{XstarQ2} & -214.5 & 725.8 & -0.2956 & 0.7722 \\
\textbf{XstarQ3} & -631.1 & 1026 & -0.6148 & 0.5493 \\
\textbf{XstarQ4} & 879.1 & 322.7 & 2.724 & 0.01738 \\
\bottomrule
\end{longtable}

\begin{longtable}[]{@{}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.21}}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.31}}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.12}}
  >{\centering\arraybackslash}p{(\columnwidth - 6\tabcolsep) * \real{0.24}}@{}}
\caption{Problem 2d: Summary of ordinary linear model with transformed
variables.}\tabularnewline
\toprule
Observations & Residual Std. Error & \(R^2\) & Adjusted \(R^2\) \\
\midrule
\endfirsthead
\toprule
Observations & Residual Std. Error & \(R^2\) & Adjusted \(R^2\) \\
\midrule
\endhead
19 & 277.7 & 0.9979 & 0.9969 \\
\bottomrule
\end{longtable}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow =} \FunctionTok{c}\NormalTok{(}\DecValTok{2}\NormalTok{, }\DecValTok{2}\NormalTok{),  }\AttributeTok{mar =} \FunctionTok{c}\NormalTok{(}\DecValTok{4}\NormalTok{, }\FloatTok{4.5}\NormalTok{, }\DecValTok{2}\NormalTok{, }\DecValTok{1}\NormalTok{))}
\FunctionTok{plot}\NormalTok{(m3)}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw08_solutions_files/figure-latex/3d_diagnostics-1.pdf}
\caption{Problem 3d: diagnostic plots.\label{f:3d_diagnostic}}
\end{figure}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{acf}\NormalTok{(}\FunctionTok{resid}\NormalTok{(m3), }\AttributeTok{main =} \StringTok{""}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw08_solutions_files/figure-latex/3d_acf-1.pdf}
\caption{Problem 3d: acf plot for ordinary linear model
residuals.\label{f:3d_acf}}
\end{figure}

\end{document}
