---
title: "Nels: Mean and Variance using rstan"
author: "Carolyn J. Anderson"
date: "10/29/2021"
output:
   html_document:
     toc:  yes
---
# Overview

I will create a series of documents that run various analysis on the nels 23 schools data. I could put them all in one document, but this would take a really, really long time. The documents in the set will do the following

* Estimate mean and sd of math scores
* Fit math ~ homework
* Fit math  ~ homework + ses
* Fit HLM:   math ~ homework + (1 |school)
* Fit HLM:   math ~ homework + (1 + homework|school)


# Set-up

## Packages 
So, let's load up the rstan package at the start and see the messages that we're given:
```{r}
library(rstan)
library(stargazer)   # use for table of descriptive
```
This package computes some useful statistics for model evaluation and comparison. One statistic is efficient approximate "leave one out" (hence it's called "loo") cross-validation for Bayesian models. The loo package will also compute waic ("widely applicable information criteria") and compare waic for multiple models using the  loo_compare function.
```{r}
library(loo)
```

## Recommendation for rstan
Note that there three recommendations here about what to run before we use rstan.  The first one, "For execution on a local, multicore CPU with excess RAM we recommend calling", only needs to be issue once per session.  This make the default number of core to use equal to the number of cores available
```{r}
options(mc.cores = parallel::detectCores())
```
The second recommendation, "To avoid recompilation of unchanged Stan programs, we recommend calling", save the compile model so that it can be re-used.  Models are compiled in C++ and this can take some time.
```{r}
rstan_options(auto_write = TRUE)
```
The third recommendation, "For improved execution time, we recommend calling", I'm not exactly sure what this does.  "For improved execution time, we recommend calling"
```{r}
Sys.setenv(LOCAL_CPPFLAGS = '-march=native -mtune=native')
```

## Directory and Data Set Up
My working directory where the data live is 
```{r}
setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/10 Ham_Stan_brms")
```

And read in data
```{r}
nels <- read.table("school23_data.txt",header=TRUE)
nels <- nels[order(nels$school,nels$student), ]
# Get information on numbers of school and students
school <- unique(nels$school)               
N <- length(school)                        
# total number of students
n <- length(nels$math)               
names(nels)
```
I like computing sample statistics 
```{r  gazer}
stargazer(nels, type="text") # default is LaTex
```

#  Estimate mean and sd of math scores
Our first stan model!  We'll start simple and work up to more complex ones.  

We could either define the model as below and run in rstan like we did with rjags OR we could put this in a file with extension .stan 
```{r}
#############################################
# stan model definition
#############################################
model1 <- "
data {                // We need to define the data
  int<lower=0> n;     // number of students
  real y[n];          // vector of data
} 
parameters {
  real<lower=0> sigma; 
  real mu;             
} 
model {
  target += cauchy_lpdf(sigma | 0, 1);
  target += normal_lpdf(mu |0, 30);
  target += normal_lpdf(y | mu, sigma);
}
generated quantities{
  real log_lik[n];
  for (i in 1:n) {
    log_lik[i] = normal_lpdf(y[i] | mu, sigma); 
          }
}
"
```
"lpdf" stands for log probability density function.

Alternatively the model part could be written as 
```{r eval=FALSE}
   sigma ~ cauchy(0,1);  
   mu ~ normal(0, 10);
   y ~ normal(mu,sigma);
}
```
"lpdf" stands for log probability density function.


Before we run rstan, we need to create a data list, as follows:
```{r}
########################################
# Make data list
########################################
y <- as.numeric(nels$math) # otherwise it's class integer
n <- nrow(nels) 
dataList <- list(n=n,y=y) 
```
Now for the part that takes a bit of time
```{r fit1}
fit.1 <- stan(
      model_code = model1, 
		  model_name = "Nels23: mean & sd",
      data = dataList, 
		  iter = 2000, 
		  chains = 4,
		  cores=4,
		  warmup = floor(2000/2),
		  verbose = FALSE) 
```

And some output but I don't want all the log likelihoods here, so I'll only ask for mean and standard deviation:

```{r}
print(fit.1,pars = c("mu","sigma"), probs = c(0.025, 0.975), 
    digits = 4)
```
and various graphics
```{r convergencechk}
plot(fit.1)
stan_trace(fit.1)
stan_dens(fit.1)
stan_hist(fit.1)
stan_ac(fit.1)
stan_scat(fit.1,pars=c("mu","sigma"))
pairs(fit.1, pars = c("mu", "sigma", "lp__"), las = 1)
```

There are some more that we could look at but I'm not exactly sure what to do with these
```{r eval=FALSE}
stan_diag(fit.2,information=c("stepsize"))
stan_diag(fit.2,information=c("divergence"))
stan_diag(fit.2,information=c("treedepth"))
```
And to use the loo package,
```{r}
log_lik.1 <- extract_log_lik(fit.1, merge_chains = FALSE)
r_eff <- relative_eff(exp(log_lik.1)) 

loo.1 <- loo(log_lik.1, r_eff = r_eff, cores = 4)
print(loo.1)

waic.1 <- waic(log_lik.1)
print(waic.1)
```

If we had two models we might want to look at 
```{r eval=FALSE}
compare(loo.1,loo.2)
```