---
title: "36-617 - Fall 2022 -- HW10 Solutions"
author: "BJ"
date: '2022-12-02'
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



```{r warning=F, message=F}
set.seed(1234567890) ## to make a reproducible simulation result

library(arm)
library(ggplot2); theme_set(theme_bw())
library(gridExtra)

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

library(bayesplot)
library(shinystan)

library(DHARMa)

print_file <- function(filename) {
  x <- paste(scan(file=filename,
                  what="",
                  sep="\n",
                  blank.lines.skip = TRUE,
                  quiet=TRUE),
             collapse="\n")
  cat(x)
}

```

# Problem 1.

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):
```{r warning=F, message=F}

set.seed(1234567890) ## to make a reproducible simulation result 

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(ggplot2); theme_set(theme_bw())
library(bayesplot)
```


```{r}
text_model <- "
   data {
     real y_mean;
   } 
   parameters {
     real y;
   } 
   model {
     y ~ normal(y_mean,1);
   }"

```
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).

## Problem 1(a).

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!).



```{r}
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))
```


```{r}
## fit or print(fit)
fit
## summary(fit)
summary(fit)
```

\it
Here's a comparison of the two outputs:
\begin{itemize}

\item \verb|print(fit)| just gives a summary of the estimated parameters, \verb|se_mean| from the Monte Carlo sumualtion, the posterior \verb|sd|, quantiles or the posterior distribution, and \verb|n_eff| \& \verb|Rhat| values, for the "kept" part of all the chains, merged together.

\item \verb|summary(fit)| gives the same overall summary, as a two-dimensional array, in the \verb|$summary| component, and then in \verb|$c_summary| it provides a 3-dimensional array containing summaries for each individual Markov chain that was run (but without the \verb|n_eff| \& \verb|Rhat| values). 

