\documentclass{article}

\usepackage{bj,txfonts,hyperref}

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

\begin{document}

\begin{center}
\Large
36-617: Applied Linear Models\\
Fall 2022\\
HW10 -- Due Wed Dec 7, 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, as usual.  

\item Reading: for the week after spring break:
\begin{itemize}
\item For last week: Lynch, Ch's 3 (read) and 4 (skim)
\item For this week \& next week: Lynch, Ch 9 (read)
\end{itemize}
All three readings will be in the week 13, week14, and hw10 folders on Canvas.

\item  There are 3 exercises
\begin{itemize}
\item The first one is just to get you going with \verb|stan()|.  {\em Instructions for installing Stan on your laptop are on the last page of this hw.}

\item The next two are about applying \verb|stan()| to estimate Bayesian versions of multilevel models.  They  are not particularly ``interesting'' exercises by themselves.  Instead they are more like ``finger exercises'' to play around with writing, debugging, and examining the output of, \verb|stan()| programs.  

\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 Install and verify the software package Stan and the library rstan, on your personal computer. Instructions for doing this are appended at the end of this hw assignment.  Then run the following short example to verify that rstan is intalled and working correctly (note that the text string assigned to the variable ``model'' runs over several lines; that's perfectly OK):
\begin{verbatim}
text_model <- "
   data {
     real y_mean;
   } 
   parameters {
     real y;
   } 
   model {
     y ~ normal(y_mean,1);
   }"
compiled_model <- stan(model_code=text_model,data=list(y_mean=0),chains=0)
fit <- stan(fit=compiled_model,data=list(y_mean=0))
fit2 <- stan(fit=compiled_model,data=list(y_mean=5))
\end{verbatim}
This Stan program is not doing anything Bayesian, it's just making draws from a normal distribution with mean=0 and sd=1 (fit) or mean 5 and sd=1 (fit2).

\begin{enumerate}

\item Turn in the output for both \verb|print(fit)| (or just \verb|fit|), and \verb|summary(fit)|.  Write a sentence or two explaining the difference between the two outputs (besides a different number of decimal places printed!).

\item Use commands like those in the appendix (from \verb|library(bayesplot)|) of lecture 13 to produce and turn in, for the fitted model \verb|fit|:
\begin{itemize}
\item Trace plots of \verb|y| (the variable being simulated) and \verb|lp__| (the value of the log-posterior (in this case, not a posterior but just the normal density) at each simulated value of \verb|y|). 
\item Autocorrelation plots of \verb|y| and \verb|lp__|.
%  \item Plots of \verb|rhat| for each of \verb|y| and \verb|lp__|.
%  \item Plots of \verb|neff| ratios for each of \verb|y| and \verb|lp__|.
\item A histogram of just \verb|y|. 
\end{itemize}
\item Use the \verb|extract| function to extract simulated values for \verb|y| from \verb|fit| and make a qqplot of the values against the Normal distribution.  Is it plausible that the simulated values for \verb|y| are normally distributed?
\item Print out and turn in output from \verb|print(fit)| and \verb|print(fit2)|.  Write a sentence noting any similarities or differences between them.
\item Use the \verb|extract| function to extract simulated values for \verb|y| from \verb|fit2| and make a qqplot of the values against the Normal distribution.  Is it plausible that the simulated values for \verb|y| are normally distributed?

\end{enumerate}

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

\clearpage

\begin{center}
\Large\bf Installing Stan  for R
\end{center}

\begin{enumerate}

\item Go to \url{https://mc-stan.org/users/interfaces/rstan} \\ (can also find by
   googling ``rstan'').

\item Click on the ``RStan Quick Start Guide'' under ``Download and Get
   Started''

   ({\em Note:} We will be installing Stan to run in R.  There are also
   installations for python, matlab, julia, stata, mathematica, scala,
   etc\ldots{}  If you are interested, see \url{https://mc-stan.org/} for
   details.)

\item Follow the instructions for configuring the \verb|C++| Toolchain on your
   computer
   
   \begin{itemize}

   \item This is the trickiest part of installing Stan in R on your
     computer.

   \item There are separate instructions for Windows, Mac, \& Linux

   \item Make sure you are using R version 4.0 or higher.  You may want to
     update Rstudio to the current version as well (you need at least
     version 1.2.5042).

   \item {\bf MAKE SURE YOU ARE} following the instructions for your version of
     R (version 4.0 or higher) not an earlier version of R such as
     version 3.6, etc.
     
   \end{itemize}
     
\item Follow the instructions under ``Installation of RStan'' and ``Verify
   Installation''
   
   \begin{itemize}

   \item The most important parts are the ``install.packages'' command, to
     actually do the installation, and the ``example'' command to verify
     that your installation is good.  If you can get these two
     commands to run successfully, you are golden.
     
     \begin{itemize}

     \item Note: the ``example'' command will take from several seconds to
        up to a minute or so to run.  It will also generate a ton of
        output that may not make much sense right now\ldots{} Be patient\ldots{}

     \end{itemize}

   \item Please note the other instructions under these two headings,
      which handle some common special cases that you may need to deal
      with on your laptop.

   \end{itemize}

\item After you have verified that your installation is good, feel free
   to browse through the longer ``How to Use RStan'' section.  As you
   get more accustomed to Stan, you'll see that the Stan modeling
   language is a kind of extension of R's usual programming language
   -- much more general and flexibe than \verb|lmer()|'s modeling language.

\item Check out \url{https://mc-stan.org/} for background info on stan.

\end{enumerate}


\end{document}
