% 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}
\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{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
\ifLuaTeX
  \usepackage{selnolig}  % disable illegal ligatures
\fi
\IfFileExists{bookmark.sty}{\usepackage{bookmark}}{\usepackage{hyperref}}
\IfFileExists{xurl.sty}{\usepackage{xurl}}{} % add URL line breaks if available
\urlstyle{same} % disable monospaced font for URLs
\hypersetup{
  pdftitle={Homework 06 Solutions},
  hidelinks,
  pdfcreator={LaTeX via pandoc}}

\title{Homework 06 Solutions}
\author{}
\date{\vspace{-2.5em}2022-10-25}

\begin{document}
\maketitle

\def\red#1{\textcolor{red}{#1}}
\def\blue#1{\textcolor{blue}{#1}}
\let\bar\overline
\let\ul\underline

\frenchspacing

\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
Solutions
\end{center}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# library(arm)  \#\# includes lme4, MASS, Matrix}

\FunctionTok{library}\NormalTok{(ggplot2); }\FunctionTok{theme\_set}\NormalTok{(}\FunctionTok{theme\_bw}\NormalTok{())}
\FunctionTok{library}\NormalTok{(gridExtra) }\DocumentationTok{\#\# to arrange ggplots...}

\DocumentationTok{\#\# library(GGally) \#\# for ggpairs...}


\FunctionTok{library}\NormalTok{(leaps) }\DocumentationTok{\#\# regsubsets(), summary(), coef()}

\DocumentationTok{\#\# library(car)   \#\# subsets(), mmps(), vif(), etc.}

\DocumentationTok{\#\# library(marginalmodelplots) \#\# because mmps doesn\textquotesingle{}t work for glm\textquotesingle{}s...}
\DocumentationTok{\#\# library(DHARMa)  \#\# for better glm residual plots...}

\DocumentationTok{\#\# library(tidyverse) \#\# at last...}

\DocumentationTok{\#\# library(foreign)  \#\# for "read.dta" for problem 3}

\FunctionTok{library}\NormalTok{(ISLR2)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning: package 'ISLR2' was built under R version 4.1.3
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(splines)}

\FunctionTok{library}\NormalTok{(boot)}

\FunctionTok{library}\NormalTok{(mgcv)}

\FunctionTok{set.seed}\NormalTok{(}\DecValTok{1000000}\NormalTok{)  }\DocumentationTok{\#\# repeatable "random" selection of CV subsets...}
\end{Highlighting}
\end{Shaded}

\hypertarget{problem-1-islr-9-p.-324.}{%
\section{Problem 1: ISLR \#9, p.~324.}\label{problem-1-islr-9-p.-324.}}

This question uses the variables dis (the weighted mean of distances to
five Boston employment centers) and nox (nitrogen oxides concentration
in parts per 10 million) from the Boston data. We will treat dis as the
predictor and nox as the response.

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

Use the poly() function to fit a cubic polynomial regression to predict
nox using dis. Report the regression output, and plot the resulting data
and polynomial fits.

\em

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{data}\NormalTok{(Boston)}

\FunctionTok{str}\NormalTok{(Boston)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 'data.frame':    506 obs. of  13 variables:
##  $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
##  $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
##  $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
##  $ chas   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
##  $ rm     : num  6.58 6.42 7.18 7 7.15 ...
##  $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
##  $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
##  $ rad    : int  1 2 2 3 3 3 5 5 5 5 ...
##  $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
##  $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
##  $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
##  $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# help(Boston)   {-}{-} to get the definitions of the variables}

\NormalTok{lm}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{poly}\NormalTok{(dis,}\DecValTok{3}\NormalTok{), }\AttributeTok{data=}\NormalTok{Boston)}
\FunctionTok{summary}\NormalTok{(lm}\FloatTok{.1}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = nox ~ poly(dis, 3), data = Boston)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.121130  -0.040619  -0.009738   0.023385   0.194904  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    0.554695   0.002759 201.021  < 2e-16 ***
## poly(dis, 3)1 -2.003096   0.062071 -32.271  < 2e-16 ***
## poly(dis, 3)2  0.856330   0.062071  13.796  < 2e-16 ***
## poly(dis, 3)3 -0.318049   0.062071  -5.124 4.27e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.003852802)
## 
##     Null deviance: 6.7810  on 505  degrees of freedom
## Residual deviance: 1.9341  on 502  degrees of freedom
## AIC: -1370.9
## 
## Number of Fisher Scoring iterations: 2
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{plot}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis, }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{fits }\OtherTok{\textless{}{-}} \FunctionTok{predict}\NormalTok{(lm}\FloatTok{.1}\NormalTok{)}
\NormalTok{ord }\OtherTok{\textless{}{-}} \FunctionTok{order}\NormalTok{(Boston}\SpecialCharTok{$}\NormalTok{dis)}
\FunctionTok{lines}\NormalTok{(fits[ord] }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis[ord], }\AttributeTok{data=}\NormalTok{Boston)}
\end{Highlighting}
\end{Shaded}

\begin{center}\includegraphics{hw06-solutions_files/figure-latex/unnamed-chunk-2-1} \end{center}

\rm

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

Plot the polynomial fits for a range of different polynomial degrees
(say, from 1 to 10), and report the associated residual sum of squares.

\em

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow=}\FunctionTok{c}\NormalTok{(}\DecValTok{4}\NormalTok{,}\DecValTok{3}\NormalTok{))}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{10}\NormalTok{) \{}
\NormalTok{  lm.fit }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{poly}\NormalTok{(dis,i), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{  rss }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(}\FunctionTok{residuals}\NormalTok{(lm.fit)}\SpecialCharTok{\^{}}\DecValTok{2}\NormalTok{)}
  \FunctionTok{plot}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis, }\AttributeTok{data=}\NormalTok{Boston, }\AttributeTok{col=}\StringTok{"grey"}\NormalTok{,}
       \AttributeTok{main=}\FunctionTok{paste0}\NormalTok{(}\StringTok{"Polynomial of degree "}\NormalTok{,i,}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{RSS = "}\NormalTok{,}\FunctionTok{round}\NormalTok{(rss,}\DecValTok{2}\NormalTok{)))}
\NormalTok{  fits }\OtherTok{\textless{}{-}} \FunctionTok{predict}\NormalTok{(lm.fit)}
\NormalTok{  ord }\OtherTok{\textless{}{-}} \FunctionTok{order}\NormalTok{(Boston}\SpecialCharTok{$}\NormalTok{dis)}
  \FunctionTok{lines}\NormalTok{(fits[ord] }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis[ord], }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{\}}
\end{Highlighting}
\end{Shaded}

\begin{center}\includegraphics{hw06-solutions_files/figure-latex/unnamed-chunk-3-1} \end{center}

\rm

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

Perform cross-validation or another approach to select the optimal
degree for the polynomial, and explain your results.

\em

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# From help(cv.glm, package=boot) we see that }
\DocumentationTok{\#\# (1) cv.glm uses mean squared prediction error as the cross{-}validation criterion}
\DocumentationTok{\#\# (2) We want $delta[1] as the raw cross{-}validation mean squared prediction error}

\DocumentationTok{\#\# First we try LOOCV...}

\NormalTok{cv.err }\OtherTok{\textless{}{-}} \ConstantTok{NULL}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{10}\NormalTok{) \{}
\NormalTok{  lm.fit }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{poly}\NormalTok{(dis,i), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{  cv.err }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(cv.err,}\FunctionTok{cv.glm}\NormalTok{(Boston,lm.fit)}\SpecialCharTok{$}\NormalTok{delta[}\DecValTok{1}\NormalTok{])  }
\NormalTok{\}}
\NormalTok{cv.err}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] 0.005523868 0.004079449 0.003874762 0.003887521 0.004164865 0.005384278
##  [7] 0.011068782 0.008121397 0.017616356 0.004430276
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"Optimal LOOCV polynomial degree ="}\NormalTok{,}\FunctionTok{which}\NormalTok{(cv.err}\SpecialCharTok{==}\FunctionTok{min}\NormalTok{(cv.err)),}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Optimal LOOCV polynomial degree = 3
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# Now let\textquotesingle{}s try 10{-}fold cross{-}validation...}

\NormalTok{cv.err }\OtherTok{\textless{}{-}} \ConstantTok{NULL}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{10}\NormalTok{) \{}
\NormalTok{  lm.fit }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{poly}\NormalTok{(dis,i), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{  cv.err }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(cv.err,}\FunctionTok{cv.glm}\NormalTok{(Boston,lm.fit, }\AttributeTok{K=}\DecValTok{10}\NormalTok{)}\SpecialCharTok{$}\NormalTok{delta[}\DecValTok{1}\NormalTok{])  }
\NormalTok{\}}
\NormalTok{cv.err}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] 0.005501407 0.004084270 0.003896280 0.003879808 0.004152568 0.005817020
##  [7] 0.010192597 0.011093052 0.017637806 0.004047557
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"Optimal 10{-}fold CV polynomial degree ="}\NormalTok{,}\FunctionTok{which}\NormalTok{(cv.err}\SpecialCharTok{==}\FunctionTok{min}\NormalTok{(cv.err)),}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Optimal 10-fold CV polynomial degree = 4
\end{verbatim}

Leave-one-out cross validation (LOOCV) does not depend on a random
selection of folds; LOOCV uses all \(n\) folds with training set of size
\(n-1\) and test set of size 1 (the one observation we leave out). On
the other hand \(K\)-fold cross-validation will depend somewhat on what
random set of folds get selected.

Here, LOOCV suggests that the optimal polynomial fit has degree 3, and
10-fold cross-validation suggests degree 4. While LOOCV will always
choose degree 3 for this data (it is not random), the degree chosen by
10-fold cross-validation may change, depending on the particular random
folds that are chosen. We will use degree 3 in the rest of this problem.

\rm

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

Use the bs() function to fit a regression spline to predict nox using
dis. Report the output for the fit using four degrees of freedom. How
did you choose the knots? Plot the resulting fit.

\em 

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lm}\FloatTok{.2} \OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\DecValTok{4}\NormalTok{),}\AttributeTok{data=}\NormalTok{Boston)}
\FunctionTok{summary}\NormalTok{(lm}\FloatTok{.2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = nox ~ bs(dis, df = 4), data = Boston)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.124622  -0.039259  -0.008514   0.020850   0.193891  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       0.73447    0.01460  50.306  < 2e-16 ***
## bs(dis, df = 4)1 -0.05810    0.02186  -2.658  0.00812 ** 
## bs(dis, df = 4)2 -0.46356    0.02366 -19.596  < 2e-16 ***
## bs(dis, df = 4)3 -0.19979    0.04311  -4.634 4.58e-06 ***
## bs(dis, df = 4)4 -0.38881    0.04551  -8.544  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.003837874)
## 
##     Null deviance: 6.7810  on 505  degrees of freedom
## Residual deviance: 1.9228  on 501  degrees of freedom
## AIC: -1371.9
## 
## Number of Fisher Scoring iterations: 2
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{plot}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis, }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{fits }\OtherTok{\textless{}{-}} \FunctionTok{predict}\NormalTok{(lm}\FloatTok{.2}\NormalTok{)}
\NormalTok{ord }\OtherTok{\textless{}{-}} \FunctionTok{order}\NormalTok{(Boston}\SpecialCharTok{$}\NormalTok{dis)}
\FunctionTok{lines}\NormalTok{(fits[ord] }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis[ord], }\AttributeTok{data=}\NormalTok{Boston)}
\end{Highlighting}
\end{Shaded}

\begin{center}\includegraphics{hw06-solutions_files/figure-latex/unnamed-chunk-5-1} \end{center}

According to help(bs,package=splines):

\begin{itemize}
\item
  The default degree of the spline is 3.
\item
  When we specify ``df=4'', the bs() function chooses df \(-\) degree
  \(=\) \(4-3=1\) internal knot, and also one knot each at the min and
  max of ``dis'', for a total of 3 knots. The internal knot is at the
  median of ``dis''.
\end{itemize}

This means we are fitting one cubic between min(dis) and median(dis),
and another cubic between median(dis) and max(dis).

\rm

\hypertarget{problem-1e}{%
\subsection{Problem 1(e)}\label{problem-1e}}

Now fit a regression spline for a range of degrees of freedom, and plot
the resulting fits and report the resulting RSS. Describe the results
obtained.

\em

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow=}\FunctionTok{c}\NormalTok{(}\DecValTok{4}\NormalTok{,}\DecValTok{3}\NormalTok{))}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{12}\NormalTok{) \{}
\NormalTok{  lm.fit }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\NormalTok{i), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{  rss }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(}\FunctionTok{residuals}\NormalTok{(lm.fit)}\SpecialCharTok{\^{}}\DecValTok{2}\NormalTok{)}
  \FunctionTok{plot}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis, }\AttributeTok{data=}\NormalTok{Boston, }\AttributeTok{col=}\StringTok{"grey"}\NormalTok{,}
       \AttributeTok{main=}\FunctionTok{paste0}\NormalTok{(}\StringTok{"Regression spline, df = "}\NormalTok{,i,}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{RSS = "}\NormalTok{,}\FunctionTok{round}\NormalTok{(rss,}\DecValTok{2}\NormalTok{)))}
\NormalTok{  fits }\OtherTok{\textless{}{-}} \FunctionTok{predict}\NormalTok{(lm.fit)}
\NormalTok{  ord }\OtherTok{\textless{}{-}} \FunctionTok{order}\NormalTok{(Boston}\SpecialCharTok{$}\NormalTok{dis)}
  \FunctionTok{lines}\NormalTok{(fits[ord] }\SpecialCharTok{\textasciitilde{}}\NormalTok{ dis[ord], }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{\}}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Warning in bs(dis, df = i): 'df' was too small; have used 3

## Warning in bs(dis, df = i): 'df' was too small; have used 3
\end{verbatim}

\begin{center}\includegraphics{hw06-solutions_files/figure-latex/unnamed-chunk-6-1} \end{center}

From the warnings above and what we learned about the bs() function in
part (d), when the formula ``df \(-\) degree'' says we should use 0
knots or a negative number of knots, the bs() function just sticks a
knot at the median again. Effectively, when we specify a df that is too
small, the bs() function just increases the df to the minimum
acceptable, which is df=3 in our case.

\rm

\hypertarget{problem-1f}{%
\subsection{Problem 1(f)}\label{problem-1f}}

Perform cross-validation or another approach in order to select the best
degrees of freedom for a regression spline on this data. Describe your
results.

\em

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# again, we can try LOOCV first...}
\NormalTok{cv.err }\OtherTok{\textless{}{-}} \ConstantTok{NULL}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{12}\NormalTok{) \{}
\NormalTok{  lm.fit }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\NormalTok{i), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{  cv.err }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(cv.err,}\FunctionTok{cv.glm}\NormalTok{(Boston,lm.fit)}\SpecialCharTok{$}\NormalTok{delta[}\DecValTok{1}\NormalTok{])  }
\NormalTok{\}}
\NormalTok{cv.err}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] 0.003874762 0.003874762 0.003874762 0.003893623 0.003704252 0.003704711
##  [7] 0.003711441 0.003699853 0.003731180 0.003692067 0.003715438 0.003701170
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"Optimal LOOCV df for bs() ="}\NormalTok{,}\FunctionTok{which}\NormalTok{(cv.err}\SpecialCharTok{==}\FunctionTok{min}\NormalTok{(cv.err)),}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Optimal LOOCV df for bs() = 10
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# ...and then 10{-}fold cv...}
\NormalTok{cv.err }\OtherTok{\textless{}{-}} \ConstantTok{NULL}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{12}\NormalTok{) \{}
\NormalTok{  lm.fit }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\NormalTok{i), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{  cv.err }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(cv.err,}\FunctionTok{cv.glm}\NormalTok{(Boston,lm.fit, }\AttributeTok{K=}\DecValTok{10}\NormalTok{)}\SpecialCharTok{$}\NormalTok{delta[}\DecValTok{1}\NormalTok{])  }
\NormalTok{\}}
\NormalTok{cv.err}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] 0.003891716 0.003879143 0.003886573 0.003913647 0.003687029 0.003698626
##  [7] 0.003684441 0.003694159 0.003737521 0.003690349 0.003712073 0.003709678
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"Optimal 10{-}fold CV df for bs() ="}\NormalTok{,}\FunctionTok{which}\NormalTok{(cv.err}\SpecialCharTok{==}\FunctionTok{min}\NormalTok{(cv.err)),}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Optimal 10-fold CV df for bs() = 7
\end{verbatim}

