\documentclass{article}

\usepackage{bj,txfonts,hyperref}

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

\begin{document}

\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
HW08 -- Due Wed Nov 16, 11:59pm 
\end{center}

\begin{itemize}

\item Please turn the homework in, as a single pdf, online in GradeScope using the link provided on the Assignment page on canvas.cmu.edu.  

\item This hw heavily depends on material in lectures 21 and 22, which I will give ``just in time'' on Mon Nov 14 and Wed Nov 16.  
\begin{itemize}
\item In order that you can start earlier than that (after all, it's due Wed Nov 16!), I have already put the lecture notes for lectures 21 and 22 in the week11 folder in the files area on canvas.  
\item Also look closely at the R files, ``\verb|21 - mlm residuals.r|'' and ``\verb|22 - estimation and model selection.r|''
\end{itemize}

\item Weekly Reading \& Quiz:
\begin{itemize}
\item There is no new reading this week.  Everything is in the lecture notes.
\item There is no quiz on Monday. 
\end{itemize}

\item Information about the final project will be available sometime on Fri Nov 11.

\item  There are four main exercises below...

\end{itemize}


\section*{Exercises}

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

\item The file allvars.csv contains CD4 percentages (CD4PCT) for a set of 254 young children with HIV who were measured several times over a period of two years (A CD4 count measures the number of CD4 cells in your blood. It's used to check the immune system function in people with HIV\footnote{https://medlineplus.gov/lab-tests/cd4-lymphocyte-count/}). The dataset also includes the ages of the children at each measurement. This is an example of {\em growth curve data} in which the Level 1 observations are CD4PCT's at each visit, and the Level 2 ``group'' is a child.  How the CD4PCT's change over time for each child is the primary question of interest.  Because of skewing, you should replace CD4PCT with sqrt(CD4PCT) for all of the following questions.

\begin{enumerate}

\item For the first 12 children only ({\tt newpid $\leq$ 12}), make a facet plot of CD4PCT vs VISIT number. Use different plotting point colors for children who have treatmnt=1 and treatmnt=0\footnote{Presumably, ``treatmnt'' stands for some sort of treatment to improve childrens' HIV status.}.  Do the same for CD4PCT vs. visage\footnote{If you are using {\tt ggplot()}, you may wish to set {\tt scales="free\_x"} in the {\tt facet\_wrap()} function.} (the age at the child on each visit).  Which is a better measure of time?

\item Build a multilevel model with random slopes and intercepts for all the children, using a measure of time (VISIT or visage, whichever you answered in part (a)) as a Level 1 predictor, and treatmnt as a Level 2 predictor.  Report the estimated fixed effects and variances from the model, and add the fitted regression lines to your plot from part (a).

\item Try to expand the model by adding an interaction between treatmnt and your measure of time\footnote{Beacuse treatmnt is a level-2 variable and the time measure is a level-1 variable, this is somethings called a {\em cross-level interaction}.}.  Comment on changes in the estimated fixed effects and variances for this model, vs. the model in part (b).  Add the fitted regression lines from this model to your previous facets plot.  Does treatmnt account for an interesting amount of the variation in slopes (or intercepts) among children in the study? 

\end{enumerate}


\item Now we will consider in more detail the model below for the CD4 data set:
\begin{verbatim}
lmer.1 <- lmer(sqrt.CD4PCT ~ 1 + VISIT * treatmnt + (1 + VISIT | newpid),
               data=cd4)
\end{verbatim}
Bring together all the residuals (except random effects residuals) that we have talked about, using library(HLMdiag):
\begin{verbatim}
> r.1 <- hlm_resid(lmer.1,level=1,include.ls=F)
> r.1s <- hlm_resid(lmer.1,level=1,include.ls=F,standardize=T)
> r.2 <- hlm_resid(lmer.1,level="newpid",include.ls=F)
> r.2s <- hlm_resid(lmer.1,level="newpid",include.ls=F,standardize=T)
> names(r.1); names(r.1s); names(r.2); names(r.2s)
\end{verbatim}
Read the help file for \texttt{hlm\_resid} if necessary to understand what all the compenents are.  (You can create all these residuals ``by hand'' as I have shown in lecture slides and in the accompanying R files, but using HLMdiag is simpler and more reliable.)

\begin{enumerate}

\item Make a facets plot of the marginal residuals, as a function of the marginal fitted values (use \texttt{scales="free\_x"} if needed to make the plot legible).  Explain in a sentence or two why a facets plot is not very useful for assessing model fit for this problem, whether we look at the first 12 children, or all 251 children).

\item Make an ungrouped (that is, no facets) scatter plot of marginal residuals as a function of marginal fitted values, using the full data set (not just the first 12 children).  Color the points for \texttt{treatmnt=1} kids and treatmnt=2 kids with different colors.  Overlay a smooth fit (\texttt{geom\_smooth} is the easist to use here). Explain, in a couple of sentences (optionally with some math):
\begin{itemize}
\item What is causing the dominant structure in this plot, and why that dominant structure is essentially irrelevant for checking the relationship between \texttt{sqrt.CD4PCT} and \texttt{VISIT}; and  
\item What in this plot makes you happy or unhappy about having a linear relationship between  \texttt{sqrt.CD4PCT} and \texttt{VISIT} in the model.
\end{itemize}

\item Make an ungrouped (that is, no facets) scatter plot of conditional residuals as a function of conditional fitted values, using the full data set (not just the first 12 children).  Color the points for \texttt{treatmnt=1} kids and \texttt{treatmnt=2} kids with different colors.  Overlay a smooth fit (\texttt{geom\_smooth} is the easist to use here). Explain, in a couple of sentences (optionally with some math):
\begin{itemize}
\item Why the dominant structure in the marginal residuals is not also present in this plot of conditional residuals
\item What might be causing the trend you see in this plot to be different from the trend in the plot of the marginal residuals.
\end{itemize}

\item Use standardized residuals and standardized random effects estimates to assess the normality of $\epsilon_i$, $\eta_{0j}$ and $\eta_{1j}$ in the fitted model, and to check for any outliers.  Include qq plots for each, and accompany each plot with a sentence or two describing what is good or bad in that plot.

\end{enumerate}

\item Continuing with the CD4 data\ldots

\begin{enumerate}

\item Make a table giving values of AIC, BIC, DIC, and cAIC (you compared two of these on the last assignment,using just AIC, BIC and DIC):
\begin{verbatim}
sqrt.CD4PCT ~ 1 + visage +            (1+visage|newpid)
sqrt.CD4PCT ~ 1 + visage + treatmnt + (1+visage|newpid)
sqrt.CD4PCT ~ 1 + visage * treatmnt + (1+visage|newpid)
sqrt.CD4PCT ~ 1 + VISIT  +            (1+VISIT|newpid)
sqrt.CD4PCT ~ 1 + VISIT  + treatmnt + (1+VISIT|newpid)
sqrt.CD4PCT ~ 1 + VISIT  * treatmnt + (1+VISIT|newpid)
\end{verbatim}
Comment briefly on any similarities or differences in how the different criteria choose fixed effects.

\item Now let's look at the random effects in the model 
\begin{verbatim}
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid)
\end{verbatim}
Again, make a table giving values of  AIC, BIC, DIC, and cAIC for the following models:
\begin{verbatim}
sqrt.CD4PCT ~ 1 + VISIT + treatmnt 
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1|newpid)
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (0+VISIT|newpid)
sqrt.CD4PCT ~ 1 + VISIT + treatmnt + (1+VISIT|newpid)
\end{verbatim}
You'll fit the first model with lm(), and the others with lmer().  Comment briefly on any similarities or differences in how the different criteria choose random effects.  (Note that only BIC and cAIC have a strong theoretical justification here).

