\documentclass{article}

\usepackage{bj,txfonts,hyperref}

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

\begin{document}

\begin{center}
\Large
36-463/663: Multilevel \& Hierarchical Models\\
Spring 2022\\
HW06 -- Due Tue Mar 29, 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.  Upload \ul{\em one} file per person.  


\item There are 2 main exercises, each with  ``parts''\ldots

\item These are not ``interesting'' exercises by themselves.  Instead they are more like ``finger exercises'' to practice writing, debugging, and examining the output of, \verb|stan()| programs.  
\begin{itemize}
\item Feel free to borrow ideas liberally from lecture notes, pdfs and R files shared in class, and/or from resources you find on the www, etc.
\end{itemize}

\end{itemize}


\section*{Exercises}

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

\item {\bf Chicks.}  The \verb|ChickWeight| data frame is included with R, and you can access it with \verb|data(ChickWeight)|.  It has 578 rows and 4 columns from an experiment on the effect of diet on early growth of chicks.  The variables (Columns) are:
\begin{itemize}
\item weight: a numeric vector giving the body weight of the chick (gm).

\item Time: a numeric vector giving the number of days since birth when the measurement was made.

\item Chick: an ordered factor with levels 18 $<$ $\ldots$ $<$ 48 giving a unique identifier for the chick. The ordering of the levels groups chicks on the same diet together and orders them according to their final weight (lightest to heaviest) within diet.

\item Diet: a factor with levels 1, ..., 4 indicating which experimental diet the chick received.
\end{itemize}
For most of this exercise we will be working with the model
\begin{eqnarray*}
  \log(\mbox{weight}_i) & = & \alpha_{0j[i]} + \alpha_{1j[i]}\cdot \mbox{Time}_i + \epsilon_i \\
  \alpha_{0j} & = & \beta_0 + \eta_{0j} \\
  \alpha_{1j} & = & \beta_1 + \eta_{1j} \\
  \epsilon_i & \stackrel{iid}{\sim} & N(0,\sigma^2) \\
  \eta_{0j}  & \stackrel{iid}{\sim} & N(0,\tau_0^2) \\
  \eta_{1j}  & \stackrel{iid}{\sim} & N(0,\tau_1^2)~,   
\end{eqnarray*}
where $i$ is the observation number and $j$ is the Chick number\footnote{Also note that in this model we are just lettng $\rho_{\eta_0,\eta_1} = 0$, to make the {\tt stan()} coding less complicated.}.  You will have to convert \verb|Chick| from a factor variable to a numerical variable, to work with \verb|stan()|.  

\begin{enumerate}