\item[] (The two-dimensional array in \verb|summary(fit)$summary| is convenient for extracting mean, sd, quantiles, etc. to make CI's. etc.)

\end{itemize}

\rm

## Problem 1(b).

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 A histogram of just \verb|y|. 
\end{itemize}

```{r fig.height=5}
library(gridExtra)
g1 <- mcmc_trace(fit)
g2 <- mcmc_acf(fit)
g3 <- mcmc_areas(fit,pars="y")
grid.arrange(arrangeGrob(g1,g2,ncol=2),g3,heights=c(2.5/4,1.4/4),ncol=1)
```


## Problem 1(c).

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?

```{r fig.height=2.75}
sim <- extract(fit)

ggplot(as.data.frame(sim),aes(sample=y)) +
  geom_qq() +
  geom_qq_line() +
  xlab("Normal Quantiles") + 
  ylab("y (from \"fit\")")
```
\it

Yes, given this plot, it's pretty plausibly that the y's are normally distributed. There's just a hint of a deviation from the normal in the tails (where there's very little sample size, so it's hard to say these are "significant" deviations), but overall the data follow the normal distribution well.

\rm

## Problem 1(d).

Print out and turn in output from \verb|print(fit)| and \verb|print(fit2)|.  Write a sentence noting any similarities or differences between them.


```{r}
fit
fit2
```

\it

Everything is pretty similar for the two outputs, except of course the estimated mean (as either a posterior mean or a posterior median) is $\approx$ 0 for \verb|fit| and it is $\approx$ 5 for \verb|fit2|.  The posterior quantiles are similarly shifted by about 5 units.

\rm

## Problem 1(e).

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?
```{r fig.height=2.75}
sim <- extract(fit2)

ggplot(as.data.frame(sim),aes(sample=y)) +
  geom_qq() +
  geom_qq_line() +
  xlab("Normal Quantiles") + 
  ylab("y (from \"fit2\")")
```

\it

Again this qq plot is pretty convincing of normality.  There's just a hint of a deviation from the normal in the tails (where there's very little sample size, so it's hard to say these are "significant" deviations).  Note that the y-intercept for the fitted line is now about 5, corresponding to the fact that we are simulating from a $N(5,1)$ instead of a $N(0,1)$ this time.

\rm



# Problem 2 -- 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 letting $\rho_{\eta_0,\eta_1} = 0$, to make the \texttt{stan()} coding less complicated.}.  You will have to convert \verb|Chick| from a factor variable to a numerical variable, to work with \verb|stan()|.

## Problem 2(a)

In the same folder as this assignment sheet there is a file \verb|chicks.stan| \textit{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.)

```{r}
data(ChickWeight)

## str(ChickWeight)
```
\it
Since the Level 1 mean \texttt{mu[i] <- a0[Chick[i]] + a1[Chick[i]]*Time[i]} and the random intercept \texttt{a0[j] <- b0 + eta0[j]} and slope \texttt{a1[j] <- b1 + eta1[j]} are already defined in the \texttt{transformed parameters} block, all we need to do in the \texttt{model} block is specify
\begin{itemize}

\item Level 1: distribution (likelihood) for $y$:
\begin{verbatim}
  for (i in 1:N) {
    y[i] ~ normal(mu[i],sigma);
  }
\end{verbatim}

\item Level 2: distributions of the $\eta$'s:
\begin{verbatim}
  for (j in 1:J) {
    eta0[j] ~ normal(0,tau0);
    eta1[j] ~ normal(0,tau1);
  }
\end{verbatim}

\item Priors for all the free parameters (fixed effects and sd's):
\begin{verbatim}
  b0 ~ normal(0,1e6);
  b1 ~ normal(0,1e6);
  sigma ~ uniform(0,50);
  tau0 ~ uniform(0,50);
  tau1 ~ uniform(0,50);
\end{verbatim}

\end{itemize}

The final, complete stan program looks like this:
\rm
```{r}
print_file("chicks.stan")

```
\it
Examining the data block, 
\begin{verbatim}
data {
  int N;  // # observations
  int J;  // # chicks
  real weight[N];
  real Time[N];
  int Chick[N]; // needs to be integers, not a factor variable
}
\end{verbatim}
all of the variables are already in the ChickWeight data frame, except for \verb|N|, \verb|J| and a version of \verb|Chick| that is an integer, not a factor, variable.  Here is code to run the model:
\rm
```{r echo=T, eval=T}
attach(ChickWeight)

chicknum <- as.numeric(Chick)

chickdata <- list(N=length(chicknum),
                  J=length(unique(chicknum)),
                  weight=weight,
                  Time=Time,
                  Chick=chicknum)

detach()

chickmodel <- stan(file="chicks.stan",data=chickdata,chains=0)

chickresult <- stan(fit=chickmodel,data=chickdata,iter=500)
```
\it
and here is the printout of the results, for $\beta_0$, $\beta_1$, $\tau_0$, $\tau_1$ and $\sigma$:
\rm
```{r echo=F}
## load(".RData")  
## load a copy of the R variables that already has the stan results saved, so 
## that I don't have to re-run stan each time I re-knit the rmarkdown file.
```

```{r}
print(chickresult,pars=c("b0","b1","tau0","tau1","sigma"))
```


## Problem 2(b)

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}
```{r fig.height=10, fig.width=7.5, warning=F, message=F}
param <- "b0"
g1 <- mcmc_hist(chickresult,pars=param)
g2 <- mcmc_acf(chickresult,pars=param)
g3 <- mcmc_trace(chickresult,pars=param)
param <- "b1"
g4 <- mcmc_hist(chickresult,pars=param)
g5 <- mcmc_acf(chickresult,pars=param)
g6 <- mcmc_trace(chickresult,pars=param)
param <- "tau0"
g7 <- mcmc_hist(chickresult,pars=param)
g8 <- mcmc_acf(chickresult,pars=param)
g9 <- mcmc_trace(chickresult,pars=param)
param <- "tau1"
g10 <- mcmc_hist(chickresult,pars=param)
g11 <- mcmc_acf(chickresult,pars=param)
g12 <- mcmc_trace(chickresult,pars=param)
param <- "sigma"
g13 <- mcmc_hist(chickresult,pars=param)
g14 <- mcmc_acf(chickresult,pars=param)
g15 <- mcmc_trace(chickresult,pars=param)
grid.arrange(g1,g2,g3,
             g4,g5,g6, 
             g7,g8,g9,
             g10,g11,g12,
             g13,g14,g15,ncol=3)
```


## Problem 2(c)

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}

