% 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 06 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}{5}
\usepackage{float}
\floatplacement{figure}{H}
\usepackage{booktabs}
\usepackage{longtable}
\usepackage{array}
\usepackage{multirow}
\usepackage{wrapfig}
\usepackage{float}
\usepackage{colortbl}
\usepackage{pdflscape}
\usepackage{tabu}
\usepackage{threeparttable}
\usepackage{threeparttablex}
\usepackage[normalem]{ulem}
\usepackage{makecell}
\usepackage{xcolor}
\ifLuaTeX
  \usepackage{selnolig}  % disable illegal ligatures
\fi

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

\begin{document}
\maketitle

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

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(tidyverse)}
\FunctionTok{library}\NormalTok{(kableExtra)}
\FunctionTok{library}\NormalTok{(GGally)}
\FunctionTok{library}\NormalTok{(grid)}
\FunctionTok{library}\NormalTok{(gridExtra)}
\FunctionTok{library}\NormalTok{(ggplotify)}
\FunctionTok{library}\NormalTok{(reshape2)}

\NormalTok{cdidata }\OtherTok{\textless{}{-}} \FunctionTok{read.table}\NormalTok{(}\StringTok{"cdi.dat"}\NormalTok{,}\AttributeTok{header=}\NormalTok{T)}
\end{Highlighting}
\end{Shaded}

\textbf{Note:} There is more detail in the solutions for problem \#1 than I would expect to see in a technical appendix. I have included extra detail so that you can see some of my thinking, and some of my approaches to problem-solving in applied regression. When you begin to convert your revised solutions to problem \#1 into a technical appendix for the project 1 IDMRAD paper, you should include technical material (and explanations) that help support the claims that you will be making in the main paper. You know you have done the technical appendix right, if you can tell the reader after each claim you make, which page(s) of the technical appendix has more evidence, calculations, etc., to back up that claim. If you do not refer to part of your technical appendix in the main body of your paper, you may not need to include that part in your technical appendix at all.

\hypertarget{part-a}{%
\subsection{Part a}\label{part-a}}

\hypertarget{make-a-table-or-tables-showing-appropriate-summary-statistics-for-each-variable-in-the-data-set.-note-that-summary-statistics-for-continuous-variables-will-be-different-from-the-summary-statistics-for-categorical-variables.}{%
\subsubsection{Make a table or tables showing appropriate summary statistics for each variable in the data set. Note that summary statistics for continuous variables will be different from the summary statistics for categorical variables.}\label{make-a-table-or-tables-showing-appropriate-summary-statistics-for-each-variable-in-the-data-set.-note-that-summary-statistics-for-continuous-variables-will-be-different-from-the-summary-statistics-for-categorical-variables.}}

First I'll just look at the ``head'' of all the variables (see Tables 1 and 2); this is a bit like looking at \texttt{str()}, which is also worth looking at (I split the ``head'' into two parts, because there are too many variables to fit horizontally on the page).

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{head}\NormalTok{(cdidata[,}\DecValTok{1}\SpecialCharTok{:}\DecValTok{10}\NormalTok{]) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{caption=}\StringTok{" "}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{table}

\caption{\label{tab:unnamed-chunk-3} }
\centering
\begin{tabular}[t]{rllrrrrrrr}
\toprule
id & county & state & land.area & pop & pop.18\_34 & pop.65\_plus & doctors & hosp.beds & crimes\\
\midrule
1 & Los\_Angeles & CA & 4060 & 8863164 & 32.1 & 9.7 & 23677 & 27700 & 688936\\
2 & Cook & IL & 946 & 5105067 & 29.2 & 12.4 & 15153 & 21550 & 436936\\
3 & Harris & TX & 1729 & 2818199 & 31.3 & 7.1 & 7553 & 12449 & 253526\\
4 & San\_Diego & CA & 4205 & 2498016 & 33.5 & 10.9 & 5905 & 6179 & 173821\\
5 & Orange & CA & 790 & 2410556 & 32.6 & 9.2 & 6062 & 6369 & 144524\\
\addlinespace
6 & Kings & NY & 71 & 2300664 & 28.3 & 12.4 & 4861 & 8942 & 680966\\
\bottomrule
\end{tabular}
\end{table}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{head}\NormalTok{(cdidata[,}\FunctionTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{,}\DecValTok{11}\SpecialCharTok{:}\DecValTok{17}\NormalTok{)]) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{caption=}\StringTok{" "}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{table}

\caption{\label{tab:unnamed-chunk-3} }
\centering
\begin{tabular}[t]{rrrrrrrl}
\toprule
id & pct.hs.grad & pct.bach.deg & pct.below.pov & pct.unemp & per.cap.income & tot.income & region\\
\midrule
1 & 70.0 & 22.3 & 11.6 & 8.0 & 20786 & 184230 & W\\
2 & 73.4 & 22.8 & 11.1 & 7.2 & 21729 & 110928 & NC\\
3 & 74.9 & 25.4 & 12.5 & 5.7 & 19517 & 55003 & S\\
4 & 81.9 & 25.3 & 8.1 & 6.1 & 19588 & 48931 & W\\
5 & 81.2 & 27.8 & 5.2 & 4.8 & 24400 & 58818 & W\\
\addlinespace
6 & 63.7 & 16.6 & 19.5 & 9.5 & 16803 & 38658 & NE\\
\bottomrule
\end{tabular}
\end{table}

We can also check to see how many unique values each variable has (Table 3; this is especially relevant for the categorical variables)

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{apply}\NormalTok{(cdidata,}\DecValTok{2}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) \{}\FunctionTok{length}\NormalTok{(}\FunctionTok{unique}\NormalTok{(x))\}) }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{col.names=}\StringTok{"unique values"}\NormalTok{,}\AttributeTok{caption=}\StringTok{" "}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{kable\_classic}\NormalTok{(}\AttributeTok{full\_width=}\NormalTok{F)}
\end{Highlighting}
\end{Shaded}

\begin{table}

\caption{\label{tab:unnamed-chunk-4} }
\centering
\begin{tabular}[t]{lr}
\toprule
  & unique values\\
\midrule
id & 440\\
county & 373\\
state & 48\\
land.area & 384\\
pop & 440\\
\addlinespace
pop.18\_34 & 149\\
pop.65\_plus & 137\\
doctors & 360\\
hosp.beds & 391\\
crimes & 437\\
\addlinespace
pct.hs.grad & 223\\
pct.bach.deg & 220\\
pct.below.pov & 155\\
pct.unemp & 97\\
per.cap.income & 436\\
\addlinespace
tot.income & 428\\
region & 4\\
\bottomrule
\end{tabular}
\end{table}

It looks like \texttt{id} is just the same as the row number for each row of the data frame -- not useful for data analysis!

\texttt{state} has 48 values -- that is a lot, and so I may set this variable aside until later in my analyses

\texttt{county} is a categorical variable with nearly as many unique values (373) as rows in the cdidata data frame (440). A little more exploration (see R code just below) shows that if I combine \texttt{county} with \texttt{state}, I get 440 unique values: some counties in different states have the same name. The upshot is that we only have one observation per \emph{unique} county (what we might have thought were multiple observations per county are really single observations from counties with the same name in different states), and so \texttt{county} is not a useful variable to include in models (though it might be useful to identify outliers and other interesting data points in our analyses)

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{county.state }\OtherTok{\textless{}{-}} \FunctionTok{with}\NormalTok{(cdidata,}\FunctionTok{paste}\NormalTok{(county,state))}
\NormalTok{tmp }\OtherTok{\textless{}{-}} \FunctionTok{as.data.frame}\NormalTok{(}\FunctionTok{matrix}\NormalTok{(}\FunctionTok{sort}\NormalTok{(county.state),}\AttributeTok{ncol=}\DecValTok{4}\NormalTok{))}
\FunctionTok{names}\NormalTok{(tmp) }\OtherTok{\textless{}{-}} \FunctionTok{paste}\NormalTok{(}\StringTok{"Counties"}\NormalTok{,}\FunctionTok{c}\NormalTok{(}\StringTok{"1{-}110"}\NormalTok{,}\StringTok{"111{-}220"}\NormalTok{,}\StringTok{"221{-}330"}\NormalTok{,}\StringTok{"331{-}440"}\NormalTok{))}
\NormalTok{tmp[}\DecValTok{1}\SpecialCharTok{:}\DecValTok{30}\NormalTok{,] }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{longtable=}\NormalTok{T,}\AttributeTok{caption=}\StringTok{" "}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{(}\AttributeTok{full\_width=}\NormalTok{F)}
\end{Highlighting}
\end{Shaded}

\begin{longtable}[t]{llll}
\caption{\label{tab:unnamed-chunk-5} }\\
\toprule
Counties 1-110 & Counties 111-220 & Counties 221-330 & Counties 331-440\\
\midrule
Ada ID & Ector TX & Lycoming PA & Rockingham NH\\
Adams CO & El\_Dorado CA & Macomb MI & Rockland NY\\
Aiken SC & El\_Paso CO & Macon IL & Rowan NC\\
Alachua FL & El\_Paso TX & Madison AL & Rutherford TN\\
Alamance NC & Elkhart IN & Madison IL & Sacramento CA\\
\addlinespace
Alameda CA & Erie NY & Madison IN & Saginaw MI\\
Albany NY & Erie PA & Mahoning OH & Salt\_Lake UT\\
Alexandria\_City VA & Escambia FL & Manatee FL & San\_Bernardino CA\\
Allegheny PA & Essex MA & Marathon WI & San\_Diego CA\\
Allen IN & Essex NJ & Maricopa AZ & San\_Francisco CA\\
\addlinespace
Allen OH & Fairfax\_County VA & Marin CA & San\_Joaquin CA\\
Anderson SC & Fairfield CT & Marion FL & San\_Luis\_Obispo CA\\
Androscoggin ME & Fairfield OH & Marion IN & San\_Mateo CA\\
Anne\_Arundel MD & Fayette KY & Marion OR & Sangamon IL\\
Arapahoe CO & Fayette PA & Martin FL & Santa\_Barbara CA\\
\addlinespace
Arlington\_County VA & Florence SC & Maui HI & Santa\_Clara CA\\
Atlantic NJ & Forsyth NC & McHenry IL & Santa\_Cruz CA\\
Baltimore MD & Fort\_Bend TX & McLean IL & Sarasota FL\\
Baltimore\_City MD & Franklin OH & McLennan TX & Saratoga NY\\
Barnstable MA & Franklin PA & Mecklenburg NC & Sarpy NE\\
\addlinespace
Bay FL & Frederick MD & Medina OH & Schenectady NY\\
Bay MI & Fresno CA & Merced CA & Schuylkill PA\\
Beaver PA & Fulton GA & Mercer NJ & Sedgwick KS\\
Bell TX & Galveston TX & Mercer PA & Seminole FL\\
Benton WA & Gaston NC & Merrimack NH & Shasta CA\\
\addlinespace
Bergen NJ & Genesee MI & Middlesex CT & Shawnee KS\\
Berks PA & Gloucester NJ & Middlesex MA & Sheboygan WI\\
Berkshire MA & Greene MO & Middlesex NJ & Shelby TN\\
Bernalillo NM & Greene OH & Midland TX & Smith TX\\
Berrien MI & Greenville SC & Milwaukee WI & Snohomish WA\\
\bottomrule
\end{longtable}

(I have only listed the first 30 rows, to give a sense of the list of county+state combinations; you can delete the \texttt{{[}1:30,{]}} in the R code to see the whole table\ldots)

Scanning through this table of county \& state names, we see that Allen is the first county name that appears in multiple states. And a little further down, both ``Baltimore MD'' and ``Baltimore\_City MD'' are listed; this makes me wonder whether these two data points are really independent (Baltimore the city is inside Baltimore the county!), but I will ignore possible dependence like this in my analyses below.

