11.1 Extending the Normal Regression Model

# Load some packages
library(bayesrules)
library(rstanarm)
library(bayesplot)
library(tidyverse)
library(broom.mixed)
library(tidybayes)

# Load the data
data(weather_WU)
weather_WU %>% 
  count(location)

## # A tibble: 2 × 2
##   location       n
##   <fct>      <int>
## 1 Uluru        100
## 2 Wollongong   100

weather_WU <- weather_WU %>% 
  select(location, 
         windspeed9am, humidity9am, pressure9am, 
         temp9am, temp3pm)

weather_WU%>%head

## # A tibble: 6 × 6
##   location windspeed9am humidity9am pressure9am temp9am temp3pm
##   <fct>           <dbl>       <int>       <dbl>   <dbl>   <dbl>
## 1 Uluru              20          23       1023.    20.9    29.7
## 2 Uluru               9          71       1013.    23.4    33.9
## 3 Uluru               7          15       1012.    24.1    39.7
## 4 Uluru              28          29       1016     26.4    34.2
## 5 Uluru              24          10       1010.    36.7    43.3
## 6 Uluru              22          32       1012.    25.1    33.5

A typical afternoon temperature is somewhere between 15 and 35 degrees Celsius.

ggplot(weather_WU, 
       aes(x = temp9am, y = temp3pm)) +
  geom_point(size = 0.2)+
  geom_smooth(method="lm")

## `geom_smooth()` using formula = 'y ~ x'

(15+35)/2

## [1] 25

\[\beta_{c0}\sim N(25,5^2)\]

mean(weather_WU$temp3pm)

## [1] 24.5585

\[\beta_{1}\sim N(0,3.1^2)\]

mean(weather_WU$temp9am);

## [1] 19.5405

sd(weather_WU$temp9am);

## [1] 6.051522

weather_model_1 <- stan_glm(
  temp3pm ~ temp9am, 
  data = weather_WU, 
  family = gaussian,
  prior_intercept = normal(25, 5),
  prior = normal(0, 2.5, autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4, iter = 5000*2, seed = 84735)

# Prior specification
prior_summary(weather_model_1)

## Priors for model 'weather_model_1' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = 25, scale = 5)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 3.1)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.13)
## ------
## See help('prior_summary.stanreg') for more details

# MCMC diagnostics
mcmc_trace(weather_model_1, size = 0.1)

mcmc_dens_overlay(weather_model_1)

mcmc_acf(weather_model_1)

neff_ratio(weather_model_1)

## (Intercept)     temp9am       sigma 
##     0.99710     1.00255     0.94530

rhat(weather_model_1)

## (Intercept)     temp9am       sigma 
##   0.9999352   0.9999668   0.9999059

posterior_interval(weather_model_1, prob = 0.80)

##                   10%      90%
## (Intercept) 2.9498083 5.448752
## temp9am     0.9802648 1.102423
## sigma       3.8739305 4.409474

pp_check(weather_model_1)