---
title: "587 Workup of Bayesian Example"
author: "Carolyn J. Anderson"
date: "3/2/2022"
output: 
  html_document:
    toc: yes
---

In this document are examples on how to use brms. We will use the nels data for $N=23$.   We start with running code to get everything set up for model fitting. I have included libraries in case you want to also do maximum likelihood estimation of the HLM models.

# Bayesian Estimation

$$ f(\mathbf{\gamma},\mathbf{T},\sigma^2 | \mathbf{data}) = \frac{ f(\mathbf{data}\ | \ \mathbf{\gamma},\mathbf{T},\sigma^2)f(\mathbf{\gamma},\mathbf{T},\sigma^2) 
}{ f(\mathbf{data})}$$

where 


*   Posterior (target) distribution is $f(\mathbf{prameters}\ | \ \mathbf{data})$.  We want this!
*   Likelihood (or data model) is $f(\mathbf{data}\ | \ \mathbf{parameters})$.  This is distribution for data given parameters.
*   The prior distribution is $f(\mathbf{parameters})$.  This reflects prior knowledge.
*   The data distribution $f(\mathbf{data}$


We will use MCMC (Markov Chain Monte Carlo to sample from the posterior), but we'll just work with 
$$f(\mathbf{parameters}\ | \ \mathbf{data}) \propto f(\mathbf{data}\ |  \ \mathbf{parameters})f(\mathbf{\gamma},\mathbf{T},\sigma^2)  $$
A link to see what we'll be doing.   

https://chi-feng.github.io/mcmc-demo/app.html

Stan (brms) uses Hamliltonian sampling.

# The brms package

The brms package is a wrapper or interface that make fitting mixed (random effects) models, which makes using Stan relatively easy.  The hardest part is installing Stan and  brms. This gives you enough so that you can fit models using Bayesian estimation but not enough so that you understand what's going on.  If you really want to use Bayesian, you should learn more about it.  

The brms package will create input that is  fed into Stan. This is easier than writing your own input and then using the rstan package, which is another way to run Stan (you can get the Stan input from the brms).



# From brms documentation:  

"Fit Bayesian generalized (non-)linear multivariate multilevel models using 'Stan' for full Bayesian inference. A wide range of distributions and link functions are supported, allowing users to fit -- among others -- linear, robust linear, count data, survival, response times, ordinal, zero-inflated, hurdle, and even self-defined mixture models all in a multilevel context. Further modeling options include non-linear and smooth terms, auto-correlation structures, censored data, meta-analytic standard errors, and quite a few more. In addition, all parameters of the response distribution can be predicted in order to perform distributional regression. Prior specifications are flexible and explicitly encourage users to apply prior distributions that actually reflect their beliefs. Model fit can easily be assessed and compared with posterior predictive     checks and leave-one-out cross-validation. References: B?rkner (2017)    <doi:10.18637/jss.v080.i01>; B?rkner (2018) <doi:10.32614/RJ-2018-017>;
    Carpenter et al. (2017) <doi:10.18637/jss.v076.i01>."

# Set up
```{r libraries, message=FALSE, warning=FALSE, include=TRUE}
library(brms)
library(rstan)
library(lme4)
library(lmerTest)
library(coda)
library(mvtnorm)
library(lattice)
library(optimx)
```


Recommended options for use with stan:
```{r RecommendedStanOptions}
# for stan
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
```

Read in the data and do a bit of set-up
```{r readData}
#setwd("D:/Dropbox/edps 590BAY/Lectures/10 Ham_Stan_brms")
setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/10 Ham_Stan_brms")
nels <- read.table("school23_data.txt",header=TRUE)
head(nels)
tail(nels)

# sort by schools:
nels <- nels[order(nels$school,nels$student) , ]

# Get information on number schools & create a variable "school.id" which has consecutive integers as Id for the schools.  
school <- unique(nels$school)                        ##  
N <- length(school)                                  ## need N and school.id
school.int <- as.data.frame(cbind(school,seq(1:N)))  ##  for dataList & model  
names(school.int) <- c("school", "school.id")        ##
nels <- merge(nels,school.int,by="school")           ##

# total number of students
n <- length(nels$math)                               ## need n for dataList & model

# You don't need this but might want it
nj <- matrix(999,nrow=N,ncol=1)
for (j in 1:N) {
  nj[j] <- nrow(subset(nels,nels$school.id==j))
}  
```

# Simple brms for mean and variance
We'll start simple (fastest of the bunch) and use default priors,

```{r MeanVariance}
options(mc.cores = parallel::detectCores())
t0<-Sys.time()

summary(model.simple <- brm(math~ 1  , data=nels))

(tfinal<-Sys.time()-t0)
plot(model.simple)
```
We'll talk about the output in lecture.  


Note that this take a while and a considerable amount of time is spent getting the model into C++ (i.e., "Compiling Stan program").


More diagnostic plots but using coda
```{r codaplots}
# Change class of model output
simple.mcmc <- as.mcmc(model.simple)

traceplot(simple.mcmc)
densplot(simple.mcmc)
gelman.diag(simple.mcmc)
geweke.diag(simple.mcmc, frac1=0.1, frac2=0.5)
autocorr.plot(simple.mcmc)
```

We can obtain posterior predictive checks see brms documentation for different types. Here are some examples (they might look better with a more reasonable model)
```{r PostPredChecks}
pp_check(model.simple)

pp_check(model.simple, type="error_hist")

pp_check(model.simple, type="ecdf_overlay")
```


## Stan input produced by brms

Before going on to more interesting models, we'll take a look at what brms is doing for you.  This can be useful for  

*     Simulation studies (faster ones)
*     Continue to learn the Stan language

To see the Stan code that brms produced 
```{r stanSimpleCode}
(simple.stan <-stancode(model.simple))
```
Take a look at this and note what the default priors are.   

We can also get the data list used with this stan code by using 
```{r stanDataList}
dataList0 <- standata(model.simple)
```

Every time we want to or need to run this model, we would have to recompile the model; however, we can take the Stan code and create a compiled version of it.  If we have such a compile version, we can re-use the model with another data set.  This would be a great time saver if you're running a simulation study or you want to fit the model to different data sets.  For, example I could use if for the anorexia data (change data list of course).

I made a few changes that I though might speed things up; in particular, getting ride of the "if" statement.  I also got ride of the temporary intercept.  You can also changes the priors (if you want or could set them before running brms using the "set_prior" command.   

To compile the model without sampling use the command "stan_model",
```{r savestanModelCode}
m.0 <- stan_model(model_code = "data {
	int<lower=1> N;  // number of observations
	vector[N] Y;  // response variable
	int prior_only;  // should the likelihood be ignored?
	}
	transformed data {
	}
	parameters {
	real b_Intercept;  // temporary intercept
	real<lower=0> sigma;  // residual SD
	}
	transformed parameters {
	}
	model {
	vector[N] mu = b_Intercept + rep_vector(0, N);
	// priors including all constants
	target += student_t_lpdf(b_Intercept | 3, 51, 13);
	target += student_t_lpdf(sigma | 3, 0, 13)
		- 1 * student_t_lccdf(0 | 3, 0, 13);
	// likelihood including all constants
		target += normal_lpdf(Y | mu, sigma);
	}
	"
)
```

I can now use and re-use this model, using the command "sampling" and include all the options I would normaly include when running stan.  For the nels data,
```{r sampingusingStan}
# Now do sampling using the nels data list
model.0 <- sampling(m.0, 
           data=dataList0,
					 chains=4,
			    core=4,
					)
stan_plot(model.0)
```

Notice that sampling is fast compared to creating the C++ version of our model.


If you want more for diagnostic graphics from Stan
```{r stanplots, eval=FALSE}
stan_trace(model.0)   # trace plots
stan_dens(model.0)    # density plots for specified parameters
stan_hist(model.0)    # histograms
stan_ac(model.0)      # auto-correlation
```


# Random intercept model

## Random intercept via lmer
Using MLE, because you are familiar with this.

```{r RIlmer}
modelRI.lmer <- lmer(math~ 1 + homew + (1 |school.id), data=nels, REML=TRUE, 
                    control = lmerControl(optimizer ="Nelder_Mead"))
summary(modelRI.lmer)
```

## Random intercept via brms
And using brms, the model formula is the same,
```{r RIbrm}
modelRI.brm <- brm(math ~ homew + (1 | school.id),
         data=nels,
         family=brmsfamily("gaussian"),
         save_all_pars=TRUE,
         cores=4, 
         silent=FALSE)
```
### Some diagnostics for Bayesian
```{r RIdiag}
plot(modelRI.brm)

modelRI <- as.mcmc(modelRI.brm)

traceplot(modelRI)
densplot(modelRI)
geweke.diag(modelRI, frac1=0.1, frac2=0.5)
autocorr.plot(modelRI)
```

### Summary table for Random Intercept
```{r summaryRI}
summary(modelRI.brm)
```

# More complex model


And lets make this more complex by adding random slope, and some more fixed effects including a cross level interaction,
```{r}
# Random intercept & slope with level 1 and level 2 predictors    

nels$public <- ifelse(nels$schtype==1,1,0)
nels$public.fac <- as.factor(nels$public)
```

##  via lmer
```{r rislmer}
modelRIS.lmer <- lmer(math~ 1 + homew + ses + public.fac  + homew*public.fac + (1 + homew|school.id),
                    data=nels, REML=TRUE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(modelRIS.lmer)
```

## via brms
And now for the Bayesian model,
```{r RISbrm}
modelRIS.brm <- brm(math~ 1 + homew + ses + public.fac  + homew*public.fac + (1 + homew|school.id), 
                data=nels, cores=4,  save_all_pars=TRUE)  # fit the model
summary(modelRIS.brm)                  # summary of model information
plot(modelRIS.brm)                     # trace and density plots
bayes_R2(modelRIS.brm)                 # Bayes version of R2
bayes_factor(modelRIS.brm, modelRI.brm)  # Computes bayes factor     
```

# Comparing Models

Possible criteria for comparing models are (in order worst to best)  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.


We can pull these out for each model,
```{r loo}
modelRI.brm <- add_criterion(modelRI.brm, criterion="loo")
modelRI.brm$criteria$loo

modelRIS.brm <- add_criterion(modelRIS.brm, criterion="loo")
modelRIS.brm$criteria$loo
```

But better yet, for a direct comparison use
```{r compareall}
loo(modelRI.brm, modelRIS.brm, 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 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.   
Smaller values of looic are better (just like other information criteria), and for elpd larger is better.  Note that looic equals $-2*elpd$.


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

# Other model families

There are lots!!  We specify them in brms pretty much how we would do these using glm; however, the range of models is greater with brms than with glm.  The model families (from the brms documentation 8/29/2019) include

*     student(link = "identity", link_sigma = "log", link_nu = "logm1")
*     bernoulli(link = "logit")
*     negbinomial(link = "log", link_shape = "log")
*     geometric(link = "log")
*     lognormal(link = "identity", link_sigma = "log")
*     shifted_lognormal(link = "identity", link_sigma = "log",
*     link_ndt = "log")skew_normal(link = "identity", link_sigma = "log",link_alpha = "identity")
*     exponential(link = "log")weibull(link = "log", link_shape = "log")
*     frechet(link = "log", link_nu = "logm1")
*     gen_extreme_value(link = "identity", link_sigma = "log",link_xi = "log1p")
*     exgaussian(link = "identity", link_sigma = "log", link_beta = "log")
*     wiener(link = "identity", link_bs = "log", link_ndt = "log",link_bias = "logit")
*     Beta(link = "logit", link_phi = "log")
*     dirichlet(link = "logit", link_phi = "log", refcat = NULL)
*     von_mises(link = "tan_half", link_kappa = "log")
*     asym_laplace(link = "identity", link_sigma = "log",link_quantile = "logit")
*     hurdle_poisson(link = "log")
*     hurdle_negbinomial(link = "log", link_shape = "log",link_hu = "logit")
*     hurdle_gamma(link = "log", link_shape = "log", link_hu = "logit")
*     hurdle_lognormal(link = "identity", link_sigma = "log",link_hu = "logit")
*     zero_inflated_beta(link = "logit", link_phi = "log",link_zi = "logit")
*     zero_one_inflated_beta(link = "logit", link_phi = "log",link_zoi = "logit", link_coi = "logit")
*     zero_inflated_poisson(link = "log", link_zi = "logit")
*     zero_inflated_negbinomial(link = "log", link_shape = "log",link_zi = "logit")
*     zero_inflated_binomial(link = "logit", link_zi = "logit")
*     categorical(link = "logit", refcat = NULL)
*     multinomial(link = "logit", refcat = NULL)
*     cumulative(link = "logit", link_disc = "log",threshold = c("flexible", "equidistant"))
*     sratio(link = "logit", link_disc = "log", threshold = c("flexible","equidistant"))
*     cratio(link = "logit", link_disc = "log", threshold = c("flexible","equidistant"))
*     acat(link = "logit", link_disc = "log", threshold = c("flexible","equidistant"))