Next let's make a table with the usual one-dimensional summary statistics (Table 5)

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cdinumeric }\OtherTok{\textless{}{-}}\NormalTok{ cdidata[,}\SpecialCharTok{{-}}\FunctionTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{,}\DecValTok{2}\NormalTok{,}\DecValTok{3}\NormalTok{,}\DecValTok{17}\NormalTok{)] }\DocumentationTok{\#\# get rid of id, county, state and (for now) region}
\FunctionTok{apply}\NormalTok{(cdinumeric,}\DecValTok{2}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{c}\NormalTok{(}\FunctionTok{summary}\NormalTok{(x),}\AttributeTok{SD=}\FunctionTok{sd}\NormalTok{(x))) }\SpecialCharTok{\%\textgreater{}\%}\NormalTok{ as.data.frame }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{t}\NormalTok{() }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{round}\NormalTok{(}\AttributeTok{digits=}\DecValTok{2}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{caption=}\StringTok{" "}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{table}

\caption{\label{tab:unnamed-chunk-6} }
\centering
\begin{tabular}[t]{lrrrrrrr}
\toprule
  & Min. & 1st Qu. & Median & Mean & 3rd Qu. & Max. & SD\\
\midrule
land.area & 15.0 & 451.25 & 656.50 & 1041.41 & 946.75 & 20062.0 & 1549.92\\
pop & 100043.0 & 139027.25 & 217280.50 & 393010.92 & 436064.50 & 8863164.0 & 601987.02\\
pop.18\_34 & 16.4 & 26.20 & 28.10 & 28.57 & 30.02 & 49.7 & 4.19\\
pop.65\_plus & 3.0 & 9.88 & 11.75 & 12.17 & 13.62 & 33.8 & 3.99\\
doctors & 39.0 & 182.75 & 401.00 & 988.00 & 1036.00 & 23677.0 & 1789.75\\
\addlinespace
hosp.beds & 92.0 & 390.75 & 755.00 & 1458.63 & 1575.75 & 27700.0 & 2289.13\\
crimes & 563.0 & 6219.50 & 11820.50 & 27111.62 & 26279.50 & 688936.0 & 58237.51\\
pct.hs.grad & 46.6 & 73.88 & 77.70 & 77.56 & 82.40 & 92.9 & 7.02\\
pct.bach.deg & 8.1 & 15.28 & 19.70 & 21.08 & 25.33 & 52.3 & 7.65\\
pct.below.pov & 1.4 & 5.30 & 7.90 & 8.72 & 10.90 & 36.3 & 4.66\\
\addlinespace
pct.unemp & 2.2 & 5.10 & 6.20 & 6.60 & 7.50 & 21.3 & 2.34\\
per.cap.income & 8899.0 & 16118.25 & 17759.00 & 18561.48 & 20270.00 & 37541.0 & 4059.19\\
tot.income & 1141.0 & 2311.00 & 3857.00 & 7869.27 & 8654.25 & 184230.0 & 12884.32\\
\bottomrule
\end{tabular}
\end{table}

and a separate table for the one categorical variable that I'll be using a lot (Table 6)

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tmp }\OtherTok{\textless{}{-}} \FunctionTok{rbind}\NormalTok{(}\FunctionTok{with}\NormalTok{(cdidata,}\FunctionTok{table}\NormalTok{(region)))}
\FunctionTok{row.names}\NormalTok{(tmp) }\OtherTok{\textless{}{-}} \StringTok{"Freq"}
\NormalTok{tmp }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{caption=}\StringTok{" "}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{(}\AttributeTok{full\_width=}\NormalTok{F)}
\end{Highlighting}
\end{Shaded}

\begin{table}

\caption{\label{tab:unnamed-chunk-7} }
\centering
\begin{tabular}[t]{lrrrr}
\toprule
  & NC & NE & S & W\\
\midrule
Freq & 108 & 103 & 152 & 77\\
\bottomrule
\end{tabular}
\end{table}

I don't see anything really exciting in the continuous variables (Table 5). There are several variables with Mean substantially larger than Median (\texttt{land.area}, \texttt{pop}, \texttt{doctors}, \texttt{hosp.beds}, \texttt{crimes}, \texttt{per.cap.income}, and \texttt{total.income}), indicating possible right-skewing. There are no variables with Mean substantially smaller than Median.

For the \texttt{region} variable (Table 6), it might be of some interest that the most counties are in the South (region `S') and the least are in the west (region `W'). The low number of counties in the West could be indicative of a lack of sampling in the West, or it could be that counties are just larger (in land area) in the West, so there are fewer counties to sample from. Similarly, the high number of counties in the South could be indicative of over-sampling, or perhap the South simply has a lot of counties that cover only small land areas.

\clearpage

\hypertarget{indicate-where-in-which-variables-there-is-missing-data-nas-if-any-how-much-there-is-in-each-variable-and-why-it-might-be-there.}{%
\subsubsection{Indicate where (in which variables) there is missing data (NA's), if any, how much there is (in each variable) and why it might be there.}\label{indicate-where-in-which-variables-there-is-missing-data-nas-if-any-how-much-there-is-in-each-variable-and-why-it-might-be-there.}}

We can check for NA's directly:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{apply}\NormalTok{(cdidata,}\DecValTok{2}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{any}\NormalTok{(}\FunctionTok{is.na}\NormalTok{(x)) )}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##             id         county          state      land.area            pop 
##          FALSE          FALSE          FALSE          FALSE          FALSE 
##      pop.18_34    pop.65_plus        doctors      hosp.beds         crimes 
##          FALSE          FALSE          FALSE          FALSE          FALSE 
##    pct.hs.grad   pct.bach.deg  pct.below.pov      pct.unemp per.cap.income 
##          FALSE          FALSE          FALSE          FALSE          FALSE 
##     tot.income         region 
##          FALSE          FALSE
\end{verbatim}

There do not appear to be any missing values in the data! (In general we might also check for old-fashioned missing value codes like ``9'', ``99'', ``98'', etc., but there's no evidence of that in Table 5 (look at the Min and Max values - no ``9's'', ``99's'', etc.))

\hypertarget{make-some-appropriate-descriptive-eda-plots-to-illustrate-any-important-features-of-the-variables-or-possible-important-relationships-among-them.}{%
\subsubsection{Make some appropriate descriptive EDA plots to illustrate any important features of the variables or possible important relationships among them.}\label{make-some-appropriate-descriptive-eda-plots-to-illustrate-any-important-features-of-the-variables-or-possible-important-relationships-among-them.}}

We want to capture (i) univariatoe distributions {[}primarily historgrams, boxplots, and the like{]} and (ii) two-variable relationships {[}primarily scatter plots{]}.

We could try \texttt{ggpairs} from \texttt{library(GGally)} to capture both (i) and (ii), but there are too many variables to make this a legible plot. Instead, we'll do (i) and (ii) separately.

For much of the rest of this problem, I want \texttt{region} to be in the data frame as well as the numerical variables, so I will create that now:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cdigood }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(cdinumeric,}\AttributeTok{region=}\NormalTok{cdidata}\SpecialCharTok{$}\NormalTok{region)}
\end{Highlighting}
\end{Shaded}

First, univariate distributions:

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# You can get almost the same thing as below with }
\DocumentationTok{\#\#}
\DocumentationTok{\#\# ggplot(gather(cdinumeric), aes(value)) + }
\DocumentationTok{\#\#    geom\_histogram(bins=30) + }
\DocumentationTok{\#\#    facet\_wrap(\textasciitilde{}key, scales = \textquotesingle{}free\_x\textquotesingle{})}

\NormalTok{hist.builder }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(df) \{  }\DocumentationTok{\#\# creates a list of graphs}
\NormalTok{    result }\OtherTok{\textless{}{-}} \ConstantTok{NULL} 
    \ControlFlowTok{for}\NormalTok{ (var }\ControlFlowTok{in} \FunctionTok{names}\NormalTok{(df)) \{}
\NormalTok{      d }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{dd=}\NormalTok{df[,var])}
      \ControlFlowTok{if}\NormalTok{(}\FunctionTok{mode}\NormalTok{(df[,var])}\SpecialCharTok{==}\StringTok{"numeric"}\NormalTok{) \{}
\NormalTok{        p }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(d,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{dd)) }\SpecialCharTok{+} \FunctionTok{geom\_histogram}\NormalTok{() }\SpecialCharTok{+}
               \FunctionTok{ggtitle}\NormalTok{(var) }\SpecialCharTok{+} \FunctionTok{xlab}\NormalTok{(}\StringTok{""}\NormalTok{) }
\NormalTok{      \} }\ControlFlowTok{else}\NormalTok{ \{}
\NormalTok{          p }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(d,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{dd)) }\SpecialCharTok{+} \FunctionTok{geom\_bar}\NormalTok{() }\SpecialCharTok{+}
                 \FunctionTok{ggtitle}\NormalTok{(var) }\SpecialCharTok{+} \FunctionTok{xlab}\NormalTok{(}\StringTok{""}\NormalTok{) }
\NormalTok{      \}}
\NormalTok{    result }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(result,}\FunctionTok{list}\NormalTok{(p))}
\NormalTok{    \}}
    \FunctionTok{return}\NormalTok{(result)}
\NormalTok{\}}

\FunctionTok{grid.arrange}\NormalTok{(}\AttributeTok{grobs=}\FunctionTok{hist.builder}\NormalTok{(cdigood))  }\DocumentationTok{\#\# from library(grid) \& library(gridExtra)}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-10-1.pdf}
\caption{\label{fig:unnamed-chunk-10}Distributions of Variables}
\end{figure}

It looks from the histograms like the variables that will really need attention (because they are severely right-skewed) are \texttt{land.area}, \texttt{pop}, \texttt{doctors}, \texttt{hosp.beds}, \texttt{crimes}, and \texttt{tot.income}, and maybe \texttt{per.cap.income}. These are the same variables that we identified from Table 5 above.

And now for bivariate relationships\ldots{} We really only care about two things:

\begin{itemize}
\item
  checking for high correlations among the predictors
\item
  checking to see if the predictors are linearly related to \texttt{per.cap.income}
\end{itemize}

When we have a large number of variables, it can be a more efficient use of the page to make a heatmap of the correlation matrix, rather than looking at the correlation matrix itself:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{corgraph }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(df) \{}
\NormalTok{  cormat }\OtherTok{\textless{}{-}} \FunctionTok{cor}\NormalTok{(df)}
\NormalTok{  melted\_cormat }\OtherTok{\textless{}{-}} \FunctionTok{melt}\NormalTok{(cormat)   }\DocumentationTok{\#\# need library(reshape2) for this...}
  \FunctionTok{ggplot}\NormalTok{(}\AttributeTok{data =}\NormalTok{ melted\_cormat, }\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{Var1, }\AttributeTok{y=}\NormalTok{Var2, }\AttributeTok{fill=}\NormalTok{value)) }\SpecialCharTok{+} 
    \FunctionTok{geom\_tile}\NormalTok{() }\SpecialCharTok{+} 
    \FunctionTok{theme}\NormalTok{(}\AttributeTok{axis.text.x =} \FunctionTok{element\_text}\NormalTok{(}\AttributeTok{angle =} \DecValTok{45}\NormalTok{,}\AttributeTok{vjust=}\FloatTok{0.9}\NormalTok{,}\AttributeTok{hjust=}\DecValTok{1}\NormalTok{)) }\SpecialCharTok{+}
    \FunctionTok{scale\_fill\_gradient2}\NormalTok{(}\AttributeTok{low=}\StringTok{"red"}\NormalTok{,}\AttributeTok{mid=}\StringTok{"white"}\NormalTok{,}\AttributeTok{high=}\StringTok{"blue"}\NormalTok{)}
\NormalTok{\}}

\FunctionTok{corgraph}\NormalTok{(cdinumeric)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-11-1.pdf}
We can make the following conclusions from the correlation matrix:

\begin{itemize}
\item
  \texttt{tot.income} and \texttt{pop} are highly correlated (no surprise there)
\item
  both are reasonably highly correlated with \texttt{crimes}, \texttt{hosp.beds} and \texttt{doctors}
\item
  the three variables \texttt{crimes}, \texttt{hosp.beds} and \texttt{doctors} seem strongly correlated with one another
\item
  \texttt{per.cap.income} isn't really highly correlated with anything, but the best possibilities seem to be \texttt{pct.hs.grad}, \texttt{pct.bach.deg} (postively correlated with \texttt{per.cap.income}) and \texttt{pct.below.pov}, \texttt{pct.unemp} (negatively correlated with \texttt{per.cap.income}); all four of these variables are moderately highly correlated with one another
\end{itemize}

These observations suggest that we may run into multi-collinearity problems when we start fitting models, but there is also some hope that we can make a good model for \texttt{per.cap.income}.

Now we turn to scatter plots, but we are just going to concentrate on relationships with \texttt{per.cap.income}:

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# I\textquotesingle{}m not sure how I\textquotesingle{}d do this directly with ggplot...}

\NormalTok{scatter.builder }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(df,}\AttributeTok{yvar=}\StringTok{"per.cap.income"}\NormalTok{) \{}
\NormalTok{    result }\OtherTok{\textless{}{-}} \ConstantTok{NULL} 
\NormalTok{    y.index }\OtherTok{\textless{}{-}} \FunctionTok{grep}\NormalTok{(yvar,}\FunctionTok{names}\NormalTok{(df))}
    \ControlFlowTok{for}\NormalTok{ (xvar }\ControlFlowTok{in} \FunctionTok{names}\NormalTok{(df)[}\SpecialCharTok{{-}}\NormalTok{y.index]) \{}
\NormalTok{      d }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{xx=}\NormalTok{df[,xvar],}\AttributeTok{yy=}\NormalTok{df[,y.index])}
      \ControlFlowTok{if}\NormalTok{(}\FunctionTok{mode}\NormalTok{(df[,xvar])}\SpecialCharTok{==}\StringTok{"numeric"}\NormalTok{) \{}
\NormalTok{        p }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(d,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{xx,}\AttributeTok{y=}\NormalTok{yy)) }\SpecialCharTok{+} \FunctionTok{geom\_point}\NormalTok{() }\SpecialCharTok{+}
               \FunctionTok{ggtitle}\NormalTok{(}\StringTok{""}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{xlab}\NormalTok{(xvar) }\SpecialCharTok{+} \FunctionTok{ylab}\NormalTok{(yvar) }
\NormalTok{      \} }\ControlFlowTok{else}\NormalTok{ \{}
\NormalTok{          p }\OtherTok{\textless{}{-}} \FunctionTok{ggplot}\NormalTok{(d,}\FunctionTok{aes}\NormalTok{(}\AttributeTok{x=}\NormalTok{xx,}\AttributeTok{y=}\NormalTok{yy)) }\SpecialCharTok{+} \FunctionTok{geom\_boxplot}\NormalTok{(}\AttributeTok{notch=}\NormalTok{F) }\SpecialCharTok{+}
                 \FunctionTok{ggtitle}\NormalTok{(}\StringTok{""}\NormalTok{) }\SpecialCharTok{+} \FunctionTok{xlab}\NormalTok{(xvar) }\SpecialCharTok{+} \FunctionTok{ylab}\NormalTok{(yvar) }
\NormalTok{      \}}
\NormalTok{    result }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(result,}\FunctionTok{list}\NormalTok{(p))}
\NormalTok{    \}}
    \FunctionTok{return}\NormalTok{(result)}
\NormalTok{\}}

\FunctionTok{grid.arrange}\NormalTok{(}\AttributeTok{grobs=}\FunctionTok{scatter.builder}\NormalTok{(cdigood))}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-12-1.pdf}
The best possibilities for predicting \texttt{per.cap.income} are the same variables we identified from the correlation matrix: \texttt{pct.hs.grad}, \texttt{pct.bach.deg},\texttt{pct.below.pov}, and \texttt{pct.unemp}. The last plot shows how \texttt{per.cap.income} varies across the four regions of the country. There is a lot of overlap in the boxplots, but the Northeast and the West seem to be doing a little better than the North Central and South regions.

