sjp.lm {sjPlot}

This document shows examples for using the sjp.lm function of the sjPlot package.


(back to table of content)

Data initialization

Please refer to this document.

# set basic theme options
          axis.title.size = .85, 
          axis.textsize = .85, 
          legend.size = .8, 
          geom.label.size = 3.5)

Fitting a linear model

First, we need to fit a linear model and create a vector with labels for the predictors.

# fit the model
fit <- lm(neg_c_7 ~ c12hour + c160age + e42dep + tot_sc_e + c172code, 
          data = efc)

Plotting coefficients of a linear model

With the sjp.lm function you can plot the beta coefficients with confidence intervals as so called “forest plots”.


As you can see in the above figure, value labels are set automatically. This works because value and variable labels are added to each variable as additional attributes. You can use e.g. set_label and set_labels to manually add label attributes to vectors. If you import data from SPSS etc., the labels are automatically attached as attributes. See this document for more details.

Plotting standardized beta coefficients

You can plot the standardized beta coefficients with type = "std".

sjp.lm(fit, type = "std")

Specific plot customizations

By default, the optimal axis limits are calculated upon the range of lowest and highest confidence interval. If you prefer a symmetrical axis, use the axis.lim parameter.

       breakLabelsAt = 30, 
       axis.lim = c(-1, 1))

Some functions like sjp.lm offer specific customizations of plots for elements that only apply to these functions (e.g. error bars). This is shown in the following figure.

sjp.setTheme(geom.errorbar.size = 1, 
             geom.errorbar.linetype = 2,
             axis.title.size = .85, 
             axis.textsize = .85)
sjp.lm(fit, geom.size = 5)

Marginal effects

Marginal effects plots can be obtained with type = "eff".

# First convert categorical predictors to factor
efc$e42dep <- to_factor(efc$e42dep)
efc$c172code <- to_factor(efc$c172code)
# now fit model again
fit2 <- lm(neg_c_7 ~ c12hour + c160age + e42dep + tot_sc_e + c172code, 
           data = efc)
# and plot marginal effects
sjp.lm(fit2, type = "eff")

As you can see in the above example, multiple plots for type = "eff" are plotted as facet grid resp. as facet wrap. Since this does not allow to set a different x-scale for each plot, x-axis are not properly labelled. However, facet.grid = FALSE produces a single plot for each predictor:

# plot marginal effects, but only for dependency this time
sjp.lm(fit2, type = "eff", facet.grid = FALSE, vars = "e42dep")

To arrange all predictors of multiple in one plot, as grid, use the plot_grid() function on multiple plot objects. plot_grid() requires multiple plots, so you have to set facet.grid = FALSE to get a plot.list-value as return value from the function (see ?sjp.lm on Return Value for more details). This allows you arrange multiple plots as grid in one plot, but with different x-axis-scales.

# get list of all plots
p <- sjp.lm(fit2, type = "eff", facet.grid = FALSE, prnt.plot = FALSE)$plot.list
# plot all marginal effects, as grid, proper x-axes
# also, set margins for this example
plot_grid(p, margin = c(0.2, 0.2, 0.2, 0.2))

Predicting values

With type = "pred", prediction from the fitted model can be obtained. The predicted values of the response, related to specific model predictors. This plot type requires the vars argument to select specific terms that should be used for the x-axis and - optional - as first or second grouping factor. Hence, vars must be a character vector with the names of one to three model predictors.

sjp.lm(fit, type = "pred", vars = "c12hour")

Predicted values can also be grouped by another predictor. In this case, you have to provide two predictor names for the vars argument.

sjp.lm(fit, type = "pred", vars = c("c12hour", "e42dep"))

With the scatter.plot-argument you can prevent the data points from being plotted.

sjp.lm(fit, type = "pred", vars = c("c12hour", "e42dep"), scatter.plot = FALSE)

Use facet.grid argument to avoid faceted plots.

sjp.lm(fit, type = "pred", 
       vars = c("c12hour", "e42dep"), facet.grid = FALSE)

With three terms specified in vars, you can add another grouping level. Predictions with three variables in vars are automatically plotted in grid layout (faceted).

sjp.lm(fit, type = "pred", vars = c("c12hour", "e42dep", "c172code"))

Investigate single predictors

If you are not sure if a predictor might be improper for linear models (because of skewness or whatever), you can investigate each single predictor of the model with the parameter type = "slope". This call requires at least the fitted model. The result are several figures (one for each predictor) of linear models where each predictor alone is fitted agains the response.

This function plots the linear fitted line, a loess fitted line and the standard error range by default. The better the linear and loess fitted line “match” (overlap), the better is the linear relationship between predictor and response.

sjp.lm(fit, type = "slope", show.loess = TRUE)

Testing model assumptions

You can assess the model assumptions of your fitted model with the parameter type = "ma". This call produces some information and several plots:

  1. An outlier test. This test removes outliers and returns an updated fitted model.
  2. The variance inflation factor for each predictor (test for multicollinearity)
  3. Non-normality of residuals and outliers
  4. Non-normality (distribution) of residuals with overlayed normal curve (in red) for comparision
  5. Test for homoscedasticity. A blue loess fitted line is also printed, which should be linear and match with the black horizontal line in best case.
sjp.lm(fit, type = "ma")
## Removed 7 cases during 1 step(s).
## R^2 / adj. R^2 of original model: 0.159176 / 0.154093
## R^2 / adj. R^2 of updated model:  0.179938 / 0.171908
## AIC of original model: 4497.751285 
## AIC of updated model:  4358.458892