---
title: "Bayesian Model Averaging: NELS and SAT"
author: "Carolyn J. Anderson"
date: '2022-09-29'
output: html_document
---
What do you do when there is uncertainity about the model?


A short explanation starts on page 2 of  https://cran.r-project.org/doc/Rnews/Rnews_2005-2.pdf and it talks specifically abou the bma package.

A very accessible and entertaining paper Hinne M, Gronau QF, van den Bergh D, Wagenmakers E-J. A Conceptual Introduction to Bayesian Model Averaging. Advances in Methods and Practices in Psychological Science. 2020;3(2):200-215. doi:10.1177/2515245919898657.  I saw this as a presentation a couple of year ago..."no demon every completely vanishes".

The comments below are limited and documentation should be consulted!

# NELs
```{r}
library(BMA)    # To do the model averaging

setwd("C:/Users/cja/Dropbox/edps 590BAY/Lectures/7 Multiple Linear Regression")

df <- read.table("nels_school_62821.txt", header=TRUE)
names(df)
head(df)

y <- df$math
x <- subset(df, select=c(sex, homework, ses, paredu))
```

Example using linear regression (glms and other models are possible.
```{r}
nels.bma <- bicreg(x,y)
summary(nels.bma)
```

Model plots: 1st color=positive coefficient, 
             2nd color=negative,
             3rd color=not included in model
```{r}
imageplot.bma(nels.bma, color=c("blue","red","grey"))
```
Density plots: Displays plots of the posterior distributions of the coefficients generated by Bayesian model averaging over linear regressions.

The posterior distribution of coefficients (model averaged) where spike at 0 shows posterior probability that the variable is NOT in the model.

```{r}
plot(nels.bma)
```

# Getting what you paid for

Another example....the data set we've been looking at all semester.

```{r}
get <- read.table("getting_what_paid_for_data.txt",header=TRUE)

names(get)

y <- get$ave_tot
x <- subset(get, select=c("exp_pp","ave_pt", "salary","taking"))

get.bma <- bicreg(x,y)
summary(get.bma)

# Model plots
imageplot.bma(get.bma, color=c("blue","red","grey"))

# Density plots
plot(get.bma)