\hypertarget{you-were-not-asked-to-do-any-transofrmations-at-this-point-but-it-is-a-good-idea-to-stop-and-think-about-transformations-for-a-minute-anyway}{%
\subsubsection{You were not asked to do any transofrmations at this point, but it is a good idea to stop and think about transformations for a minute anyway\ldots{}}\label{you-were-not-asked-to-do-any-transofrmations-at-this-point-but-it-is-a-good-idea-to-stop-and-think-about-transformations-for-a-minute-anyway}}

I am not going to work too hard on transformations for now. I am going to take the logarithms of the variables I identified earlier, to address heavy skewing and potential leverage and influence issues. But I will leave the other variables alone, even though there may be some skew in them. I can explain the logarithms to the social scientist in terms of percent-change concepts, as we have done in class. And the more untransformed variables there are, the easier it will be for the social scientist to think about the model(s) I present.

(Another plausible approach would be to observe that essentially all of the continuous variables are skewed to some extent, and so one could just log all of them for starters.)

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cdilogs }\OtherTok{\textless{}{-}}\NormalTok{ cdigood}

\NormalTok{skewed.vars }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"land.area"}\NormalTok{, }\StringTok{"pop"}\NormalTok{, }\StringTok{"doctors"}\NormalTok{, }\StringTok{"hosp.beds"}\NormalTok{, }\StringTok{"crimes"}\NormalTok{, }\StringTok{"tot.income"}\NormalTok{, }
                 \StringTok{"per.cap.income"}\NormalTok{)}

\ControlFlowTok{for}\NormalTok{ (tmp }\ControlFlowTok{in}\NormalTok{ skewed.vars) \{}
\NormalTok{  loc }\OtherTok{\textless{}{-}} \FunctionTok{grep}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"\^{}"}\NormalTok{,tmp,}\StringTok{"$"}\NormalTok{,}\AttributeTok{sep=}\StringTok{""}\NormalTok{),}\FunctionTok{names}\NormalTok{(cdilogs))}
    \DocumentationTok{\#\# note special characters "\^{}" and "$" to anchor}
    \DocumentationTok{\#\# the search pattern at the begining and end of the }
    \DocumentationTok{\#\# string}
\NormalTok{  cdilogs[,loc] }\OtherTok{\textless{}{-}} \FunctionTok{log}\NormalTok{(cdilogs[,loc])}
  \FunctionTok{names}\NormalTok{(cdilogs)[loc] }\OtherTok{\textless{}{-}} \FunctionTok{paste}\NormalTok{(}\StringTok{"log."}\NormalTok{,}\FunctionTok{names}\NormalTok{(cdilogs)[loc],}\AttributeTok{sep=}\StringTok{""}\NormalTok{)}
\NormalTok{\}}

\DocumentationTok{\#\# Now I get to re{-}use those crazy functions I wrote...}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{grid.arrange}\NormalTok{(}\AttributeTok{grobs=}\FunctionTok{hist.builder}\NormalTok{(cdilogs))}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-14-1.pdf}
\caption{\label{fig:unnamed-chunk-14}Histograms after log transformations}
\end{figure}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{corgraph}\NormalTok{(cdilogs[,}\SpecialCharTok{{-}}\FunctionTok{grep}\NormalTok{(}\StringTok{"region"}\NormalTok{,}\FunctionTok{names}\NormalTok{(cdilogs))])}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-15-1.pdf}
\caption{\label{fig:unnamed-chunk-15}Correlations after log transformations}
\end{figure}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{grid.arrange}\NormalTok{(}\AttributeTok{grobs=}\FunctionTok{scatter.builder}\NormalTok{(cdilogs,}\StringTok{"log.per.cap.income"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-16-1.pdf}
\caption{\label{fig:unnamed-chunk-16}Scatterplots after log transformations}
\end{figure}

Except for \texttt{log.pop} (which I probably won't use much anyway, since \texttt{per.cap.income} = \texttt{tot.income} / \texttt{pop} is a deterministic function of \texttt{pop}), the skewing seems to have largely been brought under control. The correlations are similar to, but a little stronger than, the correlations for the untransformed variables, and we can see from the scatter plots that the linear relationships we thought were there seem stronger, and there may be some other relationships to pay attention to.

\hypertarget{part-b}{%
\subsection{Part B}\label{part-b}}

\textbf{Build a regression model that predicts per-capita income from crime rate and region of the country. Should there be any interactions in the model? What does your model say about the relationship between per- capita income and crime rate? Do your answers change, depending on whether you use number of crimes, or ``per-capita crime'' = (number of crimes)/(population) as a crime rate measure? If so, which one best answers the question? Why? Show the fitted model results and explain your answer to these questions in terms of those results.}

For each version of the income and crime variables, there are essentially three models to think about. Since logarithms cleaned up a lot of the skewing in the data, I will only use log-transformed variables.

For log.per.cap.income and log.crimes, the three models to consider are:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{ancova}\FloatTok{.01} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ log.crimes,}\AttributeTok{data=}\NormalTok{cdilogs)}

\NormalTok{ancova}\FloatTok{.02} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ log.crimes }\SpecialCharTok{+}\NormalTok{ region,}\AttributeTok{data=}\NormalTok{cdilogs)}

\NormalTok{ancova}\FloatTok{.03} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ log.crimes }\SpecialCharTok{*}\NormalTok{ region,}\AttributeTok{data=}\NormalTok{cdilogs)}
\end{Highlighting}
\end{Shaded}

Residual plots appear in Figure 5 below. None of them are awful, so we could trust F-tests to compare these models:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{anova}\NormalTok{(ancova}\FloatTok{.01}\NormalTok{,ancova}\FloatTok{.02}\NormalTok{, ancova}\FloatTok{.03}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Variance Table
## 
## Model 1: log.per.cap.income ~ log.crimes
## Model 2: log.per.cap.income ~ log.crimes + region
## Model 3: log.per.cap.income ~ log.crimes * region
##   Res.Df    RSS Df Sum of Sq       F    Pr(>F)    
## 1    438 17.271                                   
## 2    435 14.949  3   2.32194 22.4823 1.523e-13 ***
## 3    432 14.872  3   0.07678  0.7434    0.5266    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

It looks like ancova.02 (the additive, or ``parallel lines'' ANCOVA model with no interactions) is doing the best.

In order to compare this with a model involving ``per-capita crime'' = (number of crimes)/(population), we have to construct the new variable:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{attach}\NormalTok{(cdigood)}
\NormalTok{per.cap.crime }\OtherTok{\textless{}{-}}\NormalTok{ crimes}\SpecialCharTok{/}\NormalTok{pop}
\NormalTok{log.per.cap.crimes }\OtherTok{\textless{}{-}}\FunctionTok{log}\NormalTok{(per.cap.crime)}
\FunctionTok{detach}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

Note that I could equally well have calculated

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{attach}\NormalTok{(cdigood)}
\NormalTok{log.per.cap.crimes }\OtherTok{\textless{}{-}}\FunctionTok{log}\NormalTok{(crimes) }\SpecialCharTok{{-}} \FunctionTok{log}\NormalTok{(pop)}
\FunctionTok{detach}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

Now let's look at the same three models with this new variable:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{ancova}\FloatTok{.04} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ log.per.cap.crimes,}\AttributeTok{data=}\NormalTok{cdilogs)}

\NormalTok{ancova}\FloatTok{.05} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ log.per.cap.crimes }\SpecialCharTok{+}\NormalTok{ region,}\AttributeTok{data=}\NormalTok{cdilogs)}

\NormalTok{ancova}\FloatTok{.06} \OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ log.per.cap.crimes }\SpecialCharTok{*}\NormalTok{ region,}\AttributeTok{data=}\NormalTok{cdilogs)}
\end{Highlighting}
\end{Shaded}

