\documentstyle[bj,txfonts]{article}

\frenchspacing
\def\vtrim{\vspace*{-0.8em}}

\begin{document}

\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
HW05 -- Due Weds Oct 12, 11:59pm 
\end{center}

\begin{itemize}

\item Please turn the homework in online to gradescope using the link on the assignment page in canvas.

\item Reading for this week:
\begin{itemize}
\item Sheather, appendix on nonparametric regression
\item ISLR, Ch 7
\end{itemize}

\item You will need the file nes5200\_{}processed\_{}voters\_{}realideo.dta for problem \#2 below.  It is in the folder for this assignment on Canvas. 

\item You will also need to install three R packages for this assignment:
\begin{itemize}
\item {\tt marginalmodelplots.} We will use the {\tt mmplots()} function from {\tt library(marginalmodelplots)} rather than the  {\tt mmps()} function from {\tt car}, since {\tt mmps()} doesn't work right for glm's.  To install it:
\begin{enumerate}
\item download \verb|marginalmodelplots_0.4.2.tar.gz| from canvas (a copy is in the directory for this hw assignment)
\item Instal with these commands:
\begin{verbatim}
install.packages("locfit")
install.packages("marginalmodelplots_0.4.2.tar.gz", repos=NULL)
\end{verbatim}
\end{enumerate}  
\item {\tt DHARMa.}  This provides an alternative set of residuals for all lm's and glm's based on a simulation method called the ``parametric bootstrap''. Install from {\tt cran} as usual.  Here is an example code that uses {\tt DHARMa}:
\begin{verbatim}
## Simulate some data (if you have real data you don't have to do this!)
expit <- function(x) {exp(x)/(1+exp(x))}
x <- rnorm(100)
p <- expit(2 -3*x)
y <- rbinom(100,1,p)
## End of simulation (if you have real data use that in the glm() function below!)

## fit a logistic glm
mymodel <- glm(y ~ x, family=binomial)
## use DHARMa to check residuals
simulationOutput <- simulateResiduals(fittedModel = mymodel, plot = T)
\end{verbatim}
This produces a qq plot on the left and a residuals vs fitted plot on the right.  
\begin{itemize}
\item If the model fits the data, DHARMa produces uniformly distributed residuals, instead of normally distributed residuals.  You can use the qq plot to check for how close to uniform the residuals are, check for skewing, overdispersion, etc., just like the normal qq plot from the usual residual diagnostic plots in R.
\item The plot on the right is a residuals vs fitted plot that you can interpret much as for ordinary regression.  The three guidelines help you assess nonconstant varianace.  If the lines are horizontal, the variance is consistent with the model (e.g., constant variance for ordinary regression; variance proportional to $p_i(1-p_i)$ for logistic regression).  Residuals that are outliers will be colored red in the plot.
\end{itemize}
There is much more you can do with DHARMa; if you are curious, have a look at \\ {\tt https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html} 
\item {\tt arm.} This is the library that goes with the Gelman and Hill text, and provides the {\\t binnedplot()} function, which lets you look at local averages of residuals, as in ther lecture notes.
\end{itemize}


\end{itemize}


\section*{Exercises}

% \begin{enumerate}% \itemsep=0pt\parsep=0pt\parskip=3pt

% \item 

\begin{enumerate}\setcounter{enumi}{-1}

\item Unlike the midterm, it is OK to collaborate on hw's.  Please list who you worked with on this hw.  If you didn't work with anyone, write ``I worked alone''.