\it
In addition to the mean and the sd, I tried the min and the max as test statistics (see below).  The mean, sd and min didn't show any misfit of the model to the data, but you can see that the max of the data is far outside the distribution of max's of data simulated from the model (from the posterior predictive distribution).  That is interesting, because it suggests a feature of the data that the model isn't fitting very well.  

Since the max of the real data is lower than the max of the simulated data, it suggests to me that perhaps $\log(weight)$ was too strong a transformation, and maybe something like $\sqrt{weight}$ or some other fractional power would be better.  (The fit would have to be much much better for the fractional power, though, since "log" is so much easier to explain to a client or colleague than a fractional power of weight would be.)
\rm

```{r fig.height=4, warning=F, message=F}
y <- log(ChickWeight$weight)
yrep <- extract(chickresult,pars="yrep")$yrep

g1 <- ppc_stat(y,yrep,stat="mean")
g2 <- ppc_stat(y,yrep,stat="sd")
g3 <- ppc_stat(y,yrep,stat="min")
g4 <- ppc_stat(y,yrep,stat="max")

grid.arrange(g1,g2,g3,g4,ncol=2)

```


## Problem 2(d)

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 distributed.
\item The plot of the "uniformized" conditional residuals\footnote{We could build these by hand from the \texttt{stan()} model, but \texttt{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.

```{r}
lmer.1 <- lmer(log(weight) ~ Time + (Time|Chick), data=ChickWeight, REML=F)
```
\it
First I'll verify that the parameter estimates from lmer() and stan() are about the same...
\rm
```{r}
round(summary(lmer.1)$coef,2)
VarCorr(lmer.1)

print(chickresult,pars=c("b0","b1","tau0","tau1","sigma"))
```
\it
\ldots and then here are the plots\ldots
\clearpage
\rm
```{r fig.height=3.75}
plot(lmer.1) ## "plot" plots conditional residuals for (g)lmer objects.

plot(simres <- simulateResiduals(lmer.1))
```
\it
The qqplot looks good, though the KS test for uniformity of the "uniformized" residuals is just at the edge of significance.  On the other hand, both residual vs fitted plots suggest that we missed some curvature in the model; I might try to fix this by adding \verb|I(Time^2)| or \verb|I(Time^2) + I(Time^3)| to the model.
\rm

## Problem 3 -- 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 separation, 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 letting $\rho_{\eta_0,\eta_1} = 0$, to make the \texttt{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?)

\it
(There is no $\sigma$ in the model since the likelihood is Bernoulli, not normal, and for the Bernoulli, once you know the mean $\mu_i = p_i$, you also know the SD $=$ $\sqrt{p_i(1-p_i)}$.)
\rm
```{r warning=F, message=F}
library(faraway)

data(toenail)
str(toenail)

## You can see where the skips are in the ID numbers by doing this :
##
## diff(sort(unique(toenail$ID)))
##
## (if there were no skips you would see all 1's in the output of this command)

## ...and a very simple way to renumber the ID numbers so there are no skips is this:
ID.new <- factor(toenail$ID)
ID.new <- as.numeric(ID.new)

## You can verify that things line up correctly with 
##
## View(cbind(toenail$ID,ID.new))
```


## Problem 3(a)

In the same folder as this assignment sheet there is a file \verb|toenails.stan| \textit{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}

\it
Clearly we have to define all of the parameters in the \texttt{model} block; the only tricky part is knowing what to put in the \texttt{parameters} block, vs. the \texttt{transformed parameters} block.  But we can imitate what was done in \texttt{chicks.stan}, and what we come up with is this:
\rm
```{r}
print_file("toenails.stan")
```
\it
Now, consulting the \texttt{data} block to see what goes in the data list for stan(), we can use the following code to make a stan run:
\rm
```{r echo=T, eval=T}
attach(toenail)

ID.new <- factor(ID)
ID.new <- as.numeric(ID.new)

toedata <- list(N=length(ID.new),
                J=length(unique(ID.new)),
                y=outcome,
                month=month,
                ID=ID.new,
                treatment=treatment)
detach()

toemodel <- stan(file="toenails.stan",data=toedata,chains=0)

toeresult <- stan(fit=toemodel,data=toedata,iter=500)
```

```{r echo=F}
## load(".RData")
## load a copy of the R variables that already has the stan results saved, so 
## that I don't have to re-run stan each time I re-knit the rmarkdown file.
```

\it 
\ldots and the requested estimated parameters are:
\rm
```{r}
print(toeresult,pars=c("b0","b1","b2","tau0","tau1"))
```
\it

The stan() result is not super (small \verb|n_eff|'s and some \verb|Rhat|'s bigger than 1.1).  I could run the stan model longer (say with iter=1000 or iter=2000 rather than iter=500) to get better \verb|Rhat|'s. but I'm going to go with what I already have from stan(); these \verb|Rhat|'s are OK for just poking around with data and models as we are doing here.

(I would definitely try a higher number of iterations, and generally do what is needed to avoid warnings and things, if this were a final run for a formal paper or report!)

\rm


## Problem 3(b)

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}
```{r fig.height=10, fig.width=7.5, warning=F, message=F}
param <- "b0"
g1 <- mcmc_hist(toeresult,pars=param)
g2 <- mcmc_acf(toeresult,pars=param)
g3 <- mcmc_trace(toeresult,pars=param)
param <- "b1"
g4 <- mcmc_hist(toeresult,pars=param)
g5 <- mcmc_acf(toeresult,pars=param)
g6 <- mcmc_trace(toeresult,pars=param)
param <- "b2"
g7 <- mcmc_hist(toeresult,pars=param)
g8 <- mcmc_acf(toeresult,pars=param)
g9 <- mcmc_trace(toeresult,pars=param)
param <- "tau0"
g10 <- mcmc_hist(toeresult,pars=param)
g11 <- mcmc_acf(toeresult,pars=param)
g12 <- mcmc_trace(toeresult,pars=param)
param <- "tau1"
g13 <- mcmc_hist(toeresult,pars=param)
g14 <- mcmc_acf(toeresult,pars=param)
g15 <- mcmc_trace(toeresult,pars=param)

## The graphs are on the next page.  You can see that the autocorrelations stay 
## pretty large for b1, tau0 and tau1, which is consistent with the small n_eff's 
## and large Rhat's we see for those parameters.

grid.arrange(g1,g2,g3,
             g4,g5,g6,
             g7,g8,g9,
             g10,g11,g12,
             g13,g14,g15,ncol=3)
```

## Problem 3(c)

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}

