`get_sims()`

is a function to simulate quantities of interest from a
multivariate normal distribution for "new data" from a regression model.

get_sims(model, newdata, nsim, seed)

model | a model object |
---|---|

newdata | A data frame on some quantities of interest to be simulated |

nsim | Number of simulations to be run |

seed | An optional seed to set |

`get_sims()`

returns a data frame (as a `tibble`

) with the quantities
of interest and identifying information about the particular simulation number.

This (should) be a flexible function that takes a `merMod`

object
(estimated from lme4, blme, etc.) or a `lm`

or `glm`

object and
generates some quantities of interest when paired with new data of observations of interest.
Of note: I've really only tested this function with linear models, generalized linear models,
and their mixed model equivalents. For mixed models, this approach does not offer support for
the incorporation of the random effects or the random slopes. It's just for the fixed effects,
which is typically what most people want anyway. Users who want to better incorporate the
random intercepts or slope could find that support in the merTools package.

Steven V. Miller

# Note: these models are dumb, but they illustrate how it works. M1 <- lm(mpg ~ hp, mtcars) # Note: this function requires the DV to appear somewhere, anywhere in the "new data" newdat <- data.frame(mpg = 0, hp = c(mean(mtcars$hp) - sd(mtcars$hp), mean(mtcars$hp), mean(mtcars$hp) + sd(mtcars$hp))) get_sims(M1, newdat, 100, 8675309)#> # A tibble: 300 x 2 #> y sim #> <dbl> <dbl> #> 1 23.8 1 #> 2 19.5 1 #> 3 15.2 1 #> 4 24.8 2 #> 5 20.0 2 #> 6 15.2 2 #> 7 23.9 3 #> 8 19.3 3 #> 9 14.7 3 #> 10 25.1 4 #> # … with 290 more rows# Note: this is likely a dumb model, but illustrates how it works. mtcars$mpgd <- ifelse(mtcars$mpg > 25, 1, 0) M2 <- glm(mpgd ~ hp, mtcars, family=binomial(link="logit")) # Again: this function requires the DV to be somewhere, anywhere in the "new data" newdat$mpgd <- 0 # Note: the simulations are returned on their original "link". Here, that's a "logit" # You can adjust that accordingly. `plogis(y)` will convert those to probabilities. get_sims(M2, newdat, 100, 8675309)#> # A tibble: 300 x 2 #> y sim #> <dbl> <dbl> #> 1 0.312 1 #> 2 -6.85 1 #> 3 -14.0 1 #> 4 0.623 2 #> 5 -2.07 2 #> 6 -4.76 2 #> 7 0.355 3 #> 8 -5.84 3 #> 9 -12.0 3 #> 10 -0.735 4 #> # … with 290 more rows#>M3 <- lmer(mpg ~ hp + (1 | cyl), mtcars) # Random effects are not required here since we're passing over them. get_sims(M3, newdat, 100, 8675309)#> # A tibble: 300 x 2 #> y sim #> <dbl> <dbl> #> 1 21.7 1 #> 2 19.9 1 #> 3 18.0 1 #> 4 24.9 2 #> 5 23.0 2 #> 6 21.1 2 #> 7 22.4 3 #> 8 19.9 3 #> 9 17.4 3 #> 10 19.3 4 #> # … with 290 more rows