\item Repeat part (b) but with the interaction \verb|VISIT * treatmnt| in each model instead of just the main effects \verb|VISIT + treatmnt|.

\end{enumerate}


\item IMRAD and IDMRAD.  
\begin{enumerate}

\item {\em (You do not need to turn anything in for this part.)}  Compare the papers
\begin{itemize}

\item \verb|COVID breakthrough rates in England - PIIS1473309921004606.pdf| and\\ \verb|COVID breakthrough... (appendix) 1-s2.0-S1473309921004606-mmc1.pdf|

\item \verb|menu pricing IDMRAD paper with technical appendix.pdf|
\end{itemize}
The COVID paper is a published IMRAD paper with a separate technical appendix.  The ``menu pricing'' paper is a simple IDMRAD paper (with technical appendix attached) that you can use as a model for your final project paper.  Note similarities and differences.  In particular, when you write your paper, the technical appendix should show all the R code and results that you used to answer the questions posed for the data analysis\footnote{A small portion of the results in the appendix actually makes it into the main paper.}.  This is a little different from the published COVID paper, which really just provides some supplementary details in the technical appendix that didn't fit in the main paper.

\item {\em (You do need to turn something in for this part.)} Refer again to the two papers in part (a).  The difference between an IMRAD paper and an IDMRAD paper is that the ``Data'' part of the \textbf{Methods} section of an IMRAD paper is taked out of \textbf{Methods} and made into its own \textbf{Data} section between \textbf{Introduction} and \textbf{Methods}. \\[-0.75em]