\item In the same folder as this assignment sheet there is a file \verb|chicks.stan| {\em with most of the program already written, but the model section left blank.}  Fill in the blank model section\footnote{You should not have to change any of the code already in the file; just add to it.  However if it is easier for you to change some of the code in the file, that is fine too.} and get the program to run using \verb|library(rstan)|, taking at least 500 MCMC steps per chain\footnote{If you need to take more steps to get convergence to the stationary distribution, that is fine.}, and accepting the other \verb|stan()| defaults, so that you end up with at least 1000 ``good'' MCMC samples from the four Markov chains produced.  Turn in:
\begin{itemize}
\item The completed stan program \verb|chicks.stan|.
\item A printout of the fitted stan object but including only the lines for $\beta_0$, $\beta_1$, $\tau_0$, $\tau_1$, and $\sigma$.
\end{itemize}
(You are free to use \verb|lmer()| to estimate the corresponding multi-level model, to make sure you are on the right track, but please don't turn in any \verb|lmer()| code or output for this part of the exercise.)

\item Use \verb|shinystan| or \verb|bayesplot| to generate and turn in the following plots, neatly organized:
\begin{itemize}
\item For $\beta_0$: (i) Histogram of the MCMC draws; (ii) Autocorrelation plot for at least one of the four Markov chains; (iii) trace plot of all four Markov chains on the same graph, with different colors to indicate the different Markov chains.
\item For $\beta_1$: the same three things.
\item For $\tau_0$: the same three things.
\item For $\tau_1$: the same three things.
\item For $\sigma$: the same three things.
\end{itemize}

\item Set \verb|y <- log(ChickWeight$weight)|, and use \verb|shinystan| or \verb|bayesplot| to produce graphical posterior predictive checks (ppc's) for 
\begin{itemize} 
\item \verb|T(y) = mean(y)|
\item \verb|T(y) = sd(y)|
\item At least one other test statistic \verb|T(y)| that produces an interesting result.
\end{itemize}

\item Now fit the corresponding \verb|lmer()| model (the parameter estimates should be almost identical to your estimates from the \verb|stan()| model). Turn in 
\begin{itemize}
\item The plot of the conditional residuals for the \verb|lmer()| model.  If the model fits well, these should be approximately normally distributied.
\item The plot of the ``uniformized'' conditional residuals\footnote{We could build these by hand from the {\tt stan()} model, but {\tt DHARMa} is faster and less prone to user errors.} for the \verb|lmer()| model, using \verb|library(DHARMa)|.  If the model fits, these should be approximately uniformly distributed.
\end{itemize}
Write a sentence or two identifying and commenting on any suggested improvements to the model, or any other useful information, that you can see in these plots.  


\end{enumerate}


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


\item{\bf Toenails.} The \verb|toenail| data frame comes from the \verb|faraway| library, and you can access it by 
\begin{verbatim}
install.packages("faraway")
library(faraway)
data(toenail)
\end{verbatim}
There are 1908 observations from a study comparing two oral treatments for toenail infection. Patients were evaluated for the degree of separation of the nail. Patients were randomized into two treatments and were followed over seven visits---four in the first year and yearly thereafter. The patients have not been treated prior to the first visit so this should be regarded as the baseline.  The variables (columns) are:
\begin{itemize}
\item ID:
ID of patient

\item outcome:
0=none or mild seperation, 1=moderate or severe

\item treatment:
the treatment A=0 or B=1

\item month:
time of the visit (not exactly monthly intervals hence not round numbers)

\item visit:
the number of the visit
\end{itemize}
For most of this exercise we will be working with the model
\begin{eqnarray*}
  \mbox{outcome}_i & \sim & \mbox{Bernoulli}(p_i) \\
  \mbox{logit}(p_i) & = & \alpha_{0j[i]} + \alpha_{1j[i]}\cdot \mbox{month}_i + 
  \beta_2\cdot\mbox{treatment}_{j[i]} + \epsilon_i \\
  \alpha_{0j} & = & \beta_0 + \eta_{0j} \\
  \alpha_{1j} & = & \beta_1 + \eta_{1j} \\
  \eta_{0j}  & \stackrel{iid}{\sim} & N(0,\tau_0^2) \\
  \eta_{1j}  & \stackrel{iid}{\sim} & N(0,\tau_1^2)   
\end{eqnarray*}
where $i$ is the observation number and $j$ is the patient ID number\footnote{Also note that in this model we are just lettng $\rho_{\eta_0,\eta_1} = 0$, to make the {\tt stan()} coding less complicated.}.  Note that there are some skips in the ID number sequence in the data set; it would be simplest for working with \verb|stan()| to convert these ID numbers to successive integers with no skips.
(Note that $\sigma$ is missing from this model---why?)

\begin{enumerate}
\item In the same folder as this assignment sheet there is a file \verb|toenails.stan| {\em with the parameters and transformed parameters sections left blank}.  Fill in the blank sections\footnote{You should not have to change any of the code already in the file; just add to it.  However if it is easier for you to change some of the code in the file, that is fine too.} and get the program to run using \verb|library(rstan)|, taking at least 500 MCMC steps per chain\footnote{If you need to take more steps to get convergence to the stationary distribution, that is fine.}, and accepting the other \verb|stan()| defaults, so that you end up with at least 1000 ``good'' MCMC samples from the four Markov chains produced.  Turn in:
\begin{itemize}
\item The completed stan program \verb|toenails.stan|.
\item A printout of the fitted stan object but including only the lines for $\beta_0$, $\beta_1$, $\beta_2$, $\tau_0$, and $\tau_1$.
\end{itemize}
Some notes:
\begin{itemize}
\item You are free to use \verb|glmer()| to estimate the corresponding multi-level model, to make sure you are on the right track, but please don't turn in any \verb|glmer()| code or output for this part of the exercise. Note that, even if you have done everything right the magnitudes of the \verb|glmer()| estimates may not agree very well with the magnitudes of the \verb|stan()| estimates; see part (d) below also.
\item You will struggle to get a ``great'' run from \verb|stan()| on this model.  If you can get a run where most or all of the difficulties are in \verb|yrep|, that will be good enough: the various convergence and stability diagnostics like \verb|Rhat| don't really apply to \verb|yrep| since it is simulated from the posterior predictive distribution; moreover the discreteness of \verb|yrep| makes $n_{eff}$ less meaningful.
\end{itemize}

\item Use \verb|shinystan| or \verb|bayesplot| to generate and turn in the following plots, neatly organized:
\begin{itemize}
\item For $\beta_0$: (i) Histogram of the MCMC draws; (ii) Autocorrelation plot for at least one of the four Markov chains; (iii) trace plot of all four Markov chains on the same graph, with different colors to indicate the different Markov chains.
\item For $\beta_1$: the same three things.
\item For $\beta_2$: the same three things.
\item For $\tau_0$: the same three things.
\item For $\tau_1$: the same three things.
\end{itemize}

\item Set \verb|y <- toenail$outcome|, and use \verb|shinystan| or \verb|bayesplot| to produce graphical posterior predictive checks (ppc's) for 
\begin{itemize} 
\item \verb|T(y) = mean(y)|
\item \verb|T(y) = sd(y)|
\item At least one other test statistic \verb|T(y)| that produces an interesting result.
\end{itemize}

\item Now fit the corresponding \verb|glmer()| model (the parameter estimates should be almost identical to your estimates from the \verb|stan()| model).  When I did this, I thought that fixed effect estimates from \verb|glmer()| were much less plausible than the fixed effect estimates from \verb|stan()|, which made me trust the \verb|stan()| results more.  Do you agree?  Why or why not?


\item Turn in 
\begin{itemize}
\item The plot of the conditional residuals for the \verb|glmer()| model.  If the model fits well, these should be approximately normally distributied.
\item The plot of the ``uniformized'' conditional residuals\footnote{We could build these by hand from the {\tt stan()} model, and we probably should since the {\tt stan()} results seem more plausible than the {\tt glmer()} results, but {\tt DHARMa} is faster and less prone to user errors; it will be fine for the purpose of just getting used to using and interpreting {\tt DHARMa} stuff.} for the \verb|lmer()| model, using \verb|library(DHARMa)|.  If the model fits, these should be approximately uniformly distributed.
\end{itemize}
Write a sentence or two identifying and commenting on any suggested improvements to the model, or any other useful information, that you can see in these plots (or if you think some or all of the plots are useless, explain that as well).  


\end{enumerate}


\end{enumerate}


\end{document}