Again, the residual plots in Figure 5 look OK, so we try F-tests to compare these nested models

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{anova}\NormalTok{(ancova}\FloatTok{.04}\NormalTok{,ancova}\FloatTok{.05}\NormalTok{, ancova}\FloatTok{.06}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Variance Table
## 
## Model 1: log.per.cap.income ~ log.per.cap.crimes
## Model 2: log.per.cap.income ~ log.per.cap.crimes + region
## Model 3: log.per.cap.income ~ log.per.cap.crimes * region
##   Res.Df    RSS Df Sum of Sq       F    Pr(>F)    
## 1    438 18.697                                   
## 2    435 16.952  3   1.74465 14.8407 3.263e-09 ***
## 3    432 16.928  3   0.02408  0.2048     0.893    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

Now ancova.05 (once again, the additive or ``parallel lines'' model) is doing the best.

If we want to compare these two winners (ancova.02 vs.~ancova.05), we need to use AIC or BIC, because the two winners are not nested models (you can't get one from the other by imposing one or more linear constraints).

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{AIC}\NormalTok{(ancova}\FloatTok{.02}\NormalTok{,ancova}\FloatTok{.05}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##           df       AIC
## ancova.02  6 -227.4746
## ancova.05  6 -172.1347
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{BIC}\NormalTok{(ancova}\FloatTok{.02}\NormalTok{,ancova}\FloatTok{.05}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##           df       BIC
## ancova.02  6 -202.9539
## ancova.05  6 -147.6140
\end{verbatim}

And so it seems ancova.02 is the best model, by either AIC or BIC (remember smaller is better!).

Recalling the model formula for ancova.02 and a brief summary of the model

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{formula}\NormalTok{(ancova}\FloatTok{.02}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## log.per.cap.income ~ log.crimes + region
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{coef}\NormalTok{(}\FunctionTok{summary}\NormalTok{(ancova}\FloatTok{.02}\NormalTok{)),}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)     9.19       0.08  115.13     0.00
## log.crimes      0.07       0.01    7.92     0.00
## regionNE        0.10       0.03    4.09     0.00
## regionS        -0.09       0.02   -3.68     0.00
## regionW        -0.06       0.03   -1.96     0.05
\end{verbatim}

we can interpret the model as follows (using the handout ``log xform and percent interpretation.pdf'' from week03):

\begin{itemize}
\item
  All across the US, for every 1\% increase in crimes, we expect a 0.07\% increase in per-capita income, om average (this increase is statistically significant, but is it practically significant?).
\item
  Different regions of the country have different baseline per-capita incomes however: In the NC region, the baseline salary is \texttt{exp(9.19)\ =} \$9,798.65. In the NE it is \texttt{exp(9.19+0.010)\ =} \$10,829.18, and so forth, so in the S it is \$8,955.29, and in the W it is \$9,228.02. All of these region baselines are, according to the model, significantly different from the NC baseline.
\item
  So, according to the model, the \emph{level} of salary varies with region in the US, but the \emph{way it is related to crime} does not.
\end{itemize}

There is an argument that per-capita crime is more comparable to, or at least on the same scale as, per-capita income, so we will briefly look at the second-best model, anova.05, to see how it compares to ancova.02:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{formula}\NormalTok{(ancova}\FloatTok{.05}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## log.per.cap.income ~ log.per.cap.crimes + region
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{coef}\NormalTok{(}\FunctionTok{summary}\NormalTok{(ancova}\FloatTok{.05}\NormalTok{)),}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)            9.94       0.07  143.30     0.00
## log.per.cap.crimes     0.04       0.02    1.98     0.05
## regionNE               0.11       0.03    4.15     0.00
## regionS               -0.07       0.03   -2.84     0.00
## regionW               -0.02       0.03   -0.81     0.42
\end{verbatim}

The interpretation of this model is similar to that of ancova.02:

\begin{itemize}
\item
  All across the US, for every 1\% increase in per-capita crime, there is an associated 0.04\% increase in per-capita income (so, slightly smaller effect, but still statistically significant; should that matter?).
\item
  The regional baseline salaries are: NC: \$20,743.74, NE: \$23,155.79, S: is \$19,341.34, and W: \$20,332.99. All but the W region have baselines that are, according to the model, significantly different from the NC baseline.
\item
  Again, the level of salary varies with region, but not the way it varies with crime, according to the model.
\end{itemize}

One would have to choose which of these two models one would want to feature (and interpret!) in the main body of the IDMRAD paper\ldots{}

\clearpage

ASIDE:

We could also just have compared all 6 models with AIC or BIC. In this case, we get the same result: the additive model ancova.02 has both the lowest AIC and the lowest BIC values.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{data.frame}\NormalTok{(}\AttributeTok{AIC=}\FunctionTok{AIC}\NormalTok{(ancova}\FloatTok{.01}\NormalTok{,ancova}\FloatTok{.02}\NormalTok{,ancova}\FloatTok{.03}\NormalTok{,ancova}\FloatTok{.04}\NormalTok{,ancova}\FloatTok{.05}\NormalTok{,ancova}\FloatTok{.06}\NormalTok{),}
\AttributeTok{BIC=}\FunctionTok{BIC}\NormalTok{(ancova}\FloatTok{.01}\NormalTok{,ancova}\FloatTok{.02}\NormalTok{,ancova}\FloatTok{.03}\NormalTok{,ancova}\FloatTok{.04}\NormalTok{,ancova}\FloatTok{.05}\NormalTok{,ancova}\FloatTok{.06}\NormalTok{))[,}\SpecialCharTok{{-}}\DecValTok{3}\NormalTok{] }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T,}\AttributeTok{col.names=}\FunctionTok{c}\NormalTok{(}\StringTok{"df"}\NormalTok{,}\StringTok{"AIC"}\NormalTok{,}\StringTok{"BIC"}\NormalTok{)) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{(}\AttributeTok{full\_width=}\NormalTok{F)}
\end{Highlighting}
\end{Shaded}

\begin{table}
\centering
\begin{tabular}[t]{lrrr}
\toprule
  & df & AIC & BIC\\
\midrule
ancova.01 & 3 & -169.9466 & -157.6863\\
ancova.02 & 6 & -227.4746 & -202.9539\\
ancova.03 & 9 & -223.7402 & -186.9593\\
ancova.04 & 3 & -135.0340 & -122.7737\\
ancova.05 & 6 & -172.1347 & -147.6140\\
\addlinespace
ancova.06 & 9 & -166.7601 & -129.9792\\
\bottomrule
\end{tabular}
\end{table}

RESIDUAL PLOTS:

And finally, here are the residual disgnostic plots that justified our earlier use of Anova F-tests to compare nested models.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{oldmar }\OtherTok{\textless{}{-}} \FunctionTok{par}\NormalTok{()}\SpecialCharTok{$}\NormalTok{mar}

\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow=}\FunctionTok{c}\NormalTok{(}\DecValTok{6}\NormalTok{,}\DecValTok{4}\NormalTok{))}

\FunctionTok{par}\NormalTok{(}\AttributeTok{mar=}\FunctionTok{c}\NormalTok{(}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{))}

\FunctionTok{invisible}\NormalTok{(}\FunctionTok{lapply}\NormalTok{(}\FunctionTok{list}\NormalTok{(ancova}\FloatTok{.01}\NormalTok{,ancova}\FloatTok{.02}\NormalTok{,ancova}\FloatTok{.03}\NormalTok{,ancova}\FloatTok{.04}\NormalTok{,ancova}\FloatTok{.05}\NormalTok{,ancova}\FloatTok{.06}\NormalTok{),}
                \ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{plot}\NormalTok{(x,}\AttributeTok{cex.main=}\FloatTok{0.5}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

\begin{figure}
\centering
\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-29-1.pdf}
\caption{\label{fig:unnamed-chunk-29}Residual diagnostics for all 6 ANCOVA models, in order: ancova.01, ancova.02, ancova.03, ancova.04, ancova.05, ancova.06}
\end{figure}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mar=}\NormalTok{oldmar)}
\end{Highlighting}
\end{Shaded}

\hypertarget{part-c}{%
\subsection{Part C}\label{part-c}}

\textbf{Use methods we have discussed in class and/or methods from Sheather Chapters 5, 6 \& 7 (including, as needed: transformations, interactions, variable selection, residual analysis, fit indices, etc.) to find the multiple regression model predicting per-capita income from the other variables, that makes the ``best'' tradeoff between the following criteria:}

\begin{itemize}
\item
  \textbf{Reflects the social science and the meaning of the variables}
\item
  \textbf{Satisfies modeling assumptions}
\item
  \textbf{Clearly indicated by the data}
\item
  \textbf{Can be explained to someone who is more interested in social, economic and health factors than in mathematics and statistics.}
\end{itemize}

\textbf{No matter what you do, you are likely to be unhappy with some or all of these criteria; the better you make one criterion, the worse another is likely to get. So you will have to find a compromise or tradeoff between these criteria. Explain how you decided to make the tradeoff(s) you made.}

I am reasonably happy with the set of variables in the \texttt{cdilogs} data frame, so for now I will do model selection using just these variables. As shown in the first part of the answer to part (a), \texttt{id} and \texttt{county} are not useful predictors. \texttt{state} might be, but I am going to see how far I can get without \texttt{state}, and then perhaps incorprate it at the end.

Before beginning, I need to take log.pop and log.tot.income out of consideration, since per.cap.income is a deterministic function of them (so if they are included, no other predictors can possibly matter, and so I won't learn anything about what is associated with per.cap.income).

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{omit }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\FunctionTok{grep}\NormalTok{(}\StringTok{"log.pop"}\NormalTok{,}\FunctionTok{names}\NormalTok{(cdilogs)),}\FunctionTok{grep}\NormalTok{(}\StringTok{"log.tot.income"}\NormalTok{,}\FunctionTok{names}\NormalTok{(cdilogs)))}
\NormalTok{cdilogred }\OtherTok{\textless{}{-}}\NormalTok{ cdilogs[,}\SpecialCharTok{{-}}\NormalTok{omit]}
\end{Highlighting}
\end{Shaded}

Because \texttt{region} is categorical, and almost none of our variable selection functions are smart about dealing with the group of indicator variables associated with a categorical variable, I am going to first work \emph{without} the \texttt{region} variable, and then I will try to include the \texttt{region} variable in some intelligent way.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cdilogred.cont }\OtherTok{\textless{}{-}}\NormalTok{ cdilogred[,}\SpecialCharTok{{-}}\FunctionTok{grep}\NormalTok{(}\StringTok{"region"}\NormalTok{,}\FunctionTok{names}\NormalTok{(cdilogred))]}
\DocumentationTok{\#\# a data set to work with that doesn\textquotesingle{}t have \textquotesingle{}region\textquotesingle{} in it {-}{-} only the continuous variables...}
\FunctionTok{names}\NormalTok{(cdilogred.cont)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] "log.land.area"      "pop.18_34"          "pop.65_plus"       
##  [4] "log.doctors"        "log.hosp.beds"      "log.crimes"        
##  [7] "pct.hs.grad"        "pct.bach.deg"       "pct.below.pov"     
## [10] "pct.unemp"          "log.per.cap.income"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# there are 10 predictor variables and the response variable log.per.cap.income...}
\end{Highlighting}
\end{Shaded}

ALL-SUBSETS:

I'll start with all-subsets regression.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(leaps)  }\DocumentationTok{\#\# for the regsubsets() function}
\FunctionTok{library}\NormalTok{(car)    }\DocumentationTok{\#\# for the subsets() plotting function}
\end{Highlighting}
\end{Shaded}

A little preliminary exploration convinces me of two things:

\begin{itemize}
\item
  There are too many variables for the subsets() function to be an effective way to summarize the results of regsubsets()
\item
  The documentation and a bit of googling shows that I have some other tools to work with:

  \begin{itemize}
  \tightlist
  \item
    plot() makes a \emph{different} graphical summary of regsubsets results that \emph{is useful}
  \item
    summary() of regsubsets results can be stored in a variable, and there are several named components which can be used to locate the ``best'' (i.e.~lowest-BIC) model
  \item
    coef() works on regsubsets results, if you can identify which model in the sequence of models from regsubsets you want
  \end{itemize}
\end{itemize}

So\ldots{} onward!

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{all.subsets}\FloatTok{.01} \OtherTok{\textless{}{-}} \FunctionTok{regsubsets}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ ., }\AttributeTok{data=}\NormalTok{cdilogred.cont,}\AttributeTok{nvmax=}\DecValTok{10}\NormalTok{)}
\FunctionTok{plot}\NormalTok{(all.subsets}\FloatTok{.01}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-33-1.pdf}
In the plot, the dark squares indicate which variables are in the model that has the BIC values on the left. The darker the squares, the better the model (remember, small is good for BIC). Rather than read them off the plot, though, we can use the components of summary() to get them directly.

We begin by locating the model with the lowest BIC\ldots{}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{all.subsets.}\FloatTok{01.}\NormalTok{summary }\OtherTok{\textless{}{-}} \FunctionTok{summary}\NormalTok{(all.subsets}\FloatTok{.01}\NormalTok{)}
\FunctionTok{names}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{summary)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] "which"  "rsq"    "rss"    "adjr2"  "cp"     "bic"    "outmat" "obj"
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{all.subsets.}\FloatTok{01.}\NormalTok{summary}\SpecialCharTok{$}\NormalTok{bic}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##  [1] -257.5260 -502.4302 -572.5538 -682.8532 -732.1894 -761.5908 -772.0715
##  [8] -770.5990 -766.2235 -760.4131
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{min}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{summary}\SpecialCharTok{$}\NormalTok{bic)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] -772.0715
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{print}\NormalTok{(best.model }\OtherTok{\textless{}{-}} \FunctionTok{which.min}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{summary}\SpecialCharTok{$}\NormalTok{bic))}
\end{Highlighting}
\end{Shaded}

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

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{coef}\NormalTok{(all.subsets}\FloatTok{.01}\NormalTok{,best.model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   (Intercept) log.land.area     pop.18_34   log.doctors   pct.hs.grad 
##  10.222495041  -0.035674062  -0.013900201   0.060676872  -0.004406396 
##  pct.bach.deg pct.below.pov     pct.unemp 
##   0.015385301  -0.024278371   0.010603691
\end{verbatim}

Note that the variables listed here are the same as the ones at the top of the graph we made.

Before we interpret the coefficients, let's refit the model so we can get standard errors, etc. To make life easy, we'll set up a temporary data frame with just the variables from the lowest-BIC model.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{all.subsets.}\FloatTok{01.}\NormalTok{summary}\SpecialCharTok{$}\NormalTok{which[best.model,]}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   (Intercept) log.land.area     pop.18_34   pop.65_plus   log.doctors 
##          TRUE          TRUE          TRUE         FALSE          TRUE 
## log.hosp.beds    log.crimes   pct.hs.grad  pct.bach.deg pct.below.pov 
##         FALSE         FALSE          TRUE          TRUE          TRUE 
##     pct.unemp 
##          TRUE
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tmp }\OtherTok{\textless{}{-}}\NormalTok{ cdilogred.cont[,all.subsets.}\FloatTok{01.}\NormalTok{summary}\SpecialCharTok{$}\NormalTok{which[best.model,][}\SpecialCharTok{{-}}\DecValTok{1}\NormalTok{]]}
\NormalTok{all.subsets.}\FloatTok{01.}\NormalTok{final.model }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ .,}\AttributeTok{data=}\NormalTok{tmp)}
\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model)}\SpecialCharTok{$}\NormalTok{coef}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                   Estimate   Std. Error    t value      Pr(>|t|)
## (Intercept)   10.222495041 0.0931210074 109.776465 1.127483e-317
## log.land.area -0.035674062 0.0047767371  -7.468291  4.533156e-13
## pop.18_34     -0.013900201 0.0011113007 -12.508046  7.514862e-31
## log.doctors    0.060676872 0.0040183327  15.100012  1.133432e-41
## pct.hs.grad   -0.004406396 0.0010822796  -4.071403  5.558448e-05
## pct.bach.deg   0.015385301 0.0009245509  16.640838  2.100590e-48
## pct.below.pov -0.024278371 0.0012583372 -19.294011  2.812246e-60
## pct.unemp      0.010603691 0.0021771148   4.870525  1.564524e-06
\end{verbatim}

All the predictors have coefficients significantly different from zero. However, most of the coefficients are small, and some seem to have the wrong sign (e.g.~\texttt{pct.hs.grad} and \texttt{pct.unemp}). Before diving into what might be going on, let's at least check VIFs and residual diagnostics\ldots{}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{vif}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## log.land.area     pop.18_34   log.doctors   pct.hs.grad  pct.bach.deg 
##      1.131867      1.416145      1.379671      3.763103      3.269565 
## pct.below.pov     pct.unemp 
##      2.241555      1.691280
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow=}\FunctionTok{c}\NormalTok{(}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{))}
\FunctionTok{plot}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-36-1.pdf}
None of the VIFs seem excessively large, and the diagnostic plots don't show much except that the QQ plot suggests both the left and the right tails are a bit longer than expected for the normal distribution.

We can also look at marginal model plots, to see if we are missing a transformation, interaction, etc.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{mmps}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-37-1.pdf}
The marginal model plots look very good -- the blue data-based curves line up well with the red model-based curves. We don't seem to be missing any important transformations, interactions, etc.

One last thing to try is to see if interaction with \texttt{region} helps in any way.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tmp }\OtherTok{\textless{}{-}} \FunctionTok{cbind}\NormalTok{(tmp,}\AttributeTok{region=}\NormalTok{cdilogred}\SpecialCharTok{$}\NormalTok{region)}
\NormalTok{all.subsets.}\FloatTok{01.}\NormalTok{final.with.region }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ .}\SpecialCharTok{*}\NormalTok{region,}\AttributeTok{data=}\NormalTok{tmp)}
\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = log.per.cap.income ~ . * region, data = tmp)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.250782 -0.042332 -0.002298  0.040559  0.313570 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            10.1244260  0.2826240  35.823  < 2e-16 ***
## log.land.area          -0.0364187  0.0151355  -2.406 0.016564 *  
## pop.18_34              -0.0147940  0.0026043  -5.681 2.55e-08 ***
## log.doctors             0.0544169  0.0093221   5.837 1.08e-08 ***
## pct.hs.grad            -0.0024773  0.0034110  -0.726 0.468088    
## pct.bach.deg            0.0140833  0.0029254   4.814 2.09e-06 ***
## pct.below.pov          -0.0237085  0.0036234  -6.543 1.81e-10 ***
## pct.unemp               0.0180393  0.0048923   3.687 0.000257 ***
## regionNE                0.3243992  0.3577081   0.907 0.365004    
## regionS                -0.0345856  0.3131668  -0.110 0.912116    
## regionW                 1.5043946  0.4226868   3.559 0.000416 ***
## log.land.area:regionNE -0.0037179  0.0201435  -0.185 0.853656    
## log.land.area:regionS  -0.0047582  0.0174155  -0.273 0.784825    
## log.land.area:regionW   0.0151234  0.0181871   0.832 0.406154    
## pop.18_34:regionNE     -0.0024780  0.0036873  -0.672 0.501939    
## pop.18_34:regionS      -0.0008777  0.0030680  -0.286 0.774970    
## pop.18_34:regionW       0.0014122  0.0040925   0.345 0.730220    
## log.doctors:regionNE   -0.0046251  0.0132571  -0.349 0.727359    
## log.doctors:regionS     0.0043337  0.0114401   0.379 0.705019    
## log.doctors:regionW    -0.0034863  0.0131576  -0.265 0.791173    
## pct.hs.grad:regionNE   -0.0037529  0.0044150  -0.850 0.395813    
## pct.hs.grad:regionS     0.0021198  0.0037853   0.560 0.575790    
## pct.hs.grad:regionW    -0.0190188  0.0045881  -4.145 4.13e-05 ***
## pct.bach.deg:regionNE   0.0069429  0.0040312   1.722 0.085776 .  
## pct.bach.deg:regionS   -0.0015774  0.0032000  -0.493 0.622328    
## pct.bach.deg:regionW    0.0071026  0.0036374   1.953 0.051541 .  
## pct.below.pov:regionNE -0.0014134  0.0050896  -0.278 0.781381    
## pct.below.pov:regionS   0.0072764  0.0040739   1.786 0.074827 .  
## pct.below.pov:regionW  -0.0161639  0.0054271  -2.978 0.003071 ** 
## pct.unemp:regionNE     -0.0083596  0.0073758  -1.133 0.257720    
## pct.unemp:regionS      -0.0249396  0.0065867  -3.786 0.000176 ***
## pct.unemp:regionW      -0.0201466  0.0067713  -2.975 0.003101 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0759 on 408 degrees of freedom
## Multiple R-squared:  0.8747, Adjusted R-squared:  0.8652 
## F-statistic: 91.91 on 31 and 408 DF,  p-value: < 2.2e-16
\end{verbatim}

