---
title: 'Nels23: Linear Regression via Rstan'
author: "Carolyn J. Anderson"
date: "10/27/2021"
output: 
   html_document:
     toc: yes
---

This is the 2nd document in the set and uses the Nels 23 schools to estimate a simple linear regression model.

```{r message=FALSE}
library(rstan)
library(loo)
library(stargazer)
```

```{r}
setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/10 Ham_Stan_brms")
```

```{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)
```
```{r}
stargazer(nels, type="text")
```

```{r}
#############################################
# Simple Linear Regression in stan model 
#############################################
lr.1 <- "
data{
  int<lower=0> n;
  vector[n] x;   
  vector[n] y;   
}
parameters {
  real beta0;
  real beta1;
  real<lower=0> sigma; 
} 
model {
  target += normal_lpdf(beta0 |0, 30);
  target += normal_lpdf(beta1 |0, 30);
  target += cauchy_lpdf(sigma | 0, 1);
  target += normal_lpdf(y | beta0+beta1*x, sigma);

}
"
```

or 
```{r eval=FALSE}
   beta0 ~ normal(0,10); 
   beta1 ~ normal(0,10);
   sigma ~ cauchy(0,1);
   y ~ normal(,sigma);
```

```{r}
#################################
# Data List                     #
#################################
y <- as.numeric(nels$math)
hmwk <- as.numeric(nels$homew)
n <- nrow(nels) 
dataList <- list(n=n,y=y,x=hmwk) 
```

```{r}
#################################
# Run rstan                     #
#################################

fit.lr1 <- stan(
      model_code = lr.1, 
	  model_name = "Nels23: math ~ 1 + homework",
      data = dataList, 
	  iter = 2000,
	  chains = 4,
	  cores = 4,
	  verbose = FALSE) 
```

Some checking the convergence
```{r converchk}
stan_trace(fit.lr1 )
stan_dens(fit.lr1 )
stan_hist(fit.lr1 )
stan_ac(fit.lr1 )
stan_scat(fit.lr1 ,pars=c("beta0","beta1"))
stan_scat(fit.lr1 ,pars=c("beta1","sigma"))
pairs(fit.lr1, pars = c("beta0", "beta1", "sigma", "lp__"), las = 1)
```

And if all looks good
```{r}
print(fit.lr1)
plot(fit.lr1 )
```