Yes. Sorry, here is the entire code (as a markdown)

title: “Bayesian Stacking”

date: “Oct. 11, 2020”

## output: html_document

```
knitr::opts_chunk$set(echo = TRUE)
```

```
library(rstanarm)
library(loo)
library(LaplacesDemon)
library(BMS)
```

```
PISAdata <- read.csv("pisa.2018.usa.imp.read.csv", header=T)
PISAdata <- PISAdata[1:4000,]
PISAdata$Female <- ifelse(PISAdata$Female=="Female",1,0)
```

# Four models of interest to stack.

Model 1: Baseline model

Model 2: FEMALE, ESCS, HOMEPOS, ICTRES

Model 3: JOYREAD, PISADIFF, SCREADCOMP, SCREADDIF

Model 4: METASUM, GFOFAIL, MASTGOAL, SWBP, WORKMAST, ADAPTIVITY, COMPETE

Model 5: PERFEED, ADAPTIVE, TEACHINT, BELONG

```
attach(PISAdata)
fits <- list()
fits[[1]] <- stan_glm(PV1READ ~ 1, data = PISAdata, seed = 832762) # This is the baseline model 1
fits[[2]] <- update(fits[[1]], formula = PV1READ ~ Female + ESCS + HOMEPOS + ICTRES) # model 1 variables
fits[[3]] <- update(fits[[1]], formula = PV1READ ~ JOYREAD + PISADIFF + SCREADDIFF + SCREADCOMP) # model 2 variables
fits[[4]] <- update(fits[[1]], formula = PV1READ ~ METASUM + GFOFAIL + MASTGOAL + SWBP + WORKMAST + ADAPTIVITY + COMPETE) # model 3 variables
fits[[5]] <- update(fits[[1]], formula = PV1READ ~ PERFEED + TEACHINT + BELONG) # model 3 variables
```

# Compute loo

```
loo_list <- lapply(fits, loo, cores = 4)
```

# Compute stacking weights. The first model is the intercept

```
wtsStacking <- loo_model_weights(loo_list, method = "stacking")
print(wtsStacking)
```

Notice that the weight is mostly given to model 3.

# Average predictions (sample with probability equal to weights)

```
n_draws <- nrow(as.matrix(fits[[1]]))
ypredStacking <- matrix(NA, nrow = n_draws, ncol = nobs(fits[[1]]))
for (d in 1:n_draws) {
k <- sample(1:length(wtsStacking), size = 1, prob = wtsStacking)
ypredStacking[d, ] <- posterior_predict(fits[[k]], draws = 1)
}
```

# This provides the distribution of predicted values and compares them to the actual values using the KLD.

```
ypredStacking <- colMeans(ypredStacking)
#kld1 <- KLD(ypredStacking,PV1READ)$sum.KLD.py.px
```

```
KLD(ypredStacking,PV1READ)$sum.KLD.py.px
```