A reasonable rule of thumb is: if any indicator for a categorical variable seems important (e.g.~a statistically significant coefficient), then keep the whole categorical variable. If none of them seem important, then drop the variable. The same thing works for interactions with categorical variables.

\emph{(A more formally correct thing to do would be to check to see whether to include each group of variables involving \texttt{region} with partial F tests (the anova() function), AIC, and/or BIC, in a kind of stepwise regression procedure.)}

Using the rule of thumb, we would

\begin{itemize}
\item
  Keep: \texttt{region}, \texttt{region:pct.hs.grad}, \texttt{region:pct.below.pov}, \texttt{region:pct.unemp}, and maybe \texttt{region:pct.bach.deg}
\item
  Drop: \texttt{region:log.land.area}, \texttt{region:pop.18\_34}, \texttt{region:log.doctors}
\end{itemize}

Thus we arrive at

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region }\OtherTok{\textless{}{-}} \FunctionTok{update}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.region,}
\NormalTok{                                                . }\SpecialCharTok{\textasciitilde{}}\NormalTok{ . }\SpecialCharTok{{-}}\NormalTok{ region}\SpecialCharTok{:}\NormalTok{log.land.area }\SpecialCharTok{{-}} 
\NormalTok{                                                  region}\SpecialCharTok{:}\NormalTok{pop}\FloatTok{.18}\NormalTok{\_34 }\SpecialCharTok{{-}}\NormalTok{ region}\SpecialCharTok{:}\NormalTok{log.doctors)}
\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = log.per.cap.income ~ log.land.area + pop.18_34 + 
##     log.doctors + pct.hs.grad + pct.bach.deg + pct.below.pov + 
##     pct.unemp + region + pct.hs.grad:region + pct.bach.deg:region + 
##     pct.below.pov:region + pct.unemp:region, data = tmp)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.268015 -0.043459 -0.002511  0.039967  0.313939 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            10.125260   0.251582  40.246  < 2e-16 ***
## log.land.area          -0.034569   0.005376  -6.430 3.50e-10 ***
## pop.18_34              -0.015404   0.001087 -14.170  < 2e-16 ***
## log.doctors             0.055342   0.004034  13.720  < 2e-16 ***
## pct.hs.grad            -0.002503   0.003151  -0.794 0.427456    
## pct.bach.deg            0.014208   0.002108   6.741 5.24e-11 ***
## pct.below.pov          -0.023634   0.003351  -7.054 7.30e-12 ***
## pct.unemp               0.017787   0.004783   3.719 0.000228 ***
## regionNE                0.219429   0.302526   0.725 0.468661    
## regionS                -0.062648   0.276125  -0.227 0.820627    
## regionW                 1.629351   0.357633   4.556 6.86e-06 ***
## pct.hs.grad:regionNE   -0.003640   0.003876  -0.939 0.348271    
## pct.hs.grad:regionS     0.002014   0.003539   0.569 0.569690    
## pct.hs.grad:regionW    -0.018916   0.004204  -4.499 8.85e-06 ***
## pct.bach.deg:regionNE   0.005905   0.002618   2.256 0.024611 *  
## pct.bach.deg:regionS   -0.001298   0.002321  -0.559 0.576352    
## pct.bach.deg:regionW    0.006326   0.002620   2.415 0.016183 *  
## pct.below.pov:regionNE -0.002435   0.004647  -0.524 0.600488    
## pct.below.pov:regionS   0.007137   0.003686   1.937 0.053482 .  
## pct.below.pov:regionW  -0.015224   0.005169  -2.945 0.003407 ** 
## pct.unemp:regionNE     -0.007967   0.007255  -1.098 0.272761    
## pct.unemp:regionS      -0.024668   0.006377  -3.868 0.000127 ***
## pct.unemp:regionW      -0.019757   0.006603  -2.992 0.002935 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07545 on 417 degrees of freedom
## Multiple R-squared:  0.8735, Adjusted R-squared:  0.8668 
## F-statistic: 130.9 on 22 and 417 DF,  p-value: < 2.2e-16
\end{verbatim}

This seems to be working\ldots{} All of the main effects and interactions that involve \texttt{region} have at least one significant term, and the \(R^2\) and residual SE hardly moved at all.

Let's quickly check VIFs and diagnostics, and then compare to the best BIC model that we obtained from stepwise without the \texttt{region} term.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{vif}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                              GVIF Df GVIF^(1/(2*Df))
## log.land.area        1.693937e+00  1        1.301513
## pop.18_34            1.600959e+00  1        1.265290
## log.doctors          1.642537e+00  1        1.281615
## pct.hs.grad          3.768383e+01  1        6.138716
## pct.bach.deg         2.007247e+01  1        4.480231
## pct.below.pov        1.877727e+01  1        4.333275
## pct.unemp            9.644220e+00  1        3.105515
## region               4.181419e+08  3       27.345589
## pct.hs.grad:region   3.373864e+08  3       26.384848
## pct.bach.deg:region  2.600598e+04  3        5.443090
## pct.below.pov:region 6.114866e+03  3        4.276264
## pct.unemp:region     1.364209e+04  3        4.888176
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow=}\FunctionTok{c}\NormalTok{(}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{))}
\FunctionTok{plot}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-40-1.pdf}

As usual, adding interactions usually creates collinearities (the interactions are necessarily somewhat collinear with their main effects). We worry less about this, because the collinearities do not appear to be causing the t-statistics and p-values in the model summary to behave strangely.

The diagnostic plots look about as good or as bad as the last set we looked at. It seems to me that in this analysis, the residual diagnostic plots are always going the toughest to feel totally satisfied about.

Now let's compare our model with some region terms to the model we got from all-subsets, but with region not included:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{anova}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Variance Table
## 
## Model 1: log.per.cap.income ~ log.land.area + pop.18_34 + log.doctors + 
##     pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp
## Model 2: log.per.cap.income ~ log.land.area + pop.18_34 + log.doctors + 
##     pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp + 
##     region + pct.hs.grad:region + pct.bach.deg:region + pct.below.pov:region + 
##     pct.unemp:region
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    432 2.9051                                  
## 2    417 2.3736 15   0.53148 6.2249 6.673e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{AIC}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                       df       AIC
## all.subsets.01.final.model             9  -942.274
## all.subsets.01.final.with.some.region 24 -1001.178
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{BIC}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                                       df       BIC
## all.subsets.01.final.model             9 -905.4931
## all.subsets.01.final.with.some.region 24 -903.0957
\end{verbatim}

The anova (F test) and AIC really like the model with the region terms in it! On the other hand, BIC prefers the simpler model. This often happens: BIC chooses a simpler model than other methods, because it is going for a parsimonius explanatory model, rather than a highly predictive model.

Finally, let's examine the table of estimated coefficients again for the model with some region terms in it, and offer some interpretations\ldots{}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{formula}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## log.per.cap.income ~ log.land.area + pop.18_34 + log.doctors + 
##     pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp + 
##     region + pct.hs.grad:region + pct.bach.deg:region + pct.below.pov:region + 
##     pct.unemp:region
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}\SpecialCharTok{$}\NormalTok{coef,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)               10.13       0.25   40.25     0.00
## log.land.area             -0.03       0.01   -6.43     0.00
## pop.18_34                 -0.02       0.00  -14.17     0.00
## log.doctors                0.06       0.00   13.72     0.00
## pct.hs.grad                0.00       0.00   -0.79     0.43
## pct.bach.deg               0.01       0.00    6.74     0.00
## pct.below.pov             -0.02       0.00   -7.05     0.00
## pct.unemp                  0.02       0.00    3.72     0.00
## regionNE                   0.22       0.30    0.73     0.47
## regionS                   -0.06       0.28   -0.23     0.82
## regionW                    1.63       0.36    4.56     0.00
## pct.hs.grad:regionNE       0.00       0.00   -0.94     0.35
## pct.hs.grad:regionS        0.00       0.00    0.57     0.57
## pct.hs.grad:regionW       -0.02       0.00   -4.50     0.00
## pct.bach.deg:regionNE      0.01       0.00    2.26     0.02
## pct.bach.deg:regionS       0.00       0.00   -0.56     0.58
## pct.bach.deg:regionW       0.01       0.00    2.41     0.02
## pct.below.pov:regionNE     0.00       0.00   -0.52     0.60
## pct.below.pov:regionS      0.01       0.00    1.94     0.05
## pct.below.pov:regionW     -0.02       0.01   -2.95     0.00
## pct.unemp:regionNE        -0.01       0.01   -1.10     0.27
## pct.unemp:regionS         -0.02       0.01   -3.87     0.00
## pct.unemp:regionW         -0.02       0.01   -2.99     0.00
\end{verbatim}

Here are some things we can say (this is only a few examples, there are more things to say!):

\begin{itemize}
\item
  For every 1\% increase in a county's land area, there is a 0.03\% decrease in expected per-capita income. (We might conjecture that this could be due to an urban-rural contrast: rural counties tend to be bigger than urban ones).
\item
  For every 1\% increase in the number of doctors in a county, the expected per-capita income increases by about 0.06\%. That makes sense; doctors are well-paid and could be big contributors to the per-capita average income.
\item
  For every 1 percentage point increase in the percent of the population aged 18--34, there iz an expected 2\% drop in per-capita income. (We might conjecture that this is because 18--34 year olds are not at peak earning capacity yet and so perhaps their lower incomes drags down the per-capita average).
\item
  percent of the population that are high school graduates doesn't have much effect, except in the South, where a one percentage point increase in hs graduates induces an expected 2\% decrease in per-capita income. I don't hae a good explanation for this: It might depend on whether college graduates are counted as a subset of hs graduates rather than counting them separately, or it might have something to do with some unique feature of economics in the southern region of the US.
\item
  In the main effect for region, and in several of the interactions for region, the West shows up as deviating significantly from the North Central part of the US
\item
  etc.
\end{itemize}

STEPWISE REGRESSION:

Now let's try to see how far we get with stepwise regression. I am going to take the same strategy: first try stepwise without \texttt{region}, and then see if adding \texttt{region} helps.

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

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{stepwise.base }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(log.per.cap.income }\SpecialCharTok{\textasciitilde{}}\NormalTok{ .,}\AttributeTok{data=}\NormalTok{cdilogred.cont)}
\DocumentationTok{\#\# try to duplicate all{-}subsets with BIC}
\NormalTok{step.result.}\FloatTok{01.}\NormalTok{bic }\OtherTok{\textless{}{-}} \FunctionTok{stepAIC}\NormalTok{(stepwise.base,}
                       \AttributeTok{scope=}\FunctionTok{list}\NormalTok{(}\AttributeTok{lower =} \SpecialCharTok{\textasciitilde{}} \DecValTok{1}\NormalTok{, }\AttributeTok{upper =} \SpecialCharTok{\textasciitilde{}}\NormalTok{ .), }
                       \AttributeTok{k=}\FunctionTok{log}\NormalTok{(}\FunctionTok{dim}\NormalTok{(cdilogred.cont)[}\DecValTok{1}\NormalTok{]),}
                       \AttributeTok{trace=}\NormalTok{F)}

\FunctionTok{anova}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,step.result.}\FloatTok{01.}\NormalTok{bic)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Variance Table
## 
## Model 1: log.per.cap.income ~ log.land.area + pop.18_34 + log.doctors + 
##     pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp
## Model 2: log.per.cap.income ~ log.land.area + pop.18_34 + log.doctors + 
##     pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp
##   Res.Df    RSS Df Sum of Sq F Pr(>F)
## 1    432 2.9051                      
## 2    432 2.9051  0         0
\end{verbatim}

We can see that stepwise regression using the BIC criterion actually found exactly the all-subsets model, which is great.