The 10-fold cross-validation result is somewhat unstable here and
depends on which random folds we get. In different runs, I've gotten
df=5, 6, 7, 8 and 10 as optimal. On the other hand, the LOOCV result is
stable and non-random (because it considers all \(n=506\) folds of size
\(n-1\), leaving one observation out at a time as a ``test dataset''),
and the result is df=10.

As a final check I will compare AIC and BIC for the 3rd degree
polynomial and several of the regression spline models:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lm.poly3 }\OtherTok{\textless{}{-}}\NormalTok{ lm}\FloatTok{.1}
\NormalTok{lm.bs3 }\OtherTok{\textless{}{-}}\NormalTok{ lm}\FloatTok{.2}
\NormalTok{lm.bs5 }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\DecValTok{5}\NormalTok{), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{lm.bs6 }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\DecValTok{6}\NormalTok{), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{lm.bs7 }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\DecValTok{7}\NormalTok{), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{lm.bs8 }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\DecValTok{8}\NormalTok{), }\AttributeTok{data=}\NormalTok{Boston)}
\NormalTok{lm.bs10 }\OtherTok{\textless{}{-}} \FunctionTok{glm}\NormalTok{(nox }\SpecialCharTok{\textasciitilde{}} \FunctionTok{bs}\NormalTok{(dis,}\AttributeTok{df=}\DecValTok{10}\NormalTok{), }\AttributeTok{data=}\NormalTok{Boston)}

