Generate fitted simulations from a glm model

Generate simulations from a glm model incorporating either error in fitted error. Simulations explore the possible space of what a model might predict rather than an interval for use in comparison to Bayesian posteriors for non-Bayesian models. The output format and functions draw inspiration from the tidybayes::tidybayes() library and merTools::predictInterval()

# S3 method for glm
add_predicted_sims(
  newdata,
  mod,
  n_sims = 1000,
  seed = NULL,
  weights = 1,
  type = c("predict", "fit", "linpred"),
  ...
)

Arguments

newdata

a data.frame of new data to predict

mod

An lm model to simulate from.

n_sims

number of simulation samples to construct

seed

numeric, optional argument to set seed for simulations

weights

numeric, optional argument for binomial models that need a number of trials

type

Character defining if we are looking at fit or predict intervals.

...

Unused dots for compatibility with generic functions.

Value

A tibble::tibble with information about simulate values.

See also

Other glm: add_fitted_sims.glm()

Examples

# Gamma
clotting <- data.frame(
   u = c(5,10,15,20,30,40,60,NA,100),
   lot1 = c(118,58,42,35,27,25,21,19,18),
   lot2 = c(69,35,26,21,18,16,13,12,12))

mod <- glm(lot1 ~ log(u) + lot2, data = clotting, family = Gamma)

sims_pred <- add_predicted_sims(clotting, mod)
#> Warning: NAs produced

head(sims_pred)
#> # A tibble: 6 × 5
#>       u  lot1  lot2 .sim  lot1_predict
#>   <dbl> <dbl> <dbl> <chr>        <dbl>
#> 1     5   118    69 Y            114. 
#> 2    10    58    35 Y             57.7
#> 3    15    42    26 Y             41.7
#> 4    20    35    21 Y             35.3
#> 5    30    27    18 Y             28.8
#> 6    40    25    16 Y             25.3

# Binomial
# example from Venables and Ripley (2002, pp. 190-2.)
ldose <- rep(0:5, 2)
numdead <- c(1, 4, 9, 13, 18, 20, 0, 2, 6, 10, 12, 16)
sex <- factor(rep(c("M", "F"), c(6, 6)))
SF <- cbind(numdead, numalive = 20-numdead)
budworm.lg <- glm(SF ~ sex*ldose, family = binomial)

dat <- data.frame(sex = factor(c("M", "F", "M", "F")),
                  ldose = c(0,0,5,5))
sims_pred_b <- add_predicted_sims(dat, budworm.lg, weights = 20)

head(sims_pred_b)
#> # A tibble: 6 × 4
#>   sex   ldose .sim  SF_predict
#>   <fct> <dbl> <chr>      <int>
#> 1 M         0 Y              2
#> 2 F         0 Y              2
#> 3 M         5 Y             17
#> 4 F         5 Y             16
#> 5 M         0 Y              0
#> 6 F         0 Y              0