\begin{Shaded}
\begin{Highlighting}[]
\DocumentationTok{\#\# now try with AIC}
\NormalTok{step.result.}\FloatTok{01.}\NormalTok{aic }\OtherTok{\textless{}{-}} \FunctionTok{stepAIC}\NormalTok{(stepwise.base,}
                       \AttributeTok{scope=}\FunctionTok{list}\NormalTok{(}\AttributeTok{lower =} \SpecialCharTok{\textasciitilde{}} \DecValTok{1}\NormalTok{, }\AttributeTok{upper =} \SpecialCharTok{\textasciitilde{}}\NormalTok{ .), }
                       \AttributeTok{k=}\DecValTok{2}\NormalTok{,}
                       \AttributeTok{trace=}\NormalTok{F)}

\FunctionTok{anova}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,step.result.}\FloatTok{01.}\NormalTok{aic)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## Analysis of Variance Table
## 
## Model 1: log.per.cap.income ~ log.land.area + pop.18_34 + log.doctors + 
##     pct.hs.grad + pct.bach.deg + pct.below.pov + pct.unemp
## Model 2: log.per.cap.income ~ log.land.area + pop.18_34 + pop.65_plus + 
##     log.doctors + pct.hs.grad + pct.bach.deg + pct.below.pov + 
##     pct.unemp
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1    432 2.9051                              
## 2    431 2.8748  1  0.030306 4.5437 0.03361 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}

Stepwise regression using the AIC criterion added one more predictor, namely \texttt{pop.65\_plus}. Looking at the coefficient estimates below, we can say that, although the coefficient on \texttt{pop.65\_plus} is significantly different from zero, its effect on expected per-capita income appears to be quite small. Without knowing more about the economics, it's hard to make a decision about whether to keep \texttt{pop.65\_plus} or not.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{01.}\NormalTok{aic)}\SpecialCharTok{$}\NormalTok{coef,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)      10.32        0.1  100.56     0.00
## log.land.area    -0.04        0.0   -7.65     0.00
## pop.18_34        -0.02        0.0  -11.82     0.00
## pop.65_plus       0.00        0.0   -2.13     0.03
## log.doctors       0.06        0.0   15.26     0.00
## pct.hs.grad       0.00        0.0   -4.30     0.00
## pct.bach.deg      0.02        0.0   16.46     0.00
## pct.below.pov    -0.02        0.0  -19.49     0.00
## pct.unemp         0.01        0.0    4.96     0.00
\end{verbatim}

Well, now let's go for broke and see what happens if we just explore 2-way interactions briefly.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{step.result.}\FloatTok{02.}\NormalTok{bic }\OtherTok{\textless{}{-}} \FunctionTok{stepAIC}\NormalTok{(stepwise.base,}\AttributeTok{scope=}\FunctionTok{list}\NormalTok{(}\AttributeTok{lower =} \SpecialCharTok{\textasciitilde{}} \DecValTok{1}\NormalTok{, }\AttributeTok{upper =} \SpecialCharTok{\textasciitilde{}}\NormalTok{ .}\SpecialCharTok{\^{}}\DecValTok{2}\NormalTok{), }
                      \AttributeTok{k=}\FunctionTok{log}\NormalTok{(}\FunctionTok{dim}\NormalTok{(cdilogred.cont)[}\DecValTok{1}\NormalTok{]),                         }\DocumentationTok{\#\# BIC penalty.}
                      \AttributeTok{trace=}\NormalTok{F)}

\NormalTok{step.result.}\FloatTok{02.}\NormalTok{aic }\OtherTok{\textless{}{-}} \FunctionTok{stepAIC}\NormalTok{(stepwise.base,}\AttributeTok{scope=}\FunctionTok{list}\NormalTok{(}\AttributeTok{lower =} \SpecialCharTok{\textasciitilde{}} \DecValTok{1}\NormalTok{, }\AttributeTok{upper =} \SpecialCharTok{\textasciitilde{}}\NormalTok{ .}\SpecialCharTok{\^{}}\DecValTok{2}\NormalTok{), }
                      \AttributeTok{k=}\DecValTok{2}\NormalTok{,                                                   }\DocumentationTok{\#\# AIC penalty.}
                      \AttributeTok{trace=}\NormalTok{F)}

\NormalTok{comparison }\OtherTok{\textless{}{-}} \FunctionTok{cbind}\NormalTok{(}
\FunctionTok{AIC}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,step.result.}\FloatTok{01.}\NormalTok{aic,step.result.}\FloatTok{01.}\NormalTok{bic,}
\NormalTok{    step.result.}\FloatTok{02.}\NormalTok{aic,step.result.}\FloatTok{02.}\NormalTok{bic),}
\FunctionTok{BIC}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.model,step.result.}\FloatTok{01.}\NormalTok{aic,step.result.}\FloatTok{01.}\NormalTok{bic,}
\NormalTok{    step.result.}\FloatTok{02.}\NormalTok{aic,step.result.}\FloatTok{02.}\NormalTok{bic))}
\NormalTok{comparison }\OtherTok{\textless{}{-}}\NormalTok{ comparison[,}\SpecialCharTok{{-}}\DecValTok{3}\NormalTok{]}
\FunctionTok{names}\NormalTok{(comparison) }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"df"}\NormalTok{,}\StringTok{"AIC"}\NormalTok{,}\StringTok{"BIC"}\NormalTok{)}
\NormalTok{comparison }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kbl}\NormalTok{(}\AttributeTok{booktabs=}\NormalTok{T) }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{kable\_classic}\NormalTok{()}
\end{Highlighting}
\end{Shaded}

\begin{table}
\centering
\begin{tabular}[t]{lrrr}
\toprule
  & df & AIC & BIC\\
\midrule
all.subsets.01.final.model & 9 & -942.2740 & -905.4931\\
step.result.01.aic & 10 & -944.8883 & -904.0206\\
step.result.01.bic & 9 & -942.2740 & -905.4931\\
step.result.02.aic & 27 & -1064.7253 & -954.3824\\
step.result.02.bic & 12 & -1020.6026 & -971.5613\\
\bottomrule
\end{tabular}
\end{table}

As you can see in the table, both AIC and BIC like models with some interactions. This is a quandary: now the model is getting complicated to explain to someone.

\newpage