\NormalTok{models }\OtherTok{\textless{}{-}} \FunctionTok{list}\NormalTok{(lm.poly3,lm.bs3,lm.bs5,lm.bs6,lm.bs7,lm.bs8,lm.bs10)}
\FunctionTok{names}\NormalTok{(models) }\OtherTok{\textless{}{-}} \FunctionTok{strsplit}\NormalTok{(}\StringTok{"lm.poly3,lm.bs3,lm.bs5,lm.bs6,lm.bs7,lm.bs8,lm.bs10"}\NormalTok{,}
                          \StringTok{","}\NormalTok{)[[}\DecValTok{1}\NormalTok{]]}

\NormalTok{comparison }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}
  \AttributeTok{logLik =} \FunctionTok{sapply}\NormalTok{(models,logLik),}
  \AttributeTok{df =} \FunctionTok{sapply}\NormalTok{(models, }\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{attributes}\NormalTok{(}\FunctionTok{logLik}\NormalTok{(x))}\SpecialCharTok{$}\NormalTok{df),}
  \AttributeTok{AIC =} \FunctionTok{sapply}\NormalTok{(models,AIC),}
  \AttributeTok{BIC =} \FunctionTok{sapply}\NormalTok{(models,BIC)}
\NormalTok{)}
\FunctionTok{rownames}\NormalTok{(comparison) }\OtherTok{\textless{}{-}} \FunctionTok{names}\NormalTok{(models)}

