---
title: "NELS23 Multiple Regression via rstan"
author: "Carolyn J. Anderson"
date: "10/18/2022"
output:
   html_document:
      toc: yes
---


To make this run but will drop when put into master set of examples.
```{r}
library(rstan)
library(stargazer)
```
As rstan recommends
```{r}
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
```
Get the data
```{r}
setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/10 Ham_Stan_brms")

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

# Model defintion for Multiple Regression

Rather than list out all the regression parameters, we can make use vectors (regression coefficients) and matrices (design matrix which contains the explanatory/predictor variables)
```{r}
#############################################
# stan model definition                     #
#############################################
mlreg <- "
data {
  int<lower=0> N;   // number of students
  int<lower=0> K;   // number of predictors
  matrix[N,K] x;   // predictor matrix
  vector[N] y;      // outcome vector
} 
parameters {
  real b0;              // intercept
  vector[K] beta;       // coefficients for predictors
  real<lower=0> sigma;  // error scale
} 
model {
  target += normal_lpdf(y | b0 + x*beta, sigma);
} 
"
```


In the model you need to make sure that matrices are conformable in the specification of the model for y.   

We could have added the intercept to beta and a column of ones in x.

For the data list we need to create the design matrix, x, and include the number of columns of x (i.e., number of predictors and coefficients for the parameter vector beta).  x should be of type matrix.

```{r}
N <- nrow(nels) 
y = nels$math
x <- matrix(cbind(nels$homework,nels$ses),nrow=N)
K <- ncol(x)
dat <- list(N = N, K=K, y = y, x = x); 
```

And to run rstan, nothing really new here.

```{r}
#########################################
# Get estimates                         #
#########################################
start.time <- Sys.time()
fitLRsimple <- stan(
     model_code = mlreg , 
	   model_name = "Multiple linear regression", 
     data = dat, 
	   iter = 2000, 
	   chains = 4,
	   cores = 4, 
	   warmup = floor(2000/2),
	   verbose = FALSE, 
	   sample_file = file.path(tempdir(), 'nels_mlreg.txt')) 
end.time <- Sys.time()
end.time-start.time
```

After the model runs, check for convergence and ask for what ever statistics (summaries) and graphics to help you here.

## Output and various graphics

```{r}
print(fitLRsimple, probs=c(.025,.50,.975), digits=2)
```

```{r convergencechk2}
stan_trace(fitLRsimple)
stan_dens(fitLRsimple)
stan_hist(fitLRsimple)
stan_ac(fitLRsimple)
```


And if all looks good
```{r}
print(fitLRsimple)
plot(fitLRsimple)
```

# Center predictors
```{r centering}
nels$chomework <- nels$homework - mean(nels$homework)
nels$cses <- nels$ses - mean(nels$ses)
```
```{r}
N <- nrow(nels) 
y = nels$math
x <- matrix(cbind(nels$chomework,nels$cses),nrow=N)
K <- ncol(x)
cdat <- list(N = N, K=K, y = y, x = x); 
```
```{r Centered}
start.time <- Sys.time()
fitLRcenter <- stan(
     model_code = mlreg , 
	   model_name = "Multiple linear regression", 
     data = cdat, 
	   iter = 2000, 
	   chains = 4,
	   cores = 4, 
	   warmup = floor(2000/2),
	   verbose = FALSE, 
	   sample_file = file.path(tempdir(), 'nels_mlreg.txt')) 
end.time <- Sys.time()
end.time-start.time
```

```{r centerchk}
stan_trace(fitLRcenter)
stan_dens(fitLRcenter)
stan_hist(fitLRcenter)
stan_ac(fitLRcenter)
```


```{r centered}
print(fitLRcenter, probs=c(.025,.50,.975), digits=2)
```

# QR parameterization

see https://mc-stan.org/docs/2_28/stan-users-guide/QR-reparameterization.html

```{r qr}
mlreg3 <- "
data {
  int<lower=0> N;   // number of students
  int<lower=0> K;   // number of predictors
  matrix[N,K] x;   // predictor matrix
  vector[N] y;      // outcome vector
} 
transformed data {
   matrix[N,K] Q_ast;
   matrix[K,K] R_ast;
   matrix[K,K] R_ast_inverse;
   
   Q_ast = qr_thin_Q(x) * sqrt(N-1);   // thin and scale QR decomposition
   R_ast = qr_thin_R(x) / sqrt(N-1);
   R_ast_inverse = inverse(R_ast);
}
parameters {
  real alpha;             
  vector[K] theta;      
  real<lower=0> sigma;  
} 
model {
  y  ~ normal(Q_ast * theta + alpha, sigma);
} 
generated quantities {
  vector[K] beta;
  beta = R_ast_inverse * theta; 
}
"
```

```{r}
#########################################
# model 3L QR                           #
#########################################
start.time <- Sys.time()
fitLR3 <- stan(
           model_code = mlreg3, 
	         model_name = "Multiple linear regression", 
           data = dat, 
	         iter = 2000, 
	         chains = 4,
	         cores = 4, 
	         warmup = floor(2000/2),
	         verbose = FALSE
           ) 
end.time <- Sys.time()
end.time - start.time
```

After the model runs, check for convergence and ask for what ever statistics (summaries) and graphics to help you here.


```{r convergencechk3}
stan_trace(fitLR3)
stan_dens(fitLR3)
stan_hist(fitLR3)
stan_ac(fitLR3)
```


And if all looks good
```{r}
print(fitLR3)
plot(fitLR3)
```

# QR and Centering

We have to center and put centered value in data list.  We could use the rescale function but I prefer this.  

```{r}
N <- nrow(nels) 
y = nels$math
# center  
nels$cx1 <- nels$homework - mean(nels$homework)
nels$cx2 <- nels$ses - mean(nels$ses)
cx <- matrix(cbind(nels$cx1, nels$cx2), nrow=N)
K <- ncol(cx)
dat <- list(N = N, K=K, y = y, x = cx); 
```

By keeping the same name for x we can reuse code.  The model is



```{r qrcentered}
mlreg4 <- "
data {
  int<lower=0> N;   // number of students
  int<lower=0> K;   // number of predictors
  matrix[N,K] x;   // predictor matrix
  vector[N] y;      // outcome vector
} 
transformed data {
   matrix[N,K] Q_ast;
   matrix[K,K] R_ast;
   matrix[K,K] R_ast_inverse;
   
   Q_ast = qr_thin_Q(x) * sqrt(N-1);   // thin and scale QR decomposition
   R_ast = qr_thin_R(x) / sqrt(N-1);
   R_ast_inverse = inverse(R_ast);
}
parameters {
  real alpha;             
  vector[K] theta;      
  real<lower=0> sigma;  
} 
model {
  y  ~ normal(Q_ast * theta + alpha, sigma);
} 
generated quantities {
  vector[K] beta;
  beta = R_ast_inverse * theta; 
}
"
```

Get our samples,

```{r QRandCenterd}
start.time <- Sys.time()
fitLR4 <- stan(
           model_code = mlreg4, 
	         model_name = "Multiple linear regression: Center & QR", 
           data = dat, 
	         iter = 2000, 
	         chains = 4,
	         cores = 4, 
	         warmup = floor(2000/2),
	         verbose = FALSE
           )
end.time <- Sys.time()
end.time-start.time
```

```{r convergencechk4}
stan_trace(fitLR4)
stan_dens(fitLR4)
stan_hist(fitLR4)
stan_ac(fitLR4)
```


And if all looks good
```{r}
print(fitLR4)
plot(fitLR4)
```