In christophergandrud/simGLM: Simulate and Plot Quantities of Interest from Generalised Linear Models

simGLM

Simulate and plot quantities of interest from generalised linear models using King, Tomz, and Wittenburg (2000). Currently only supports normal linear and logistic regression models.

Christopher Gandrud

Example: normal linear regression

library(car) # contains data
library(simGLM)
library(ggplot2) # only needed for adding additional arguments outside of sim_glm

# Estimate model
m1 <- lm(prestige ~ education + type, data = Prestige)

# Create fitted values
fitted_prestige <- expand.grid(education = 6:16, typewc = 1)

# Simulate and plot
sim_glm(obj = m1, newdata = fitted_prestige, x_coef = 'education') +
        ylab('Predicted Job Prestige\n') + xlab('\nYears of Education')

Note: when you create the data frame for the fitted values for your simulations, each column needs to be assigned a name that exactly matches one of the coefficient names in the model summary.

Example: logistic regression

# Download data
URL <- 'http://www.ats.ucla.edu/stat/data/binary.csv'
Admission <- read.csv(URL)
Admission$rank <- as.factor(Admission$rank)

# Estimate model
m2 <- glm(admit ~ gre + gpa + rank, data = Admission, family = 'binomial')

# Create fitted values
fitted_admit <- expand.grid(gre = seq(220, 800, by = 10), gpa = c(2, 4), 
                            rank4 = 1)

# Simulate and plot
sim_glm(obj = m2, newdata = fitted_admit, model = 'logit', x_coef = 'gre', 
        group_coef = 'gpa')

Examples: bayesglm

sim_glm also works with estimates made using the bayesglm funciton in the arm package. This function uses minimal prior information suggested by Gelman et al. (2008) to avoid well known problems of unrealistic logistic regression coefficient sizes and, in the extreme case, complete separation

library(arm)

# Estimate model
m3 <- bayesglm(admit ~ gre + gpa + rank, data = Admission, 
               family = binomial(link = 'logit'))

# Simulate and plot
sim_glm(obj = m3, newdata = fitted_admit, model = 'logit', x_coef = 'gre', 
        group_coef = 'gpa')

Example: Interactions and polynomials

If you model multi-term effects with interactions and polynomials then you need to specify fitted values for the interaction/polynomial terms, not just the base term. For example:

# Estimate model
m4 <- glm(admit ~ gre * gpa + rank, data = Admission, family = 'binomial')

fitted_admit$`gre:gpa` <- fitted_admit$gre * fitted_admit$gpa

# Simulate and plot
sim_glm(obj = m4, newdata = fitted_admit, model = 'logit', x_coef = 'gre', 
        group_coef = 'gpa')

Custom Quantities of Interest

Experimental: in order to allow the user to specify other quantites of interest for other GLM type models, the model argument of sim_glm allows you to specify a function with which to calculate a custom quantity of interest. The function must take as its input a vector of values from your simulated point estimates and fitted values (e.g. alpha + beta1 * x1 + beta2 * x2) and return a numeric vector with your custom quantity of interest.

Install

To install the development version of simGLM use:

devtools::install_github('christophergandrud/simGLM')

christophergandrud/simGLM documentation built on May 13, 2019, 7:03 p.m.