posthoc distribution from linear models

predict_intervals_lm.R

#-------------------------------
# Let's explore simulating from lm 
# fits! 
#
# Jarrett Byrnes
# 2020-09-08
#-------------------------------


#some useful libraries
library(arm)
library(palmerpenguins)
library(tidyr)
library(dplyr)

#let's fit a model from the palmer penguins ####
mod <- lm(flipper_length_mm ~ body_mass_g, data = penguins)

#get some simulated coefficients using arm::sim() ####
sims <- sim(mod)

#now the model frame
X <- model.frame(mod)

#now get those simulated fitted values
yhat_arm <- coef(sims) %*% t(model.matrix(mod))

#and make it into a nice tidy data frame
yhat_arm <- as.data.frame(t(yhat_arm)) %>%
  gather(sim, flipper_length_mm)


#what's it look like versus the original data?
ggplot() +
  geom_density(data = penguins, aes(x = flipper_length_mm), 
               lwd=3, color = "blue", alpha = 0.1) +
  geom_density(data = yhat_arm, aes(x = flipper_length_mm, group = sim), 
               alpha = 0.2, lwd=0.1) +
  theme_classic()


#can we make a nice little function? ####
predictInterval_lm <- function(mod, n = 100){
  sims <- arm::sim(mod, n = 100)
  model_mat <- t(model.matrix(mod))
  mod_coefs <- coef(sims)
  
  y_sim <- coef(sims) %*% t(model.matrix(mod))
}



# Wait, does stats::simulate work? ####
yhat_simulate <- simulate(mod, nsim = 100)

yhat_simulate <- yhat_simulate %>%
  gather(sim, flipper_length_mm)


#compare - yhat_arm is in blue
ggplot() +
  geom_density(data = yhat_arm, aes(x = flipper_length_mm, group = sim), 
               alpha = 0.2, lwd=0.1, color = "blue") +
  
  geom_density(data = yhat_simulate, aes(x = flipper_length_mm, group = sim), 
               alpha = 0.2, lwd=0.1) +
  theme_classic()


#no - that was just showing error around fitted values
#what about coreSim? ####
library(coreSim)

#get those simulations, add simulation numbers, and rename
yhat_coresim <- qi_builder(mod, newdata = penguins, ci = 1,
                slim = FALSE, nsim = 100) %>%
  mutate(sim = rep(1:100, 344)) %>%
  rename(flipper_length_mm = qi_) %>%
  filter(!is.na(flipper_length_mm))


#plot it - with coreSim in red
ggplot() +
  geom_density(data = yhat_arm, aes(x = flipper_length_mm, group = sim), 
               alpha = 0.2, lwd=0.1, color = "blue") +
  
  geom_density(data = yhat_simulate, aes(x = flipper_length_mm, group = sim), 
               alpha = 0.2, lwd=0.1) +
  
  geom_density(data = yhat_coresim, aes(x = flipper_length_mm, group = sim), 
               alpha = 0.2, lwd=0.1, color = "red") +
  theme_classic()