Let's look at the two models with interactons and see what we can say:

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{02.}\NormalTok{bic)}\SpecialCharTok{$}\NormalTok{coef,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                          Estimate Std. Error t value Pr(>|t|)
## (Intercept)                 10.83       0.13   83.19     0.00
## log.land.area               -0.03       0.00   -7.53     0.00
## pop.18_34                   -0.02       0.00  -13.09     0.00
## pop.65_plus                 -0.02       0.00   -8.81     0.00
## log.doctors                  0.05       0.00   13.29     0.00
## pct.hs.grad                 -0.01       0.00   -5.32     0.00
## pct.bach.deg                 0.01       0.00    3.01     0.00
## pct.below.pov               -0.02       0.00  -20.42     0.00
## pct.unemp                   -0.02       0.01   -1.83     0.07
## pop.65_plus:pct.bach.deg     0.00       0.00    8.33     0.00
## pct.hs.grad:pct.unemp        0.00       0.00    2.73     0.01
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{R2 = "}\NormalTok{,}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{02.}\NormalTok{bic)}\SpecialCharTok{$}\NormalTok{r.squared)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## R2 =  0.8721672
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{R2adj = "}\NormalTok{,}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{02.}\NormalTok{bic)}\SpecialCharTok{$}\NormalTok{adj.r.squared)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## R2adj =  0.8691874
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{02.}\NormalTok{aic)}\SpecialCharTok{$}\NormalTok{coef,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                             Estimate Std. Error t value Pr(>|t|)
## (Intercept)                    10.24       0.79   12.98     0.00
## log.land.area                   0.21       0.08    2.50     0.01
## pop.18_34                      -0.03       0.01   -4.73     0.00
## pop.65_plus                    -0.02       0.01   -1.92     0.06
## log.doctors                    -0.17       0.07   -2.53     0.01
## log.hosp.beds                  -0.11       0.10   -1.05     0.29
## log.crimes                      0.03       0.02    1.36     0.17
## pct.hs.grad                     0.00       0.01   -0.06     0.95
## pct.bach.deg                    0.06       0.01    4.60     0.00
## pct.below.pov                  -0.03       0.02   -1.52     0.13
## pct.unemp                      -0.01       0.01   -1.03     0.30
## pop.65_plus:pct.bach.deg        0.00       0.00    6.89     0.00
## log.land.area:pct.below.pov     0.00       0.00   -2.59     0.01
## log.land.area:pct.hs.grad       0.00       0.00   -2.39     0.02
## pct.hs.grad:pct.unemp           0.00       0.00    2.17     0.03
## pct.bach.deg:pct.below.pov      0.00       0.00   -5.81     0.00
## pop.18_34:log.doctors           0.01       0.00    2.74     0.01
## pct.hs.grad:pct.bach.deg        0.00       0.00   -4.10     0.00
## pct.bach.deg:pct.unemp          0.00       0.00   -1.80     0.07
## log.hosp.beds:pct.below.pov     0.01       0.00    4.52     0.00
## pop.18_34:pct.below.pov         0.00       0.00    3.29     0.00
## log.crimes:pct.below.pov        0.00       0.00   -2.82     0.01
## log.hosp.beds:pct.hs.grad       0.00       0.00    2.69     0.01
## pop.18_34:log.hosp.beds         0.00       0.00   -1.68     0.09
## log.land.area:pop.65_plus       0.00       0.00   -1.89     0.06
## pop.65_plus:log.crimes          0.00       0.00    1.60     0.11
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{R2 = "}\NormalTok{,}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{02.}\NormalTok{aic)}\SpecialCharTok{$}\NormalTok{r.squared)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## R2 =  0.8919859
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{R2adj = "}\NormalTok{,}\FunctionTok{summary}\NormalTok{(step.result.}\FloatTok{02.}\NormalTok{aic)}\SpecialCharTok{$}\NormalTok{adj.r.squared)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## R2adj =  0.8854633
\end{verbatim}

In the BIC-based model, there are two interactions that appear to be statistically significant, but their practical effect is almost zero on per-capita income.

The story is similar in the AIC model, which has many more interaction terms, but only two with an effect as large as 0.01.

Although both interactions models produced big jumps in AIC and BIC (much bigger than 10!), the improvement in \(R^2\) and \(R^2_{adj}\) is pretty small, for all the terms that have been added to the models.

For these reasons I might be willing to discuss these interactions with the social scientist, but I am disinclined to include them in a final model.

If I stick with the model found by all-subsets and stepAIC with a BIC penalty, then my conclusions about adding interactions with \texttt{region} will also be the same, and I will once again be led to \texttt{all.subsets.01.final.with.some.region}, which has some interesting and mostly-interpretable structure.

LASSO:

The last variable selection method I'll consider is the lasso.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(glmnet)}
\FunctionTok{library}\NormalTok{(arm)}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{loc }\OtherTok{\textless{}{-}} \FunctionTok{grep}\NormalTok{(}\StringTok{"log.per.cap.income"}\NormalTok{,}\FunctionTok{names}\NormalTok{(cdilogred.cont))}

\NormalTok{y }\OtherTok{\textless{}{-}}\NormalTok{ cdilogred.cont[,loc]}

\NormalTok{X }\OtherTok{\textless{}{-}} \FunctionTok{apply}\NormalTok{(}\FunctionTok{as.matrix}\NormalTok{(cdilogred.cont[,}\SpecialCharTok{{-}}\NormalTok{loc]),}\DecValTok{2}\NormalTok{,}\ControlFlowTok{function}\NormalTok{(x) }\FunctionTok{rescale}\NormalTok{(x,}\StringTok{"full"}\NormalTok{))}

\NormalTok{Xnames }\OtherTok{\textless{}{-}} \FunctionTok{dimnames}\NormalTok{(X)[[}\DecValTok{2}\NormalTok{]]}

\NormalTok{lasso.result }\OtherTok{\textless{}{-}} \FunctionTok{glmnet}\NormalTok{(X,y)}

\FunctionTok{plot}\NormalTok{(lasso.result,}\AttributeTok{xvar=}\StringTok{"lambda"}\NormalTok{,}\AttributeTok{xlim=}\FunctionTok{c}\NormalTok{(}\SpecialCharTok{{-}}\DecValTok{9}\NormalTok{,}\DecValTok{0}\NormalTok{))}

\FunctionTok{abline}\NormalTok{(}\AttributeTok{h=}\DecValTok{0}\NormalTok{,}\AttributeTok{lty=}\DecValTok{2}\NormalTok{)}

\FunctionTok{legend}\NormalTok{(}\StringTok{\textquotesingle{}bottomright\textquotesingle{}}\NormalTok{,}\AttributeTok{lty=}\DecValTok{1}\NormalTok{,}\AttributeTok{col=}\DecValTok{1}\SpecialCharTok{:}\FunctionTok{length}\NormalTok{(Xnames),}\AttributeTok{legend=}\NormalTok{Xnames,}\AttributeTok{cex=}\FloatTok{0.75}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-51-1.pdf}
There is not such an obvious place to cut the shrinkage plot, to help determine what variables should be kept in the model.

Let's see if choosing lambda by cross-valiation works:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{cv.lasso.result }\OtherTok{\textless{}{-}} \FunctionTok{cv.glmnet}\NormalTok{(X,y)}

\FunctionTok{plot}\NormalTok{(cv.lasso.result)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-52-1.pdf}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{c}\NormalTok{(}\AttributeTok{lambda.1se=}\NormalTok{cv.lasso.result}\SpecialCharTok{$}\NormalTok{lambda}\FloatTok{.1}\NormalTok{se,}\AttributeTok{lambda.min=}\NormalTok{cv.lasso.result}\SpecialCharTok{$}\NormalTok{lambda.min)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##   lambda.1se   lambda.min 
## 0.0049081710 0.0004795332
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{tmp }\OtherTok{\textless{}{-}} \FunctionTok{cbind}\NormalTok{(}\FunctionTok{coef}\NormalTok{(cv.lasso.result,}\AttributeTok{s=}\NormalTok{cv.lasso.result}\SpecialCharTok{$}\NormalTok{lambda.min),}
             \FunctionTok{coef}\NormalTok{(cv.lasso.result,}\AttributeTok{s=}\NormalTok{cv.lasso.result}\SpecialCharTok{$}\NormalTok{lambda}\FloatTok{.1}\NormalTok{se)}
\NormalTok{            )}
\FunctionTok{dimnames}\NormalTok{(tmp)[[}\DecValTok{2}\NormalTok{]] }\OtherTok{\textless{}{-}} \FunctionTok{c}\NormalTok{(}\StringTok{"lambda(minMSE)"}\NormalTok{,}\StringTok{"lambda(minMSE+1se)"}\NormalTok{)}

\NormalTok{tmp}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 11 x 2 sparse Matrix of class "dgCMatrix"
##               lambda(minMSE) lambda(minMSE+1se)
## (Intercept)     9.806955e+00        9.806954586
## log.land.area  -6.233501e-02       -0.058650987
## pop.18_34      -1.266292e-01       -0.103473858
## pop.65_plus    -2.099172e-02        .          
## log.doctors     1.178374e-01        0.137118673
## log.hosp.beds   2.479267e-02        .          
## log.crimes     -4.850651e-05        .          
## pct.hs.grad    -5.975737e-02       -0.005665686
## pct.bach.deg    2.352630e-01        0.187649870
## pct.below.pov  -2.281218e-01       -0.191409203
## pct.unemp       4.988665e-02        0.034590648
\end{verbatim}

It's interesting to note that the model which minimizes 10-fold cross-validation error contains all 10 predictors (or sometimes 9---log.crimes can get left out---depending on what random folds we get), while the model that is 1 SE above that is again the same model that we first saw with all-subsets regression.

Although the magnitudes of the coefficients are a little different from what we saw before (because they arise from mimimizing RSS penalized by the L1 penalty, rather than minimizing straight RSS), the signs are the same.

Once again, if we were to settle on this model, and try to bring in \texttt{region}, chances are we would end up with the same model as \texttt{all.subsets.01.final.with.some.region}.

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{round}\NormalTok{(}\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}\SpecialCharTok{$}\NormalTok{coef,}\DecValTok{2}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)               10.13       0.25   40.25     0.00
## log.land.area             -0.03       0.01   -6.43     0.00
## pop.18_34                 -0.02       0.00  -14.17     0.00
## log.doctors                0.06       0.00   13.72     0.00
## pct.hs.grad                0.00       0.00   -0.79     0.43
## pct.bach.deg               0.01       0.00    6.74     0.00
## pct.below.pov             -0.02       0.00   -7.05     0.00
## pct.unemp                  0.02       0.00    3.72     0.00
## regionNE                   0.22       0.30    0.73     0.47
## regionS                   -0.06       0.28   -0.23     0.82
## regionW                    1.63       0.36    4.56     0.00
## pct.hs.grad:regionNE       0.00       0.00   -0.94     0.35
## pct.hs.grad:regionS        0.00       0.00    0.57     0.57
## pct.hs.grad:regionW       -0.02       0.00   -4.50     0.00
## pct.bach.deg:regionNE      0.01       0.00    2.26     0.02
## pct.bach.deg:regionS       0.00       0.00   -0.56     0.58
## pct.bach.deg:regionW       0.01       0.00    2.41     0.02
## pct.below.pov:regionNE     0.00       0.00   -0.52     0.60
## pct.below.pov:regionS      0.01       0.00    1.94     0.05
## pct.below.pov:regionW     -0.02       0.01   -2.95     0.00
## pct.unemp:regionNE        -0.01       0.01   -1.10     0.27
## pct.unemp:regionS         -0.02       0.01   -3.87     0.00
## pct.unemp:regionW         -0.02       0.01   -2.99     0.00
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{R2 = "}\NormalTok{,}\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}\SpecialCharTok{$}\NormalTok{r.squared)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## R2 =  0.8734897
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{cat}\NormalTok{(}\StringTok{"}\SpecialCharTok{\textbackslash{}n}\StringTok{R2adj = "}\NormalTok{,}\FunctionTok{summary}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}\SpecialCharTok{$}\NormalTok{adj.r.squared)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## R2adj =  0.8668153
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{par}\NormalTok{(}\AttributeTok{mfrow=}\FunctionTok{c}\NormalTok{(}\DecValTok{2}\NormalTok{,}\DecValTok{2}\NormalTok{))}
\FunctionTok{plot}\NormalTok{(all.subsets.}\FloatTok{01.}\NormalTok{final.with.some.region)}
\end{Highlighting}
\end{Shaded}

\includegraphics{hw06_solutions_files/figure-latex/unnamed-chunk-53-1.pdf}

Here are a few of the models plus's and minus's, and some tradeoffs I considered:

\begin{itemize}
\tightlist
\item
  \textbf{Plus's:}

  \begin{itemize}
  \tightlist
  \item
    The model is fairly parsimonious, and most of the estimated coefficients have the expected sign.
  \item
    The model is confirmed by stepwise and lasso procedures.\\
  \item
    Those procedures also found more complex models with somewhat better fit, but improvements in fit seemed small compared to the added complexity of the model.
  \item
    All of the variables are either in their original scale, or the have been replaced with their logarithm. This facilitates explaining the models to anyone who is interested in and knowledgeable about the social science \& economics but less knowledgeable about technical matters.
  \end{itemize}
\item
  \textbf{Minus's:}

  \begin{itemize}
  \tightlist
  \item
    The coefficient on \texttt{pct.unemp} seems to go the wrong way, and the coefficient on 'pct.hs.grad` is quite small, statistically and practically (it remains in the model because there is a noticable interaction that it participates in).
  \item
    The residual diagnostic plots are just OK.
  \item
    The fact that stepwise regression found some well-fitting models with interactions between continuous variables suggests exploring those more complex models in the future.
  \item
    I did not have time to explore the state variable at all. Some of the relationship between these demographic variables and per capita income might be explainable in terms of varying economic policy from one state to the next. (If one includes \texttt{state} in the model, one could take out \texttt{region} because the two are perfectly collinear (states are entirely nested within regions).
  \end{itemize}
\end{itemize}

Finally, it would be very useful to have additional data to compare some of the models we found. We are using reasonable methods for variable selection, but since it is all within-sample (our entire data set is our training sample), there is ample room for overfitting noise in the data. Some of our inferences about which variables to leave in or take out may be based on overly optimistic standard error estimates, for example. If we were able to cross-validate on some new or hold-out data, we might be able to better distinguish what is the best model, at least in terms of prediction error.

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

This will be much quicker and will mostly consist of code and comments written by a TA from an earlier version of this course.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{beauty }\OtherTok{\textless{}{-}} \FunctionTok{read\_csv}\NormalTok{(}\StringTok{"../ProfEvaltnsBeautyPublic.csv"}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  dplyr}\SpecialCharTok{::}\FunctionTok{select}\NormalTok{(}\SpecialCharTok{{-}}\NormalTok{profevaluation, }\SpecialCharTok{{-}}\NormalTok{profnumber, }\SpecialCharTok{{-}}\NormalTok{multipleclass, }
                \SpecialCharTok{{-}}\FunctionTok{starts\_with}\NormalTok{(}\StringTok{"class"}\NormalTok{))}
\end{Highlighting}
\end{Shaded}

\hypertarget{part-a-1}{%
\subsection{Part a}\label{part-a-1}}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(leaps)}
\NormalTok{all\_subsets }\OtherTok{\textless{}{-}} \FunctionTok{regsubsets}\NormalTok{(courseevaluation }\SpecialCharTok{\textasciitilde{}}\NormalTok{.,}
                          \AttributeTok{data =}\NormalTok{ beauty,}
                          \AttributeTok{nvmax =} \DecValTok{10}\NormalTok{, }\AttributeTok{intercept =}\NormalTok{ T,}
                          \AttributeTok{method =} \StringTok{"exhaustive"}\NormalTok{)}

\NormalTok{summ }\OtherTok{\textless{}{-}} \FunctionTok{summary}\NormalTok{(all\_subsets)}

\NormalTok{best\_model\_index }\OtherTok{\textless{}{-}} \FunctionTok{which.min}\NormalTok{(summ}\SpecialCharTok{$}\NormalTok{bic)}

\NormalTok{dropped }\OtherTok{\textless{}{-}}\NormalTok{ summ}\SpecialCharTok{$}\NormalTok{outmat }\SpecialCharTok{\%\textgreater{}\%} \FunctionTok{apply}\NormalTok{(}\DecValTok{2}\NormalTok{, }\ControlFlowTok{function}\NormalTok{(col) \{}\FunctionTok{mean}\NormalTok{(col }\SpecialCharTok{==} \StringTok{" "}\NormalTok{) }\SpecialCharTok{==} \DecValTok{1}\NormalTok{\})}
\NormalTok{all\_subset\_info }\OtherTok{\textless{}{-}} \FunctionTok{cbind}\NormalTok{(summ}\SpecialCharTok{$}\NormalTok{outmat[,}\SpecialCharTok{!}\NormalTok{dropped], }\AttributeTok{BIC =}\NormalTok{ summ}\SpecialCharTok{$}\NormalTok{bic) }
\FunctionTok{rownames}\NormalTok{(all\_subset\_info) }\OtherTok{\textless{}{-}} \DecValTok{1}\SpecialCharTok{:}\FunctionTok{nrow}\NormalTok{(all\_subset\_info)}
\NormalTok{all\_subset\_info }\SpecialCharTok{\%\textgreater{}\%}\NormalTok{ data.frame }\SpecialCharTok{\%\textgreater{}\%}
  \FunctionTok{rownames\_to\_column}\NormalTok{(}\AttributeTok{var =} \StringTok{"num variables"}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  knitr}\SpecialCharTok{::}\FunctionTok{kable}\NormalTok{() }\SpecialCharTok{\%\textgreater{}\%} 
\NormalTok{  kableExtra}\SpecialCharTok{::}\FunctionTok{kable\_styling}\NormalTok{(}\AttributeTok{latex\_options =} \StringTok{"scale\_down"}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%}
  \FunctionTok{row\_spec}\NormalTok{(best\_model\_index, }\AttributeTok{background =} \StringTok{"gray"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{table}
\centering
\resizebox{\linewidth}{!}{
\begin{tabular}{l|l|l|l|l|l|l|l|l|l|l|l|l}
\hline
num variables & minority & beautyf2upper & beautyfupperdiv & btystdave & btystdf2u & btystdfu & female & nonenglish & onecredit & percentevaluating & blkandwhite & BIC\\
\hline
1 &  &  &  &  &  &  &  &  & * &  &  & -14.0820700135937\\
\hline
2 &  &  &  & * &  &  &  &  & * &  &  & -30.0187608382119\\
\hline
3 &  &  & * &  &  &  & * &  & * &  &  & -42.0422726206511\\
\hline
4 &  &  & * &  &  &  & * &  & * & * &  & -52.6352707419912\\
\hline
5 & * &  & * &  &  &  & * &  & * & * &  & -61.8604866501093\\
\hline
\cellcolor{gray}{6} & \cellcolor{gray}{} & \cellcolor{gray}{} & \cellcolor{gray}{*} & \cellcolor{gray}{} & \cellcolor{gray}{} & \cellcolor{gray}{} & \cellcolor{gray}{*} & \cellcolor{gray}{*} & \cellcolor{gray}{*} & \cellcolor{gray}{*} & \cellcolor{gray}{*} & \cellcolor{gray}{-68.7273809819797}\\
\hline
7 & * &  & * &  &  &  & * & * & * & * & * & -68.168278525645\\
\hline
8 &  & * & * &  & * &  & * & * & * & * & * & -67.592998887792\\
\hline
9 & * & * & * &  & * &  & * & * & * & * & * & -67.4303455606005\\
\hline
10 & * & * & * &  & * & * & * & * & * & * & * & -64.9400133420853\\
\hline
\end{tabular}}
\end{table}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{var\_names }\OtherTok{\textless{}{-}} \FunctionTok{colnames}\NormalTok{(all\_subset\_info)[all\_subset\_info[best\_model\_index,] }\SpecialCharTok{==} \StringTok{"*"}\NormalTok{]}

\NormalTok{all\_subsets\_model }\OtherTok{\textless{}{-}} \FunctionTok{lm}\NormalTok{(}\FunctionTok{as.formula}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"courseevaluation \textasciitilde{}"}\NormalTok{,}
                                         \FunctionTok{paste0}\NormalTok{(var\_names, }\AttributeTok{collapse =} \StringTok{" + "}\NormalTok{))),}
                        \AttributeTok{data =}\NormalTok{ beauty)}

\FunctionTok{summary}\NormalTok{(all\_subsets\_model)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = as.formula(paste("courseevaluation ~", paste0(var_names, 
##     collapse = " + "))), data = beauty)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.97607 -0.30533  0.04791  0.37959  1.04118 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.355273   0.113922  29.452  < 2e-16 ***
## beautyfupperdiv    0.058282   0.012559   4.641 4.54e-06 ***
## female            -0.254260   0.048290  -5.265 2.16e-07 ***
## nonenglish        -0.388186   0.097125  -3.997 7.49e-05 ***
## onecredit          0.500846   0.099146   5.052 6.35e-07 ***
## percentevaluating  0.005502   0.001415   3.889 0.000115 ***
## blkandwhite        0.251170   0.063629   3.947 9.15e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.495 on 456 degrees of freedom
## Multiple R-squared:  0.2143, Adjusted R-squared:  0.204 
## F-statistic: 20.73 on 6 and 456 DF,  p-value: < 2.2e-16
\end{verbatim}

\hypertarget{part-b-1}{%
\subsection{Part b}\label{part-b-1}}

We'll show forward regression.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{forward\_step }\OtherTok{\textless{}{-}} \FunctionTok{step}\NormalTok{(}\FunctionTok{lm}\NormalTok{(courseevaluation }\SpecialCharTok{\textasciitilde{}}\DecValTok{1}\NormalTok{, }\AttributeTok{data =}\NormalTok{ beauty),}
  \AttributeTok{scope  =} \FunctionTok{as.formula}\NormalTok{(}\FunctionTok{paste}\NormalTok{(}\StringTok{"courseevaluation \textasciitilde{}"}\NormalTok{,}
                      \FunctionTok{paste0}\NormalTok{(}\FunctionTok{names}\NormalTok{(beauty)[}\FunctionTok{names}\NormalTok{(beauty) }\SpecialCharTok{!=} \StringTok{"courseevaluation"}\NormalTok{], }
                             \AttributeTok{collapse =} \StringTok{" + "}\NormalTok{))),}
                          \AttributeTok{direction =} \StringTok{"forward"}\NormalTok{, }\AttributeTok{trace =} \DecValTok{0}\NormalTok{)}

\FunctionTok{summary}\NormalTok{(forward\_step)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## lm(formula = courseevaluation ~ onecredit + btystdave + percentevaluating + 
##     female + minority + blkandwhite + nonenglish + beautyfupperdiv + 
##     btystdfu + btystdf2u + tenuretrack + formal + age, data = beauty)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.86588 -0.29589  0.04495  0.33963  1.06133 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -0.836838   1.639295  -0.510 0.609962    
## onecredit          0.503749   0.107838   4.671 3.96e-06 ***
## btystdave         -0.167839   0.075372  -2.227 0.026456 *  
## percentevaluating  0.004850   0.001451   3.342 0.000902 ***
## female            -0.302259   0.052308  -5.778 1.41e-08 ***
## minority          -0.155301   0.074788  -2.077 0.038412 *  
## blkandwhite        0.285158   0.069381   4.110 4.70e-05 ***
## nonenglish        -0.381223   0.104535  -3.647 0.000297 ***
## beautyfupperdiv    0.939194   0.323894   2.900 0.003918 ** 
## btystdfu          -1.770340   0.656800  -2.695 0.007294 ** 
## btystdf2u          0.092159   0.042790   2.154 0.031790 *  
## tenuretrack       -0.132773   0.062020  -2.141 0.032828 *  
## formal             0.143150   0.069643   2.055 0.040410 *  
## age               -0.004442   0.002837  -1.566 0.118145    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4867 on 449 degrees of freedom
## Multiple R-squared:  0.2523, Adjusted R-squared:  0.2306 
## F-statistic: 11.65 on 13 and 449 DF,  p-value: < 2.2e-16
\end{verbatim}

\hypertarget{part-c-1}{%
\subsection{Part c}\label{part-c-1}}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{library}\NormalTok{(glmnet)}
\NormalTok{cv\_lasso\_selection }\OtherTok{\textless{}{-}} \FunctionTok{cv.glmnet}\NormalTok{(}\AttributeTok{x =} \FunctionTok{as.matrix}\NormalTok{(beauty[,}\FunctionTok{names}\NormalTok{(beauty) }\SpecialCharTok{!=} \StringTok{"courseevaluation"}\NormalTok{]),}
          \AttributeTok{y =}\NormalTok{ beauty}\SpecialCharTok{$}\NormalTok{courseevaluation, }\AttributeTok{alpha =} \DecValTok{1}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{lasso }\OtherTok{\textless{}{-}} \FunctionTok{glmnet}\NormalTok{(}\AttributeTok{x =} \FunctionTok{as.matrix}\NormalTok{(beauty[,}\FunctionTok{names}\NormalTok{(beauty) }\SpecialCharTok{!=} \StringTok{"courseevaluation"}\NormalTok{]),}
                \AttributeTok{y =}\NormalTok{ beauty}\SpecialCharTok{$}\NormalTok{courseevaluation,}
                \AttributeTok{alpha =} \DecValTok{1}\NormalTok{, }\AttributeTok{lambda =}\NormalTok{ cv\_lasso\_selection}\SpecialCharTok{$}\NormalTok{lambda}\FloatTok{.1}\NormalTok{se)}


\NormalTok{coef\_lasso }\OtherTok{\textless{}{-}} \FunctionTok{as.matrix}\NormalTok{(}\FunctionTok{coef}\NormalTok{(lasso))}
\NormalTok{coef\_lasso\_df }\OtherTok{\textless{}{-}} \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{columns =} \FunctionTok{names}\NormalTok{(coef\_lasso[coef\_lasso[,}\DecValTok{1}\NormalTok{] }\SpecialCharTok{!=} \DecValTok{0}\NormalTok{,]),}
         \AttributeTok{beta =}\NormalTok{ coef\_lasso[coef\_lasso[,}\DecValTok{1}\NormalTok{] }\SpecialCharTok{!=} \DecValTok{0}\NormalTok{,]) }
\FunctionTok{rownames}\NormalTok{(coef\_lasso\_df) }\OtherTok{\textless{}{-}} \ConstantTok{NULL}

\NormalTok{coef\_lasso\_df}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##              columns          beta
## 1        (Intercept)  3.533028e+00
## 2           minority -4.888278e-02
## 3      beautyf2upper  5.133109e-03
## 4    beautyfupperdiv  2.973839e-02
## 5          btystdf2u  3.833213e-06
## 6             female -1.166587e-01
## 7           fulldept  7.914660e-02
## 8         nonenglish -1.675603e-01
## 9          onecredit  3.553839e-01
## 10 percentevaluating  3.304740e-03
## 11       tenuretrack -2.928520e-03
## 12       blkandwhite  1.190623e-01
\end{verbatim}

\hypertarget{part-d}{%
\subsection{Part d}\label{part-d}}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{\#basically coding up a document term matrix function for some reason}

\NormalTok{dtm\_coef }\OtherTok{\textless{}{-}} \ControlFlowTok{function}\NormalTok{(names\_list)\{}
\NormalTok{  all\_words }\OtherTok{\textless{}{-}} \FunctionTok{unique}\NormalTok{(}\FunctionTok{unlist}\NormalTok{(names\_list))}
\NormalTok{  matrix\_info }\OtherTok{\textless{}{-}} \FunctionTok{matrix}\NormalTok{(}\AttributeTok{nrow =} \FunctionTok{length}\NormalTok{(names\_list), }\AttributeTok{ncol =} \FunctionTok{length}\NormalTok{(all\_words))}
  
\NormalTok{  r\_idx }\OtherTok{\textless{}{-}} \DecValTok{1}
  \ControlFlowTok{for}\NormalTok{ (col\_names }\ControlFlowTok{in}\NormalTok{ names\_list)\{}
\NormalTok{    matrix\_info[r\_idx,] }\OtherTok{\textless{}{-}}\NormalTok{ all\_words }\SpecialCharTok{\%in\%}\NormalTok{ col\_names}
\NormalTok{    r\_idx }\OtherTok{\textless{}{-}}\NormalTok{ r\_idx }\SpecialCharTok{+} \DecValTok{1}
\NormalTok{  \}}
  
  \FunctionTok{colnames}\NormalTok{(matrix\_info) }\OtherTok{\textless{}{-}}\NormalTok{ all\_words}
  \FunctionTok{rownames}\NormalTok{(matrix\_info) }\OtherTok{\textless{}{-}} \FunctionTok{names}\NormalTok{(names\_list)}
  
  \FunctionTok{return}\NormalTok{(matrix\_info)}
\NormalTok{\}}

\NormalTok{variables\_in\_model }\OtherTok{\textless{}{-}} \FunctionTok{list}\NormalTok{(}
  \StringTok{"best subset"} \OtherTok{=} \FunctionTok{names}\NormalTok{(}\FunctionTok{coef}\NormalTok{(all\_subsets\_model)),}
  \StringTok{"lasso"} \OtherTok{=}\NormalTok{ coef\_lasso\_df}\SpecialCharTok{$}\NormalTok{columns,}
  \StringTok{"forward stepwise"} \OtherTok{=} \FunctionTok{names}\NormalTok{(forward\_step}\SpecialCharTok{$}\NormalTok{coefficients),}
  \StringTok{"homework 05, 2c"} \OtherTok{=} \FunctionTok{c}\NormalTok{(}\FunctionTok{names}\NormalTok{(beauty[,}\FunctionTok{names}\NormalTok{(beauty) }\SpecialCharTok{!=} \StringTok{"courseevaluation"}\NormalTok{]), }
                        \StringTok{"profevaluation"}\NormalTok{, }\StringTok{"(Intercept)"}\NormalTok{)}
\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{Shaded}
\begin{Highlighting}[]
\FunctionTok{t}\NormalTok{(}\FunctionTok{dtm\_coef}\NormalTok{(variables\_in\_model)) }\SpecialCharTok{\%\textgreater{}\%}
  \FunctionTok{data.frame}\NormalTok{(}\AttributeTok{check.names =}\NormalTok{ F) }\SpecialCharTok{\%\textgreater{}\%} 
  \FunctionTok{rownames\_to\_column}\NormalTok{(}\AttributeTok{var =} \StringTok{"features"}\NormalTok{) }\SpecialCharTok{\%\textgreater{}\%}
  \FunctionTok{mutate}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{best subset}\StringTok{\textasciigrave{}} \OtherTok{=} \FunctionTok{ifelse}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{best subset}\StringTok{\textasciigrave{}}\NormalTok{,  }
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{best subset}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"black"}\NormalTok{),}
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{best subset}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"red"}\NormalTok{)),}
         \StringTok{\textasciigrave{}}\AttributeTok{lasso}\StringTok{\textasciigrave{}} \OtherTok{=} \FunctionTok{ifelse}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{lasso}\StringTok{\textasciigrave{}}\NormalTok{,  }
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{lasso}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"black"}\NormalTok{),}
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{lasso}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"red"}\NormalTok{)),}
         \StringTok{\textasciigrave{}}\AttributeTok{forward stepwise}\StringTok{\textasciigrave{}} \OtherTok{=} \FunctionTok{ifelse}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{forward stepwise}\StringTok{\textasciigrave{}}\NormalTok{,  }
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{forward stepwise}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"black"}\NormalTok{),}
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{forward stepwise}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"red"}\NormalTok{)),}
         \StringTok{\textasciigrave{}}\AttributeTok{homework 05, 2c}\StringTok{\textasciigrave{}} \OtherTok{=} \FunctionTok{ifelse}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{homework 05, 2c}\StringTok{\textasciigrave{}}\NormalTok{,  }
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{homework 05, 2c}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"black"}\NormalTok{),}
                                \FunctionTok{cell\_spec}\NormalTok{(}\StringTok{\textasciigrave{}}\AttributeTok{homework 05, 2c}\StringTok{\textasciigrave{}}\NormalTok{, }\StringTok{"latex"}\NormalTok{, }\AttributeTok{color =} \StringTok{"red"}\NormalTok{))) }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  knitr}\SpecialCharTok{::}\FunctionTok{kable}\NormalTok{(}\AttributeTok{format =} \StringTok{"latex"}\NormalTok{, }\AttributeTok{caption =} \StringTok{"Features include in each model."}\NormalTok{,}
               \AttributeTok{escape =}\NormalTok{ F) }\SpecialCharTok{\%\textgreater{}\%}