\it
Here are the plots for the mean and the SD.  No indications of misfit here.  I also tried the min and the max; nothing interesting there either.
\rm
```{r fig.height=2, warning=F, message=F}
y <- toenail$outcome
yrep <- extract(toeresult,pars="yrep")$yrep

g1 <- ppc_stat(y,yrep,stat="mean")
g2 <- ppc_stat(y,yrep,stat="sd")

grid.arrange(g1,g2,ncol=2)
```

\it

As for coming up with one other test statistic, you may have had some creative idea that I didn't think of, but here's my idea:

After some thought, I decided to check to see if there were too many or too few positive residuals in different ranges of the fitted values $\hat p_i$, where $p_i=P[outcome_i=1]$.  The idea is this
\begin{itemize}
\item Divide the range of fitted values into three groups
\begin{itemize}
\item The lower 1/3 of the range of fitted values $\hat p_i$
\item The middle 1/3 of the range of fitted values $\hat p_i$
\item The upper 1/3 of the range of fitted values $\hat p_i$
\end{itemize}
\item In each group:
\begin{itemize}
\item Calculate $T^{(obs)}$ $=$ number of positive residuals ($y_i-\hat p_i>0$) in the real data
\item For each replicated data set $y^{rep,m}$, $m=1,\ldots M$, calculate $T^{(rep,m)}$ $=$ number of positive residuals ($y^{rep,m}_i-\hat p^{rep,m}_i>0$) in the $m^{th}$ replicated data set.
\item Compare $T^{(obs)}$ to a histogram of $T^{(rep,m)}$'s, $m=1,\ldots M$.
\end{itemize}
\end{itemize}
If for example there are two many positive residuals in the lower and upper thirds, but too few in the middle third, that suggests a U-shape for the residuals, and perhaps I would add \verb|I(month^2)| to the model.