Find the \textbf{Methods} section of the COVID paper, several pages in, and within the \textbf{Methods} section find the subsections that describe the data (rather than describing methods of analysis).  Copy and paste these sections (including subsection titles) into a separate \textbf{Data} section, and turn that section in as an answer to this question.


\end{enumerate}


\end{enumerate}


\end{document}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% great problem below!


\item The file \texttt{olympics1932.txt} has the judging data from the 1932 Olympics pairs figure skating competition (Yes, there is more current data, but this was easy to obtain!).  You can read the data from the file into R, but you will have to skip the first few lines, and a complete description of the scoring system and the awarding of medals is given at \url{https://www.olympedia.org/results/12760}, but in this problem we will use MLM's to separate skating quality from judging harshness, and to explore how close the skaters' performances were to one another.

After you read the data in, you will note that each judge gives two scores to each pair: a Performance and a Program score.  Arrange the data into a ``tall'' data frame with the following variables:
\begin{description}
\item[Pair:] Factor with labels 1 to 7, indicating which Skating Pair is being scored
\item[Judge:] Fsctor with labels 1 to 7, indicating which Judge is giving a score
\item[Category:] Either Performance or Program, indicating in which category the score is being given
\item[Score:] The score itself
\end{description}
Make sure that your new data frame is consistent with the data frame you read in from \texttt{olympics1932.txt}. 

In this problem we have two sets of individuals contributing to scores: Pairs, and Judges. Social scientists have developed some ``philosophical'' arguments for when individuals should be modeled as a fixed effect, and when they should be modeled as a random effect.  Briefly,
\begin{itemize}
\item If the individuals can be thought of as a complete set of all individuals that could be in that situation, then model the individuals as fixed effects.
\item If the individuals can be thought of as a sample from a larger population of individuals, then model the individuals as random effects.
\end{itemize}
This leads to four possible models for this data, depending on which individuals we think of as a complete set (fixed effects) or a random sample (random effects):
\begin{verbatim}
M1 <- lm(Score ~ Pair + Judge)
M2 <- lmer(Score ~ Pair + (1|Judge))
M3 <- lmer(Score ~ Judge + (1|Pair))
M4 <- lmer(Score ~ 1 + (1|Judge) + (1|Pair))
\end{verbatim}
There is no clear way to decide among these models on the basis of those philosophical considerations, but from a more practical point of view, models \texttt{M3} and \texttt{M4} have the advantage that they directly provide seven scores for the Pairs, adjusted for any Judge effects, in the form of the Pair random effects.

\begin{enumerate}

\item Fit models M3 and M4, produce display()'s or summary()'s of each, and comment on similarities or differences between them.

\item How would you rank skating pairs based on the Pair random effects from M3? How about based on the Pair random effects from M4?  How does this compare with the rankings you would get from the total scores (Points column) from \url{https://www.olympedia.org/results/12760}? \textit{(Medals were not determined from total scores in the 1932 Olympics, but rather from a majority-ranking scheme based on each judge's individual totals.)}

\item Now we will explore the uncertainty in the scores, using model M4 (we would get similar results from M3).  Use the \texttt{sims()} function from \texttt{library(arm)}, as shown in lecture 08, to simulate 1000 new sets of parameters from model M4.  Provide a table with the following, for each Pair:
\begin{itemize}
\item Mean, and $\pm$2 SD interval, from the 1000 simulated random effects for each Pair
\item Median, and 0.025$^{th}$ and 0.975$^{th}$ quantiles, from the 1000 simulated random effects from each Pair
\end{itemize}
Each of the above provides a point estimate $\hat\eta$ and an approximate 95\% interval for $\eta_j$, for each skating pair, $j=1,\ldots,7$.  If the simulated random effects are approximately normal, the two sets of intervals will agree.  If the distribution of the simulated effects is not approximately normal, the first set is invalid, but the second set can still be used.
\item Make a plot of the point estimates and 95\% intervals based on medians and quantiles.  Comment on which pairs of skaters are distnguishable in the sense that their interavals do not overlap.  is it possible, on the basis of non-overlapping intervals, to distinguish between the three medal winners, or between them and the other skating pairs?
\end{enumerate}





1932 Olympics data

\item Variable selection for london schools data
\begin{itemize}
\item Consider all existing covartiates
\item Consider transformation of LRT (e.g. poly(3,LRT)) in both fixed and random effects
\item Consider both AIC and cAIC
\end{itemize}


\end{enumerate}

\end{document}
