---
title: "NELS 23 HLM via brms and rstan"
author: "Carolyn J. Anderson"
date: "11/4/2021"
output:
   html_document:
      toc: yes
---

```{r homedesktop}
#.libPaths("C:/Program Files/R/R-4.1.1/library")
```

# Set up
To make this run but will drop when put into master set of examples.
```{r}
library(rstan)
library(loo)
library(brms)
library(coda)

setwd("D:/Dropbox/edps 590BAY/Lectures/10 Ham_Stan_brms")

nels <- read.table("school23_data.txt",header=TRUE)
names(nels)

# sort by schools (may speed things up for large data sets):
nels <- nels[order(nels$school,nels$student) , ]


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


# Data

Need information on number schools & create school.id
```{r data}
school <- unique(nels$school)                         
N <- length(school)                                   
school.int <- as.data.frame(cbind(school,seq(1:N)))     
names(school.int) <- c("school", "school.id")        
nels <- merge(nels,school.int,by="school")           

# total number of students
n <- length(nels$math)                    
```

# Random Intercept (null/empty HLM)

## brm
```{r, brms1}
hlm1 <- brm(math ~ 1 + (1 | school.id),
            data=nels,
            warmup=2000,
            iter=4000,
            chains=4,
            cores=4,
            seed=11121
            )
```

## Stan data and code
Corresponding stan data and code
```{r hlm1stan}
dataList <- standata(hlm1)

hlm1.stan <- stancode(hlm1)
```

This provide trace plots and densities,
```{r convergencechk1}
plot(hlm1)             
```

If you want the diagnostics from coda (ones we used with Gibbs sampling), you need to extract the samples from brm as an mcmc object.  Note that as.mcmc is deprecated, but I haven't found the function that can be used instead.

```{r hlmcoda}
hlm1.mcmc <- as.mcmc(hlm1)
                     
traceplot(hlm1.mcmc)       
densplot(hlm1.mcmc)
gelman.diag(hlm1.mcmc)
gelman.plot(hlm1.mcmc)
autocorr.plot(hlm1.mcmc)
geweke.diag(hlm1.mcmc)
effectiveSize(hlm1.mcmc)
```

And if all looks good
```{r}
summary(hlm1)
```

# Random Intercept and Slope

```{r, brms2}
hlm2 <- brm(math ~ 1 + homew + ses + (1 + homew| school.id),
            data=nels,
            warmup=2000,
            iter=4000,
            chains=4,
            cores=4,
            seed=11121
            )
```

```{r}
plot(hlm2)              
summary(hlm2)

plot(conditional_effects(hlm2), ask = FALSE)

#or
# plot(conditional_effects(hlm2), points=TRUE, ask = FALSE)

pp_check(hlm2)  # shows dens_overlay plot by default
```
There are many possible posterior predictive checks.  I suggest checking the brms manual.  Just some of them are below.

``{r eval=FALSE}
pp_check(hlm2, type = "error_hist", ndraws = 11)
pp_check(hlm2, type = "scatter_avg", ndraws = 100)
pp_check(hlm2, type = "stat_2d")
pp_check(hlm2, type = "loo_pit")
```
# Level 2: Multiple Predictors, Factors and Cross-level Interactions

brm handles factors like lm or glm (i.e., dummy codes).  We need to declare the variables as factors.
```{r makefactors}
nels$urban <- as.factor(nels$urban)
nels$schtyp <- as.factor(nels$schtype)
nels$classtr <- as.factor(nels$classtr)
```

Then just are above,
```{r, brms3}
hlm3 <- brm(math ~ 1 + homew + ses +  urban + ses*urban + (1 + homew| school.id),
            data=nels,
            warmup=2000,
            iter=4000,
            chains=4,
            cores=4,
            seed=11121
            )
```

```{r}
plot(hlm3)              # trace plots & density
```

```{r}
summary(hlm3)
plot(conditional_effects(hlm3), ask = FALSE)

pp_check(hlm3)  # shows dens_overlay plot by default
```


real normal_lccdf(reals y | reals mu, reals sigma)

The log of the complementary cumulative normal distribution of y given location mu and scale sigma; normal_lccdf will overflow to 0 for y−μσ below -37.5 and underflow to −∞ for y−μσ above 8.25; see above for discussion of Phi_approx as an alternative.


# Specifying QR Decomposition 

By default centering is used.  If you want QR decomposition (which is especially good for correlated predictors), this can be requested using brmsformula.  I am including this because we talked about it for Stan.

The model formula can be specified separately and then put into the brm call. There are a number of options and two that I added here are center (which is TRUE by default) and QR decomposition.  