\comment{

\item Let $X$ be a continuous random variable with cumulative distribution function (CDF) $F(x) = P[X\leq x]$, and assume that the inverse function $F^{-1}(u)$ exists.  Now let $Y$ be the transformed variable $Y=F(X)$, and let $G(y) = P[Y\leq y]$.  Show that 
\[
     G(y) = \left\{\begin{array}{rl}
             0 & \mbox{, if $y\leq 0$}\\
             y & \mbox{, if $0<y<1$} \\
             1 & \mbox{, if $1\leq y$}
             \end{array}\right.
\]  
i.e., $Y\sim Unif(0,1)$.  This is a key idea behind the DHARMa residuals approach.  ({\em Hint: No calculus required.  Just write down $G(y) = P[Y\leq y]$, substitute in the definition of $Y$, and follow your nose!})

}

\item Please do ISLR, \#11, pp. 325--326. ({\em Hint: You have seen me do a variation of part (a) many times in class, for demonstrations! Some examples are in the lecture notes on omitted variable bias, for example.}) Among other things, this shows that backfitting (a method used to estimate GAM's) is similar to iteratively adjusting added variable plots(!)\ldots

\item The {\tt warpbreaks} data set included in R gives the results of an experiment to determine the the effect of wool type (A or B) and tension (low, medium or high) on the number of warp breaks per loom. Data was collected for nine looms for each combination of settings. You can get more information on this data set from {\tt help(warpbreaks)}.  Use {\tt View(warpbreaks)} and \verb|xtabs(breaks ~ wool + tension, data=warpbreaks)| to familiarize yourself with the data.  (You don't have to turn in any of this initial exploration.)
\begin{enumerate}
\item Fit the two Poisson regression models \verb|breaks ~ wool + tension| and \verb|breaks ~ wool * tension| (dont't forget {\tt family=poisson}!).  Consider {\tt summary()}'s, {\tt plot()}'s of the fitted glm's, {\tt mmplot()}'s and also {\tt binnedplot()}'s (from {\tt library(arm)}) and/or plots from {\tt DHARMa}, and any other plots you think are useful.  Comment on the fits of each model, using the graphical evidence you have obtained {\em (if you think some plot(s) are \ul{\bf not} useful for this, please tell which ones \& why)}.

\item We would like to know whether the interaction between {\tt wool} and {\tt tension} is needed.  
\begin{enumerate}
\item Compare AIC's for the two models, and then compare BIC's for the two models.  Do AIC and BIC choose the same model? 
\item Conduct a likelihood ratio test for including the interaction. Does the LR test choose the same model as either AIC or BIC? {\em (Hint/recipe: For the likelihood ratio test, use the {\tt logLik()} function to obtain the log-likelihood of the two models,  compute $-2$ times the difference between the null and alternative model log-liklihoods, and obtain a $p$-value from the tail of an approporiate $\chi^2$ distribution.)} 
\end{enumerate}
\item For the best model from part (b), we would like to know if overdispersion is a problem. Use the Pearson residuals to conduct a test of overdispersion (hint: follow the recipe from lecture \#17.1).  Does this agree with the conclusion of {\tt testDispersion()} from {\tt library(DHARMa)}?   
\item Refit the model from part (c) using {\tt family=quasipoisson}, and note any differences in overdispersion parameter, coefficient estimates, SE's, etc.  Are there any differences in which predictors are significant?

\end{enumerate}

% \vspace*{\fill}
%
% \begin{center}
% {\sl [More exercises on the next page]}
% \end{center}
%
% \clearpage

\item{} [Based on some exercises in Gelman \& Hill.] The data set {\tt nes5200\_{}processed\_{}voters\_{}realideo.dta} has data from the National Election Survey in every congressional election year in the US from 1948 to 2002.  You can read it into R using the function {\tt read.dta()} from {\tt library(foreign)}.  We are interested in learning whether some of the survey data is useful in predicting respondents' vote for president.

\begin{enumerate}

\item\label{extract} Extract from the full NES data set above a data frame with just data from the presidential election year 1976, and the following variables: {\tt ideo7, black, female, rep\_{}pres\_{}intent, presapprov, age, educ1, urban, income, union} and {\tt perfin1}.  The variable {\tt rep\_{}pres\_{}intent} is 1 if the respondent intended to vote for the Republican candidate (Gerald Ford, who assumed the presidency when Richard Nixon resigned amid an impeachment investigation), and 0 if the respondent intended to vote for the Democratic candidate (Jimmy Carter).  Familiarize yourself with all of the variables by using {\tt summary()} on your 1976 data frame\footnote{A command like {\tt attr(nes\_{}data,"var.labels")} on the full data set will give you brief definitions of all of the variables.} to look at the categories. 
You don't have to turn anything in for this, but you should do it.

\item\label{reg1} Fit a logistic regression predicting {\tt rep\_{}pres\_{}intent} from the other variables, for 1976.  Assess the fit using appropriate graphical methods.  Interpret the results: which variables seem to matter for predicting voters' preference for president?  How do they affect the odds of voting for the Republican candidate?

\item\label{reg2} Make a new data frame by converting all of the variables in the 1976 data frame to numeric variables, using {\tt as.numeric()}.  Repeat part (\ref{reg1}), using this new data frame.

\item\label{compare} Compare the results from parts (\ref{reg1}) and (\ref{reg2}).  How are they consistent? Inconsistent?  In this case, would you recommend converting all of the variables to numeric?

\item NOT TO TURN IN.  If you have time, it is interesting to repeat parts (\ref{extract}) through (\ref{compare}) for the presidential election year 2000.  In this year, the democratic candidate was Al Gore, who served as vice-president to Bill Clinton, who was impeached but not removed from office.  The Republican candidate was George Bush.  How are the results from 1976 and 2000 similar?  How are they different?

%\item Repeat parts (\ref{extract}) through (\ref{compare}) for the presidential election year 2000.  In this year, the democratic candidate was Al Gore, who served as vice-president to Bill Clinton, who was impeached but not removed from office.  The Republican candidate was George Bush.

% \item Finally, compare your results from 1976 and 2000.  Interpret any similarities or differences in terms of voter preferences.  

\end{enumerate}
Of course there are very many other variables, as well as transformations and interactions, to explore in this data set, and the models you fitted above are not ``optimal'' in any sense.  However I hope you see that you can still get interesting signal out of the data, even with non-optimal models.

% it would be nice to return to the nes data, have the students 
% add ``presvote_intent'' to their 1976 and 2000 data sets, 
% and for each data set have them:
%
% library(nnet)
% summary(mu.3 <- multinom(presvote_intent ~ black + female + age,data=nes1976c))
% 
% then
% (a) interpret this model
% (b) try to add more variables from the 1976 or 2000 data sets, respectively, 
%     to make a better model 
%
% polr from library(MASS) is also an option but it's not clear the categories
% in presvote_intent are plausibly ordered...



\end{enumerate}

% \end{enumerate}

\end{document}