\FunctionTok{round}\NormalTok{(comparison,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##          logLik df      AIC      BIC
## lm.poly3 690.44  5 -1370.88 -1349.75
## lm.bs3   691.93  6 -1371.85 -1346.50
## lm.bs5   703.04  7 -1392.07 -1362.49
## lm.bs6   703.89  8 -1391.78 -1357.97
## lm.bs7   704.45  9 -1390.91 -1352.87
## lm.bs8   706.24 10 -1392.49 -1350.22
## lm.bs10  709.67 12 -1395.34 -1344.63
\end{verbatim}

We see that AIC likes the regression spline with df \(=\) 10 best, and
BIC likes the df \(=\) 5 regression spline model best. Referring back to
the regression spline plots, the curves for df \(=\) 5,6 or 7 seem
preferable to the somewhat wiggly df \(=\) 8 or 10 plots. (This is my
usual ``interpretability'' bias: smooth curves are generally easier to
tell a story about than wiggly ones. If your goal is prediction,
however, you may indeed prefer the df \(=\) 10 plot, which did better
with AIC and LOOCV.)

\rm

\hypertarget{problem-2-islr-10-pp.-324325.}{%
\section{Problem 2: ISLR \#10,
pp.~324--325.}\label{problem-2-islr-10-pp.-324325.}}

This question relates to the College data set.

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

Split the data into a training set and a test set. Using out-of-state
tuition as the response and the other variables as the predictors,
perform forward stepwise selection on the training set in order to
identify a satisfactory model that uses just a subset of the predictors.

\em

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{data}\NormalTok{(College)}
\FunctionTok{str}\NormalTok{(College)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 'data.frame':    777 obs. of  18 variables:
##  $ Private    : Factor w/ 2 levels "No","Yes": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Apps       : num  1660 2186 1428 417 193 ...
##  $ Accept     : num  1232 1924 1097 349 146 ...
##  $ Enroll     : num  721 512 336 137 55 158 103 489 227 172 ...
##  $ Top10perc  : num  23 16 22 60 16 38 17 37 30 21 ...
##  $ Top25perc  : num  52 29 50 89 44 62 45 68 63 44 ...
##  $ F.Undergrad: num  2885 2683 1036 510 249 ...
##  $ P.Undergrad: num  537 1227 99 63 869 ...
##  $ Outstate   : num  7440 12280 11250 12960 7560 ...
##  $ Room.Board : num  3300 6450 3750 5450 4120 ...
##  $ Books      : num  450 750 400 450 800 500 500 450 300 660 ...
##  $ Personal   : num  2200 1500 1165 875 1500 ...
##  $ PhD        : num  70 29 53 92 76 67 90 89 79 40 ...
##  $ Terminal   : num  78 30 66 97 72 73 93 100 84 41 ...
##  $ S.F.Ratio  : num  18.1 12.2 12.9 7.7 11.9 9.4 11.5 13.7 11.3 11.5 ...
##  $ perc.alumni: num  12 16 30 37 2 11 26 37 23 15 ...
##  $ Expend     : num  7041 10527 8735 19016 10922 ...
##  $ Grad.Rate  : num  60 56 54 59 15 55 63 73 80 52 ...
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# help(College)  {-}{-} to get the definitions of the variables}
\end{Highlighting}
\end{Shaded}

We have \(n=777<1000\) which is a smallish sample size for
cross-validation, so I will use a 50-50 split (see lecture 8 from
class).

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{n }\OtherTok{\textless{}{-}} \FunctionTok{dim}\NormalTok{(College)[}\DecValTok{1}\NormalTok{] }\DocumentationTok{\#\# should be n=777}
\NormalTok{train }\OtherTok{\textless{}{-}} \FunctionTok{sample}\NormalTok{(}\DecValTok{1}\SpecialCharTok{:}\NormalTok{n,}\AttributeTok{size=}\FunctionTok{floor}\NormalTok{(n}\SpecialCharTok{/}\DecValTok{2}\NormalTok{))}
\NormalTok{train }\OtherTok{\textless{}{-}} \FunctionTok{ifelse}\NormalTok{(}\DecValTok{1}\SpecialCharTok{:}\NormalTok{n }\SpecialCharTok{\%in\%}\NormalTok{ train,T,F)}
\NormalTok{test }\OtherTok{\textless{}{-}} \SpecialCharTok{!}\NormalTok{train}
\NormalTok{train.data }\OtherTok{\textless{}{-}}\NormalTok{ College[train,]}
\NormalTok{test.data }\OtherTok{\textless{}{-}}\NormalTok{ College[test,]}

\NormalTok{n.train }\OtherTok{\textless{}{-}} \FunctionTok{dim}\NormalTok{(train.data)[}\DecValTok{1}\NormalTok{]}
\NormalTok{n.test }\OtherTok{\textless{}{-}} \FunctionTok{dim}\NormalTok{(test.data)[}\DecValTok{1}\NormalTok{]}


\NormalTok{forward }\OtherTok{\textless{}{-}} \FunctionTok{regsubsets}\NormalTok{(Outstate }\SpecialCharTok{\textasciitilde{}}\NormalTok{ ., }\AttributeTok{data=}\NormalTok{train.data,}\AttributeTok{method=}\StringTok{"forward"}\NormalTok{,}\AttributeTok{nvmax=}\DecValTok{18}\NormalTok{)}

\NormalTok{tmp }\OtherTok{\textless{}{-}} \FunctionTok{summary}\NormalTok{(forward)}
\NormalTok{p.max }\OtherTok{\textless{}{-}} \FunctionTok{dim}\NormalTok{(train.data)[}\DecValTok{2}\NormalTok{]}
\NormalTok{p }\OtherTok{\textless{}{-}} \DecValTok{1}\SpecialCharTok{:}\NormalTok{(p.max}\DecValTok{{-}1}\NormalTok{)}

\NormalTok{results }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(tmp}\SpecialCharTok{$}\NormalTok{which,}\AttributeTok{rss=}\NormalTok{tmp}\SpecialCharTok{$}\NormalTok{rss,}\AttributeTok{adjr2=}\NormalTok{tmp}\SpecialCharTok{$}\NormalTok{adjr2,}
                    \AttributeTok{bic=}\NormalTok{n.train}\SpecialCharTok{*}\FunctionTok{log}\NormalTok{(tmp}\SpecialCharTok{$}\NormalTok{rss)}\SpecialCharTok{+}\FunctionTok{log}\NormalTok{(n.train)}\SpecialCharTok{*}\NormalTok{(p}\SpecialCharTok{+}\DecValTok{2}\NormalTok{),}
                    \AttributeTok{aic=}\NormalTok{n.train}\SpecialCharTok{*}\FunctionTok{log}\NormalTok{(tmp}\SpecialCharTok{$}\NormalTok{rss)}\SpecialCharTok{+}\DecValTok{2}\SpecialCharTok{*}\NormalTok{(p}\SpecialCharTok{+}\DecValTok{2}\NormalTok{),}
                    \AttributeTok{caic=}\NormalTok{n.train}\SpecialCharTok{*}\FunctionTok{log}\NormalTok{(tmp}\SpecialCharTok{$}\NormalTok{rss) }\SpecialCharTok{+} \DecValTok{2}\SpecialCharTok{*}\NormalTok{(p}\SpecialCharTok{+}\DecValTok{2}\NormalTok{) }\SpecialCharTok{+}
                      \DecValTok{2}\SpecialCharTok{*}\NormalTok{(p}\SpecialCharTok{+}\DecValTok{2}\NormalTok{)}\SpecialCharTok{*}\NormalTok{(p}\SpecialCharTok{+}\DecValTok{1}\NormalTok{)}\SpecialCharTok{/}\NormalTok{(n.train}\SpecialCharTok{{-}}\NormalTok{p}\DecValTok{{-}1}\NormalTok{))}

\DocumentationTok{\#\# results}

\FunctionTok{names}\NormalTok{(results)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] "X.Intercept." "PrivateYes"   "Apps"         "Accept"       "Enroll"      
##  [6] "Top10perc"    "Top25perc"    "F.Undergrad"  "P.Undergrad"  "Room.Board"  
## [11] "Books"        "Personal"     "PhD"          "Terminal"     "S.F.Ratio"   
## [16] "perc.alumni"  "Expend"       "Grad.Rate"    "rss"          "adjr2"       
## [21] "bic"          "aic"          "caic"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{names}\NormalTok{(results)[}\DecValTok{2}\NormalTok{] }\OtherTok{\textless{}{-}} \StringTok{"Private"}  

\NormalTok{minimize }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(res,col) \{}
\NormalTok{  obj }\OtherTok{\textless{}{-}}\NormalTok{ res[,col]}
\NormalTok{  k }\OtherTok{\textless{}{-}}\NormalTok{ (}\DecValTok{1}\SpecialCharTok{:}\FunctionTok{length}\NormalTok{(obj))[obj}\SpecialCharTok{==}\FunctionTok{min}\NormalTok{(obj)]}
  \FunctionTok{return}\NormalTok{(res[k,])}
\NormalTok{\}}


\NormalTok{sel.BIC }\OtherTok{\textless{}{-}} \FunctionTok{minimize}\NormalTok{(results,}\StringTok{"bic"}\NormalTok{)[}\DecValTok{2}\SpecialCharTok{:}\NormalTok{p.max]}
\NormalTok{sel.AIC }\OtherTok{\textless{}{-}} \FunctionTok{minimize}\NormalTok{(results,}\StringTok{"aic"}\NormalTok{)[}\DecValTok{2}\SpecialCharTok{:}\NormalTok{p.max]}

\NormalTok{fla.BIC }\OtherTok{\textless{}{-}} \FunctionTok{paste0}\NormalTok{(}\StringTok{"Outstate \textasciitilde{} "}\NormalTok{,}\FunctionTok{paste}\NormalTok{(}\FunctionTok{names}\NormalTok{(sel.BIC)[}\FunctionTok{unlist}\NormalTok{(sel.BIC)],}\AttributeTok{collapse=}\StringTok{" + "}\NormalTok{))}
\NormalTok{fla.AIC }\OtherTok{\textless{}{-}} \FunctionTok{paste0}\NormalTok{(}\StringTok{"Outstate \textasciitilde{} "}\NormalTok{,}\FunctionTok{paste}\NormalTok{(}\FunctionTok{names}\NormalTok{(sel.AIC)[}\FunctionTok{unlist}\NormalTok{(sel.AIC)],}\AttributeTok{collapse=}\StringTok{" + "}\NormalTok{))}

\NormalTok{fwd.BIC }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(fla.BIC, }\AttributeTok{data=}\NormalTok{train.data)}
\NormalTok{fwd.AIC }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(fla.AIC, }\AttributeTok{data=}\NormalTok{train.data)}
\end{Highlighting}
\end{Shaded}

Just for fun, I compare the best AIC and BIC models with a partial F
test, using the anova() function.

From the anova() output below we see that the AIC model adds the
variables Top25perc, Books and GradRate to the BIC model.

Unsurprisingly, AIC and the F test prefer the more complex model for
this split of the data. However, I will work with the simpler BiC model
in what follows.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{anova}\NormalTok{(fwd.BIC,fwd.AIC)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Variance Table
## 
## Model 1: Outstate ~ Private + Accept + Enroll + Room.Board + Personal + 
##     PhD + perc.alumni + Expend
## Model 2: Outstate ~ Private + Accept + Enroll + Top25perc + Room.Board + 
##     Books + Personal + PhD + perc.alumni + Expend + Grad.Rate
##   Res.Df        RSS Df Sum of Sq      F  Pr(>F)  
## 1    379 1553264546                              
## 2    376 1512317067  3  40947479 3.3935 0.01808 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

Note that the details of the results above depend strongly on the
particular split of the data into training and test subsets. A different
split, with a different set.seed() seed, could produce different (and
possibly non-nested) ``best'' AIC and BIC models; the result of the
anova() test (when the models turn out to be nested) can also change.

\rm

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

Fit a GAM on the training data, using out-of-state tuition as the
response and the features selected in the previous step as the
predictors. Plot the results, and explain your findings.

\em

Since ``Private'' is in the model, but ``Private'' is a 0/1 variable, it
doesn't get a transformation. But all the other variables are
quantitative and so get a transformation

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# Note that I am using cubic smoothing splines (bs=\textquotesingle{}cr\textquotesingle{} in the code below) }
\DocumentationTok{\#\# to be consistent with what we did in class, but you could equally well }
\DocumentationTok{\#\# use mgcv\textquotesingle{}s default "thin plate" splines...}

\NormalTok{fla.gam }\OtherTok{\textless{}{-}} \FunctionTok{paste0}\NormalTok{(}\StringTok{"Outstate \textasciitilde{} Private + s("}\NormalTok{,}
                  \FunctionTok{paste}\NormalTok{(}\FunctionTok{names}\NormalTok{(sel.BIC)[}\SpecialCharTok{{-}}\DecValTok{1}\NormalTok{][}\FunctionTok{unlist}\NormalTok{(sel.BIC)[}\SpecialCharTok{{-}}\DecValTok{1}\NormalTok{]],}
                  \AttributeTok{collapse=}\StringTok{",bs=\textquotesingle{}cr\textquotesingle{}) + s("}\NormalTok{),}\StringTok{",bs=\textquotesingle{}cr\textquotesingle{})"}\NormalTok{)}

\NormalTok{gam}\FloatTok{.1} \OtherTok{\textless{}{-}} \FunctionTok{gam}\NormalTok{(}\FunctionTok{formula}\NormalTok{(fla.gam), }\AttributeTok{data=}\NormalTok{train.data,}\AttributeTok{method=}\StringTok{"GCV.Cp"}\NormalTok{)}

\FunctionTok{summary}\NormalTok{(gam}\FloatTok{.1}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## Outstate ~ Private + s(Accept, bs = "cr") + s(Enroll, bs = "cr") + 
##     s(Room.Board, bs = "cr") + s(Personal, bs = "cr") + s(PhD, 
##     bs = "cr") + s(perc.alumni, bs = "cr") + s(Expend, bs = "cr")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8720.5      271.9  32.071  < 2e-16 ***
## PrivateYes    2618.3      346.8   7.549 3.65e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                  edf Ref.df      F  p-value    
## s(Accept)      3.666  4.275  8.276 2.70e-06 ***
## s(Enroll)      5.840  6.524  5.509 9.60e-06 ***
## s(Room.Board)  4.821  5.801  4.829 0.000168 ***
## s(Personal)    5.120  6.052  2.393 0.025545 *  
## s(PhD)         1.000  1.000 12.364 0.000494 ***
## s(perc.alumni) 2.150  2.733 10.397 8.82e-06 ***
## s(Expend)      4.244  4.925 32.353  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.796   Deviance explained = 81.1%
## GCV = 3.5935e+06  Scale est. = 3.3264e+06  n = 388
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{plot}\NormalTok{(gam}\FloatTok{.1}\NormalTok{,}\AttributeTok{pages=}\DecValTok{1}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{center}\includegraphics{hw06-solutions_files/figure-latex/unnamed-chunk-13-1} \end{center}

It appears that PhD and perc.alumni not need transformations, while
Accept and Expend clearly do. In particular it looks like

\begin{itemize}
\item
  Out of state tuition increases with institutional expenditures per
  student, up to a point at which it levels off (perhaps this is the
  highest tuition that the market will bear)
\item
  Out of state tuition also increases with the number of students
  accepted (which seems a bit suprising to me) but the influence of more
  students accepted, on out of state tuition, lessens after roughly
  Accept=4000.
\end{itemize}

\rm

\hypertarget{problem-2c}{%
\subsection{Problem 2(c)}\label{problem-2c}}

Evaluate the model obtained on the test set, and explain the results
obtained.

\em

First we do mean squared prediction error for the two models

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{mse.BIC }\OtherTok{\textless{}{-}} \FunctionTok{mean}\NormalTok{((test.data}\SpecialCharTok{$}\NormalTok{Outstate }\SpecialCharTok{{-}} \FunctionTok{predict}\NormalTok{(fwd.BIC,}\AttributeTok{newdata=}\NormalTok{test.data))}\SpecialCharTok{\^{}}\DecValTok{2}\NormalTok{)}
\NormalTok{mse.gam }\OtherTok{\textless{}{-}} \FunctionTok{mean}\NormalTok{((test.data}\SpecialCharTok{$}\NormalTok{Outstate }\SpecialCharTok{{-}} \FunctionTok{predict}\NormalTok{(gam}\FloatTok{.1}\NormalTok{,}\AttributeTok{newdata=}\NormalTok{test.data))}\SpecialCharTok{\^{}}\DecValTok{2}\NormalTok{)}

\NormalTok{mse.results }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\StringTok{"Prediction MSE"}\OtherTok{=}\FunctionTok{c}\NormalTok{(}\AttributeTok{fwd.BIC=}\NormalTok{mse.BIC,}\AttributeTok{gam.1=}\NormalTok{mse.gam))}
\FunctionTok{round}\NormalTok{(mse.results)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##         Prediction.MSE
## fwd.BIC        4524190
## gam.1          3520081
\end{verbatim}

Now we do a goodness-of-fit comparison:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{test.BIC }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(}\FunctionTok{formula}\NormalTok{(fla.BIC),}\AttributeTok{data=}\NormalTok{test.data)}
\NormalTok{test.gam }\OtherTok{\textless{}{-}} \FunctionTok{gam}\NormalTok{(}\FunctionTok{formula}\NormalTok{(fla.gam),}\AttributeTok{data=}\NormalTok{test.data)}
\NormalTok{models }\OtherTok{\textless{}{-}} \FunctionTok{list}\NormalTok{(test.BIC,test.gam)}
\FunctionTok{names}\NormalTok{(models) }\OtherTok{\textless{}{-}} \FunctionTok{strsplit}\NormalTok{(}\StringTok{"test.BIC,test.gam"}\NormalTok{,}\StringTok{","}\NormalTok{)[[}\DecValTok{1}\NormalTok{]]}

\NormalTok{comparison }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}
  \AttributeTok{logLik =} \FunctionTok{sapply}\NormalTok{(models,logLik),}
  \AttributeTok{df =} \FunctionTok{sapply}\NormalTok{(models, }\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{attributes}\NormalTok{(}\FunctionTok{logLik}\NormalTok{(x))}\SpecialCharTok{$}\NormalTok{df),}
  \AttributeTok{AIC =} \FunctionTok{sapply}\NormalTok{(models,AIC),}
  \AttributeTok{BIC =} \FunctionTok{sapply}\NormalTok{(models,BIC)}
\NormalTok{)}
\FunctionTok{rownames}\NormalTok{(comparison) }\OtherTok{\textless{}{-}} \FunctionTok{names}\NormalTok{(models)}

\FunctionTok{round}\NormalTok{(comparison,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##            logLik    df     AIC     BIC
## test.BIC -3502.34 10.00 7024.68 7064.32
## test.gam -3442.53 20.45 6925.96 7007.01
\end{verbatim}

All three comparisons, prediction MSE, AIC, and BIC, prefer the gam
model in this example.

Explanation: The extra flexibility of the nonparametric transformations
of the predictors provide better predictions, reducing the mean squared
prediction error by
(\ensuremath{3.5200809\times 10^{6}})/(\ensuremath{4.5241897\times 10^{6}})
= 0.78, as well as dramatic improvements in fit as reflected by
reductions of 60--100 in AIC and BIC.

\rm

\hypertarget{problem-2d}{%
\subsection{Problem 2(d)}\label{problem-2d}}

For which variables, if any, is there evidence of a non-linear
relationship with the response?

\em

From the transformation plots above, it seems clear that Accept and
Expand have a non-linear relationship with the response. There really
isn't a clear case for any of the others, though one could make an
argument for Enroll and Personal.

One could more formally assess whether a nonlinear term is needed for
Personal, for example, by comparing the model with s(Personal) as a
predictor, with the model that just uses Personal as a predictor.

Once again, these results may vary, depending on what split into
training and test data sets you are working with.

\rm

\hypertarget{problem-3-based-on-gelman-hill-chapter-9-4.}{%
\section{Problem 3 (Based on Gelman \& Hill, Chapter 9,
\#4).}\label{problem-3-based-on-gelman-hill-chapter-9-4.}}

The table below describes a hypothetical experiment on 2400 persons.
Each row of the table specifies a category of person, as defined by his
or her pre-treatment predictor \(x\), treatment indicator \(T\), and
potential outcomes \(y^0\) and \(y^1\). (For simplicity, we
assume---unrealistically---that all people in this experiment fit into
one of these eight categories).

\begin{center}
\begin{tabular}{cccccc}
Category & \# persons in category & $x$ & $T$ & $y^0$ & $y^1$ \\
\hline
\hline
1 & 300 & 0 & 0 & 4 & 6 \\
2 & 300 & 1 & 0 & 4 & 6 \\
3 & 500 & 0 & 1 & 4 & 6 \\
4 & 500 & 1 & 1 & 4 & 6 \\
5 & 200 & 0 & 0 & 10 & 12 \\
6 & 200 & 1 & 0 & 10 & 12 \\
7 & 200 & 0 & 1 & 10 & 12 \\
8 & 200 & 1 & 1 & 10 & 12 \\
\hline
\end{tabular}
\end{center}

In making this table we are assuming omniscience, so that we know both
\(y^0\) and \(y^1\) for all observations. But the (non-omniscient)
investigator woud only observe \(x\), \(T\), and \(y=y^T\) for each
unit. (For example, a person in category 1 would have \(x=0\), \(T=0\)
and \(y=4\), and a person in category 3 would have \(x=0\), \(T=1\), and
\(y=6\).)

\em

Let's start by putting the data above into R\ldots{}

\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}

\rm

\hypertarget{problem-3a}{%
\subsection{Problem 3(a)}\label{problem-3a}}

What is the true ACE (average causal effect), if we could observe
\(y^1\) and \(y^0\) for every person in this population of 2400 persons?

\em

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

\rm

\hypertarget{problem-3b}{%
\subsection{Problem 3(b)}\label{problem-3b}}

Another population quantity is the mean of \(y\) for those who received
the treatment, minus the mean of \(y\) for those who did not. What is
the relationship between this quantity and the quantity you calculated
in part (a)?

\em

\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! Below in part (d) we explore the
data further and discover that there is another confounding variable not
accounted for in the data.

\rm

\hypertarget{problem-3c}{%
\subsection{Problem 3(c)}\label{problem-3c}}

Suppose we draw a person randomly from this population or 2400 people.
What is the probability that \(T=1\) for this person? If I tell you that
\(x=1\) for this person, does that change the probability that \(T=1\)?
I.e., is \(P[T=1|x=1]=P[T=1]\)?

\em

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

\rm

\hypertarget{problem-3d}{%
\subsection{Problem 3(d)}\label{problem-3d}}

Is it plausible to believe this data came from a randomized experiment?
Defend your answer.

\em

From part (c) we know that \(P[T=1]=P[T=1|x=1]\), and more generally,
\(P[T=t,x=x] = P[T=t]P[x=x]\), i.e.~\(T\) and \(X\) are independent. So,
if we think \(x\) is the only pre-treatment covariate (as suggested in
the beginning of the problem), it is plausible that this is data from a
randomized experiemnt (since randomization forces \(T\) to be indpendent
of any pre-treatment covariates).

However, looking back at the data table, it seems that there is a jump
of 6 in the level of response (for both \(y^0\) and \(y^1\)) going from
the first four categories, vs.~the last four categories. For example, if
\(y\) is measuring a health outcome after a hospital stay, and the first
four categories represent people who would go to hospital \#1 and the
second four categories represent people who would go to hospital \#2, it
appears that the first group of people are generally less healthy than
the second group.

If we define another grouping variable

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df}\SpecialCharTok{$}\NormalTok{x2 }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

then (i) we can see that \texttt{x2} and \texttt{T} are not independent:
\texttt{P{[}T=1\textbar{}x2=1{]}\ =\ 1000/1600\ =} 0.625 \(\neq\)
\texttt{P{[}T=1\textbar{}x2=0{]}\ =\ 400/800\ =} 0.5, and (ii) in the
regression \texttt{y\ \textasciitilde{}\ T\ +\ x2}, the coefficient
\(\hat\beta_T\) would plausibly be an unbiased estimate of ACE. The
simpler regression \texttt{y\ \textasciitilde{}\ T} (which is
essentially what we were doing in part (b) above) would produce a biased
\(\hat\beta_T\).

\rm

\hypertarget{problem-3e}{%
\subsection{Problem 3(e)}\label{problem-3e}}

For this hypothetical data, figure out (by hand, or by using R) the
estimate and standard error of the coefficient of \(T\) in a regression
of \(y\) on \(T\) and \(x\). Is this a reasonable way to estimate the
true ACE?

\em

Below are two 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 (b). 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{splitting-the-data-into-individual-cases}{%
\subsubsection{Splitting the data into individual
cases}\label{splitting-the-data-into-individual-cases}}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{new.df }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{,}\DecValTok{0}\NormalTok{)  }\DocumentationTok{\#\# a dummy line so we can rbind to it}
\FunctionTok{names}\NormalTok{(new.df) }\OtherTok{\textless{}{-}} \FunctionTok{names}\NormalTok{(df)}
\ControlFlowTok{for}\NormalTok{ (i }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\DecValTok{8}\NormalTok{) \{}
  \ControlFlowTok{for}\NormalTok{ (j }\ControlFlowTok{in} \DecValTok{1}\SpecialCharTok{:}\NormalTok{df[i,}\StringTok{"num"}\NormalTok{]) \{}
\NormalTok{    new.df }\OtherTok{\textless{}{-}} \FunctionTok{rbind}\NormalTok{(new.df,df[i,])}
\NormalTok{  \}}
\NormalTok{\}}
\NormalTok{new.df }\OtherTok{\textless{}{-}}\NormalTok{ new.df[}\SpecialCharTok{{-}}\DecValTok{1}\NormalTok{,]  }\DocumentationTok{\#\# remove the dummy line}

\NormalTok{new.df}\SpecialCharTok{$}\NormalTok{y }\OtherTok{\textless{}{-}} \FunctionTok{ifelse}\NormalTok{(new.df}\SpecialCharTok{$}\NormalTok{T}\SpecialCharTok{==}\DecValTok{0}\NormalTok{,new.df}\SpecialCharTok{$}\NormalTok{y\_0,new.df}\SpecialCharTok{$}\NormalTok{y\_1)}

\FunctionTok{summary}\NormalTok{(lm.splitting }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(y }\SpecialCharTok{\textasciitilde{}}\NormalTok{ T }\SpecialCharTok{+}\NormalTok{ x,}\AttributeTok{data=}\NormalTok{new.df))}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = y ~ T + x, data = new.df)
## 
## 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}

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

A weighted regression solution is faster than expanding out the full
dataset,

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df}\SpecialCharTok{$}\NormalTok{y }\OtherTok{\textless{}{-}} \FunctionTok{ifelse}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{T}\SpecialCharTok{==}\DecValTok{0}\NormalTok{,df}\SpecialCharTok{$}\NormalTok{y\_0,df}\SpecialCharTok{$}\NormalTok{y\_1)}

\NormalTok{lm.weighted }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(y }\SpecialCharTok{\textasciitilde{}}\NormalTok{ T }\SpecialCharTok{+}\NormalTok{ x, }\AttributeTok{data =}\NormalTok{ df, }\AttributeTok{weights =}\NormalTok{ num)}
\FunctionTok{summary}\NormalTok{(lm.weighted)  }\DocumentationTok{\#\# uncorrected summary}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = y ~ T + x, data = df, weights = num)
## 
## Weighted Residuals:
##      1      2      3      4      5      6      7      8 
## -41.57 -41.57 -38.33 -38.33  50.91  50.91  60.61  60.61 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    6.400      2.316   2.764   0.0397 *
## T              1.314      2.547   0.516   0.6279  
## x              0.000      2.512   0.000   1.0000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 61.52 on 5 degrees of freedom
## Multiple R-squared:  0.05055,    Adjusted R-squared:  -0.3292 
## F-statistic: 0.1331 on 2 and 5 DF,  p-value: 0.8784
\end{verbatim}

but requires a correction for the model's residual degrees of freedom.
See this
\href{https://www.r-bloggers.com/2015/09/linear-models-with-weighted-observations/}{blogpost}
for more comments.

The correct residual degrees of freedom are

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{df.res }\OtherTok{\textless{}{-}} \FunctionTok{sum}\NormalTok{(df}\SpecialCharTok{$}\NormalTok{num) }\SpecialCharTok{{-}} \FunctionTok{length}\NormalTok{(lm.weighted}\SpecialCharTok{$}\NormalTok{coefficients)}
\NormalTok{df.res}
\end{Highlighting}
\end{Shaded}

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

and then new \(s^2\) and \(SE(\hat\beta_T)\) can be calculated with this
new residual df.

For now we will just observe that the point estimates are the same: Note
that we get the same estimate either way:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{rbind}\NormalTok{(}\StringTok{"b\_T from lm.splitting"}\OtherTok{=}\FunctionTok{summary}\NormalTok{(lm.splitting)}\SpecialCharTok{$}\NormalTok{coef[}\DecValTok{2}\NormalTok{,}\DecValTok{1}\NormalTok{],}
      \StringTok{"b\_T from lm.weighted"}\OtherTok{=}\FunctionTok{summary}\NormalTok{(lm.weighted)}\SpecialCharTok{$}\NormalTok{coef[}\DecValTok{2}\NormalTok{,}\DecValTok{1}\NormalTok{]),}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                       [,1]
## b_T from lm.splitting 1.31
## b_T from lm.weighted  1.31
\end{verbatim}

This is \emph{not} a reasonable way to estimate the ACE: Since we are
not accounting for the variable \texttt{x2}, we get the same biased
estimate as in part (b).

If we include \texttt{x2} also, we can get a better answer:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{corrected.lm.splitting }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(y }\SpecialCharTok{\textasciitilde{}}\NormalTok{ T }\SpecialCharTok{+}\NormalTok{ x }\SpecialCharTok{+}\NormalTok{ x2,}\AttributeTok{data=}\NormalTok{new.df)}
\FunctionTok{summary}\NormalTok{(corrected.lm.splitting)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = y ~ T + x + x2, data = new.df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -3.731e-12 -9.000e-16  2.900e-15  3.900e-15  1.803e-12 
## 
## Coefficients:
##               Estimate Std. Error    t value Pr(>|t|)    
## (Intercept)  1.000e+01  4.833e-15  2.069e+15  < 2e-16 ***
## T            2.000e+00  4.397e-15  4.549e+14  < 2e-16 ***
## x           -1.535e-14  4.304e-15 -3.566e+00  0.00037 ***
## x2          -6.000e+00  4.598e-15 -1.305e+15  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.054e-13 on 2396 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 5.977e+29 on 3 and 2396 DF,  p-value: < 2.2e-16
\end{verbatim}

We can see that the coefficient \(\hat\beta_T\) is now the correct
value, 2, to many decimal places (and with a very tight confidence
interval!).

\rm

\hypertarget{problem-4-based-on-gelman-hill-chapter-9-6.}{%
\section{Problem 4 (Based on Gelman \& Hill, Chapter 9,
\#6).}\label{problem-4-based-on-gelman-hill-chapter-9-6.}}

You are consulting for a researcher who has performed a randomized trial
where the treatment was a series of 26 weekly therapy sessions, the
control was no therapy, and the outcome was self-report of emotional
state one year later. However, most people in the treatment group did
not attend every therapy session. In fact there was a good deal of
variation in the number of therapy sessions actually attended. The
researcher is concerned that her results represent ``watered down'\,'
estimates because of this variation and suggests adding in another
predictor to the model: number of therapy sessions attended.

\hypertarget{problem-4a}{%
\subsection{Problem 4(a)}\label{problem-4a}}

What would you advise her about her suggestion? Carefully justify your
answer in terms of our discussion about controlled experiments,
observational studies, confounders, and so forth.

\em

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.

\rm

\hypertarget{problem-4b}{%
\subsection{Problem 4(b)}\label{problem-4b}}

Can you suggest another kind of analysis that might give her a better
estimate of the treatment effect? Carefully explain why this might work.
\em (Hint: Think about the more sophisticated analyses presented in lecture 16 in class.)

\rm

\em

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.

\rm

\end{document}