Here is the code to implement this idea, and the resulting plots:

\rm
```{r fig.height=5}

##### Get the data elements needed

y <- toenail$outcome  
## the real data

yrep <- extract(toeresult,pars="yrep")$yrep
## the replicated data sets; each row of this matrix is one replication of y

logit_p <- extract(toeresult,pars="logit_p")$logit_p
## the logits sampled from the posterior distribution again, each row of this 
## matrix contains the logits for one replication yrep

p <- invlogit(logit_p)
## the p's used to generate the yrep's each row of this matrix contains the p's
## used to generate Bernoulli draws for one replication yrep

p_hat <- apply(p,2,mean)
## posterior mean estimate of p[outcome=1]

##### now, make the plots...

par(mfrow=c(3,1))

T.num.pos <- function(x) sum(x>0)   ## the test statistic

##### Plot for the lower 1/3 of the range of p_hat

lo <- -1e6
hi <- min(p_hat) + (max(p_hat)-min(p_hat))/3
sel <- (p_hat<=hi) & (p_hat>lo)

T.obs <- T.num.pos((y-p_hat)[sel])
T.rep <-  apply((yrep-p)[,sel],1,T.num.pos)

hist(T.rep,xlim=c(min(c(T.obs,T.rep)),max(c(T.obs,T.rep))),
     main="T=Number of Positive Residuals\n(Lower 1/3 of range of p_hat)",
     xlab="T(yrep)")
abline(v=T.obs,col="red",lwd=2)
text(x=T.obs,y=25,labels="T(y)",pos=4)

##### Plot for the middle 1/3 of the range of p_hat

lo <- min(p_hat) + (max(p_hat)-min(p_hat))/3
hi <- min(p_hat) + 2*(max(p_hat)-min(p_hat))/3
sel <- (p_hat<=hi) & (p_hat>lo)

T.obs <- T.num.pos((y-p_hat)[sel])
T.rep <-  apply((yrep-p)[,sel],1,T.num.pos)

hist(T.rep,xlim=c(min(c(T.obs,T.rep)),max(c(T.obs,T.rep))),
     main="T=Number of Positive Residuals\n(Middle 1/3 of range of p_hat)",
     xlab="T(yrep)")
abline(v=T.obs,col="red",lwd=2)
text(x=T.obs,y=75,labels="T(y)",pos=4)

##### Plot for the upper 1/3 of the range of p_hat

lo <- min(p_hat) + 2*(max(p_hat)-min(p_hat))/3
hi <- +1e6
sel <- (p_hat<=hi) & (p_hat>lo)

T.obs <- T.num.pos((y-p_hat)[sel])
T.rep <-  apply((yrep-p)[,sel],1,T.num.pos)

hist(T.rep,xlim=c(min(c(T.obs,T.rep)),max(c(T.obs,T.rep))),
     main="T=Number of Positive Residuals\n(Upper 1/3 of range of p_hat)",
     xlab="T(yrep)")
abline(v=T.obs,col="red",lwd=2)
text(x=T.obs,y=50,labels="T(y)",pos=2)
```
\it

