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

```
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')
```

```
## typeprof fitted at 0.
```

*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`

.

```
# 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')
```

```
## rank2 fitted at 0.
```

```
## rank3 fitted at 0.
```

```
## Predicted probabilities calculated outside of the [0, 1] interval.
## These are being removed from the simulated distribution. . .
```

`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')
```

```
## rank2 fitted at 0.
```

```
## rank3 fitted at 0.
```

```
## Predicted probabilities calculated outside of the [0, 1] interval.
## These are being removed from the simulated distribution. . .
```

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')
```

```
## Friendly reminder:
## Unless you are including interactions, to make meaningful plots only x_var and group_var fitted values should vary.
```

```
## rank2 fitted at 0.
```

```
## rank3 fitted at 0.
```

```
## Predicted probabilities calculated outside of the [0, 1] interval.
## These are being removed from the simulated distribution. . .
```

**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.

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.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.