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_fitted_sims(newdata, mod, n_sims = 1000, seed = NULL, weights = 1, ...)

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

...

Unused dots for compatibility with generic functions.

Value

A tibble::tibble() with information about simulate values.

See also

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_fit <- add_fitted_sims(clotting, mod)

head(sims_fit)
#> # A tibble: 6 × 5
#>       u  lot1  lot2 .sim  lot1_fit
#>   <dbl> <dbl> <dbl> <chr>    <dbl>
#> 1     5   118    69 1         119.
#> 2     5   118    69 2         116.
#> 3     5   118    69 3         114.
#> 4     5   118    69 4         123.
#> 5     5   118    69 5         115.
#> 6     5   118    69 6         120.

# 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_fit_b <- add_fitted_sims(dat, budworm.lg)

head(sims_fit_b)
#> # A tibble: 6 × 4
#>   sex   ldose .sim  SF_fit
#>   <fct> <dbl> <chr>  <dbl>
#> 1 M         0 1     0.126 
#> 2 M         0 2     0.117 
#> 3 M         0 3     0.0817
#> 4 M         0 4     0.0588
#> 5 M         0 5     0.0424
#> 6 M         0 6     0.0253