```{r}
# Options for the forumla
fm3 <- brmsformula(math ~ 1 + homew + ses +  urban + ses*urban + (1 + homew| school.id),
          center=TRUE,
          decomp="QR"
          )
```

The rest is repeating what we did above, except using fm3 as the formula.
```{r, brms3alt}
hlm3.alt <- brm(fm3,
            data=nels,
            warmup=2000,
            iter=4000,
            chains=4,
            cores=4,
            seed=11121
            )
```

```{r}
plot(hlm3.alt)              # trace plots & density
```

```{r}
summary(hlm3.alt)
plot(conditional_effects(hlm3.alt), ask = FALSE)
```

# Comparing Models

In terms of criteria to compare models in order worst to best are Bayes Factor, DIC, waic, loo, and psis-loo.

The criterion "waic" (Widely accessible information criterion) is asymptotically equivalent to loo (Leave-one-out-validation), but loo is faster and easier to compute.  Some work has shown that loo is the preferred for model comparison.  For these data and models, R throws warnings when I ask for "waic" and suggests using loo instead..  

loo is noisy and a solution suggested for this is to use Pareto smoothed importance sampling LOO or PSIS-loo  (Verhari, Gelman & Gabry, "Practical Bayesian model evaluation using leave-one-out cross-validation, http://www.stat.columbia.edu/~gelman/research/unpublished/loo_stan.pdf)  
Since I didn't get is at first, I can add this.  The command below computes PSIS-LOO and compares values for our 3 models.

For more information a good starting place is https://mc-stan.org/loo/reference/loo-glossary.html


Smaller values of looic are better (just like other information criteria), and for elpd larger is better. Note that looic equals −2∗elpd

```{r loo}
hlm1 <- add_criterion(hlm1, criterion="loo")
hlm1$criteria$loo

hlm2 <- add_criterion(hlm2, criterion="loo")
hlm2$criteria$loo

hlm3 <- add_criterion(hlm3, criterion="loo")
hlm3$criteria$loo
```

If you wan to compare, this is a better
```{r loocompare}
loo(hlm1, hlm2, hlm3, compare=TRUE, moment_match=TRUE, save_psis=TRUE)
```

Note that 

*      elpd_loo is the Expected Log Pointwise Predictive Density, which is the sum of N individual pointwise log predictive densities.  
*      p_loo is the defective number of parameters (i.e., the difference between elpd_loo and the non-cross validity log posterior predictive density).
*      elpd_diff is the difference in elpl_loo for two models.  If more than 2 models are compared then it is the difference relative to the model with the highest elpd_loo.   


For more information a good starting place is https://mc-stan.org/loo/reference/loo-glossary.html

# Changing Stan code
Although most of the time you can modify your models within brms, you might want to modify the Stan code and run it using rstan (e.g., a simulation study where only data changes, or this is a starting point for your own model).  Also this is a good way to learn Stan coding.   

For the data input to Stan by brms,
```{r datalist}
data_1 <- standata(hlm1)
```

and for the model code,
```{r stancode1}
stan_code <- stancode(hlm1)
stan_code
```

# Change Priors

If you don't want to look at Stan code, you may want to know what the priors are that were used to fit a model.  For example,

```{r prior2}
prior_summary(hlm2)
```

Note that for sigma this is actually a half non-central t-distribution.  Other distributions can be put here. Gelman recommends a half Cauchy distribution (a truncated Cauchy, 

$$ f(y| \mu, \sigma) = \left\{ \begin{array}{ll}
   \frac{2}{\pi\sigma} \frac{1}{1+(y-\mu)^2/\sigma^2} & y \ge \mu \\
   0 & \mbox{otherwise} \end{array}\right. $$


How to change your priors?
You can define priors and ask brms to create the Stan code.  

```{r setallb}
prior1 <- prior(normal(0, 5), class = "Intercept")
prior2 <- prior(cauchy(0, 1), class = "sd")
hlm1priors <- c(prior1, prior2)

hlm1.stan <- make_stancode(math ~ 1 + (1 | school.id), data=nels, prior=hlm1priors)

data1 <-standata(hlm1)
```

```{r rstan1}
hlm1.rstan <- stan(model_code = hlm1.stan, 
	    model_name = "Nels23: math ~ 1 + (1 |school.id)",
      data = data1,
	    iter = 4000,
	    chains = 4,
	    cores = 4,
	    verbose = FALSE
	   ) 
```

Based on the following diagnostics, my changes make things pretty bad. Some of the chains are not mixing well, some large auto-correlations, and strange looking posteriors.  

```{r lastone}
stan_trace(hlm1.rstan)

stan_dens(hlm1.rstan)

stan_hist(hlm1.rstan)

stan_ac(hlm1.rstan)

print(hlm1.rstan)
```