\NormalTok{  kableExtra}\SpecialCharTok{::}\FunctionTok{kable\_styling}\NormalTok{(}\AttributeTok{latex\_options =} \FunctionTok{c}\NormalTok{(}\StringTok{"scale\_down"}\NormalTok{),)}
\end{Highlighting}
\end{Shaded}

\begin{table}

\caption{\label{tab:tab-feature}Features include in each model.}
\centering
\resizebox{\linewidth}{!}{
\begin{tabular}[t]{l|l|l|l|l}
\hline
features & best subset & lasso & forward stepwise & homework 05, 2c\\
\hline
(Intercept) & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
beautyfupperdiv & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
female & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
nonenglish & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
onecredit & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
percentevaluating & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
blkandwhite & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
minority & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
beautyf2upper & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdf2u & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
fulldept & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
tenuretrack & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
btystdave & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
btystdfu & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
formal & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
age & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE} & \textcolor{black}{TRUE}\\
\hline
tenured & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
beautyflowerdiv & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
beautym2upper & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
beautymlowerdiv & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
beautymupperdiv & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdfl & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdm2u & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdml & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdmu & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
didevaluation & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
lower & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
students & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdvariance & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdavepos & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
btystdaveneg & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
profevaluation & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{red}{FALSE} & \textcolor{black}{TRUE}\\
\hline
\end{tabular}}
\end{table}

The number of features (excluding the intercept) in our models increase from 6 in the best subset model, lasso with 10, forward stepwise regression selecting 14 and our full model with 30 features in homework 05. If we take out the lasso model these models are nested inside each other. \textbf{It's hard to evaluate prediction models as we left nothing to evaluate them on.} We observe that the lasso had actually higher training mse than the best subset model - but that's not to crazy because the penalty means that we're get more training error allowed.

Following the advice of Gelman and Hill from the end of lecture 11, it might also be useful to compare tables of estimated coefficients and their SE's, to see whether the direction and size of influence of each predictor on teacher ratings is (a) plausible and (b) consistent (or inconsistent) across models. As usual, however, since we are selecting variables and fitting the model on the same data, we can be misled by overly optimistic standard errors, overfitting to noise in the data, etc. So this approach ultimately requires some care.

\end{document}