These plots show interesting misfit, but not the misfit I was expecting.  Instead, there are too few positive residuals in the lower range, about the right number in the middle, and too many positive residuals in the upper range.  That suggests that the residuals have an increasing trend, relative to the fitted model.  Since the residual is $y_i-\hat p_i$ and $y_i$ can only be 0 or 1, this says that the model under-predicts the number of cases of $y_i=0$ when $\hat p_i$ is low, and under-predicts the number of cases of $y_i=1$ when $\hat p_i$ is high: the data is both zero-inflated and one-inflated, relative to the model predictions.  There may be some other variable that we could collect and include in the model, that might help the model to correctly predict these 0's and 1'. 

\rm

## Problem 3(d)

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?

```{r}
glmer.1 <-  glmer(outcome ~ month + treatment + (month||ID),data=toenail,family=binomial)
## note use of "||" in the random effect notation to force the correlation between eta's
## to be zero, just like the stan() model...

display(glmer.1)
```

```{r echo=F}
b0 <- summary(glmer.1)$coef[1,1]
b1 <- summary(glmer.1)$coef[2,1]
b2 <- summary(glmer.1)$coef[3,1]
tau0 <- attr(VarCorr(glmer.1)$ID,"stddev")
tau1 <- attr(VarCorr(glmer.1)$ID.1,"stddev")
```

\it
The glmer() results just aren't very believable.  The intercept is estimated to be $\hat\beta_0$ $=$ `r round(b0,2)`; the corresponding probability is $1/(1+\exp(-\hat \beta_0))$ $=$`r invlogit(b0)` so that, without the random effects, the probability of a good outcome is essentially zero; $\hat \beta_1$ $=$ `r round(b1,2)` is just as large, and the treatment effect $\hat beta_2$ $=$ `r round(b2,2)` is minuscule and very different from the estimate from stan().  The SD's of the random effects are huge however ($\hat\tau_0$ $=$ `r round(tau0,2)` and $\hat\tau_1$ $=$ `r round(tau1,2)`), which suggests to me that the random effects, rather than the fixed effects, are picking up the interesting variation in $p_i(outcome=1)$ and in the treatment effect.  

```{r echo=F}
res <- summary(toeresult)$summary[c("b0","b1","b2","tau0","tau1"),1]
b0 <- res[1]
b1 <- res[2]
b2 <- res[3]
tau0 <- res[4]
tau1 <- res[5]
```
The estimates of the fixed effects from the stan() model ($\hat \beta_0$ $=$ `r round(b0,2)`, $\hat \beta_1$ $=$ `r round(b1,2)` and $\hat \beta_2$ $=$ `r round(b2,2)`) and SD's of random effects ($\hat\tau_0$ $=$ `r round(tau0,2)` and $\hat\tau_1$ $=$ `r round(tau1,2)`) are just much more believable. 


\rm

## Problem 3(e)

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 \texttt{stan()} model, and we probably should since the \texttt{stan()} results seem more plausible than the \texttt{glmer()} results, but \texttt{DHARMa} is faster and less prone to user errors; it will be fine for the purpose of just getting used to using and interpreting \texttt{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).

```{r, fig.height=5}
plot(glmer.1)

plot(simres <- simulateResiduals(glmer.1))
```
\it
The plot of the conditional residuals is virtually useless, just as we expect for untransformed residuals from a generalized linear model.  

The plots of the uniformized residuals tell a somewhat different story.  

In the qqplot on the left, there is some indication that may be a little long-tailed to the right (the plot curves upward from a straight line), the KS test for whether the "uniformized" residuals really follow the normal distribution is has a very low p-value (suggesting that the transformed data doesn't really follow the normal distribution, as it would if the model fit well), and the dispersion test also has a low p-value, indicating that the residuals are over-dispersed (this over-dispersion could be related to the 0-inflation and 1-inflation we discovered with posterior predictive checks [ppc's]).  

In the residuals vs fitted plot on the right, the smooth curves that are supposed to be estimating the horizontal lines at the 25$^{th}$, 50$^{th}$ and 75$^{th}$ percentiles, seem to be increasing (rather than horizontal) over most of their range.  This too is consistent with what we discovered with ppc's.

\rm

