sjp.lm: Plot estimates, predictions or effects of linear models

Description

Depending on the type, this function plots coefficients (estimates) of linear regressions (including panel models fitted with the plm-function from the plm-package and generalized least squares models fitted with the gls-function from the nlme-package) with confidence intervals as dot plot (forest plot), model assumptions for linear models or slopes and scatter plots for each single coefficient. See type for details.

Usage

sjp.lm(fit, type = "lm", vars = NULL, group.estimates = NULL,
  remove.estimates = NULL, sort.est = TRUE, poly.term = NULL,
  title = NULL, legend.title = NULL, axis.labels = NULL,
  axis.title = NULL, resp.label = NULL, geom.size = NULL,
  geom.colors = "Set1", wrap.title = 50, wrap.labels = 25,
  axis.lim = NULL, grid.breaks = NULL, show.values = TRUE,
  show.p = TRUE, show.ci = TRUE, show.legend = FALSE,
  show.loess = FALSE, show.loess.ci = FALSE, show.summary = FALSE,
  show.scatter = TRUE, point.alpha = 0.2, point.color = NULL,
  jitter.ci = FALSE, digits = 2, vline.type = 2, vline.color = "grey70",
  coord.flip = TRUE, y.offset = 0.15, facet.grid = TRUE,
  complete.dgns = FALSE, prnt.plot = TRUE, ...)

Arguments

fit

Fitted linear regression model (of class lm, gls or plm).

type

Type of plot. Use one of following:

"lm": (default) for forest-plot of estimates. If the fitted model only contains one predictor, slope-line is plotted.
"pred": to plot predicted values (marginal effects) for specific model terms. See 'Details'.
"eff": to plot marginal effects of all terms in fit. Note that interaction terms are excluded from this plot.
"std": for forest-plot of standardized beta values.
"std2": for forest-plot of standardized beta values, however, standardization is done by dividing by two sd (see 'Details').
"resid": to plot regression lines for each single predictor of the fitted model, against the residuals (linear relationship between each model term and residuals). May be used for model diagnostics.
"ma": to check model assumptions.
"vif": to plot Variance Inflation Factors.

vars

Numeric vector with column indices of selected variables or a character vector with variable names of selected variables from the fitted model, which should be used to plot - depending on type - estimates, fixed effects slopes or predicted values (mean, probabilities, incidents rates, ...). See 'Examples'.

group.estimates

Numeric or character vector, indicating a group identifier for each estimate. Dots and confidence intervals of estimates are coloured according to their group association. See 'Examples'.

remove.estimates

Character vector with coefficient names that indicate which estimates should be removed from the plot. remove.estimates = "est_name" would remove the estimate est_name. Default is NULL, i.e. all estimates are printed.

sort.est

Logical, determines whether estimates should be sorted according to their values. If group.estimates is not NULL, estimates are sorted according to their group assignment.

poly.term

name of a polynomial term in fit as string. Needs to be specified, if type = "poly", in order to plot marginal effects for polynomial terms. See 'Examples'.

title

character vector, used as plot title. Depending on plot type and function, will be set automatically. If title = "", no title is printed. For effect-plots, may also be a character vector of length > 1, to define titles for each sub-plot or facet.

legend.title

Character vector, used as title for the plot legend. Note that only some plot types have legends (e.g. type = "pred" or when grouping estimates with group.estimates).

axis.labels

character vector with labels used as axis labels. Optional argument, since in most cases, axis labels are set automatically.

axis.title

Character vector of length one or two (depending on the plot function and type), used as title(s) for the x and y axis. If not specified, a default labelling is chosen. To set multiple axis titles (e.g. with type = "eff" for many predictors), axis.title must be a character vector of same length of plots that are printed. In this case, each plot gets an own axis title (applying, for instance, to the y-axis for type = "eff"). Note: Some plot types do not support this argument. In such cases, use the return value and add axis titles manually with labs, e.g.: $plot.list[[1]] + labs(x = ...)

resp.label

Name of dependent variable, as string. Only used if fitted model has only one predictor and type = "lm".

geom.size

size resp. width of the geoms (bar width, line thickness or point size, depending on plot type and function). Note that bar and bin widths mostly need smaller values than dot sizes.

geom.colors

User defined color palette for geoms. If group.estimates is not specified, must either be vector with two color values or a specific color palette code (see 'Details' in sjp.grpfrq). Else, if group.estimates is specified, geom.colors must be a vector of same length as groups. See 'Examples'.

wrap.title

numeric, determines how many chars of the plot title are displayed in one line and when a line break is inserted.

wrap.labels

numeric, determines how many chars of the value, variable or axis labels are displayed in one line and when a line break is inserted.

axis.lim

Numeric vector of length 2, defining the range of the plot axis. Depending on plot type, may effect either x- or y-axis, or both. For multiple plot outputs (e.g., from type = "eff" or type = "slope" in sjp.glm), axis.lim may also be a list of vectors of length 2, defining axis limits for each plot (only if non-faceted).

grid.breaks

numeric; sets the distance between breaks for the axis, i.e. at every grid.breaks'th position a major grid is being printed.

show.values

Logical, whether values should be plotted or not.

show.p

Logical, adds significance levels to values, or value and variable labels.

show.ci

Logical, if TRUE, depending on type, a confidence interval or region is added to the plot. For frequency plots, the confidence interval for the relative frequencies are shown.

show.legend

logical, if TRUE, and depending on plot type and function, a legend is added to the plot.

show.loess

logical, if TRUE, and depending on type, an additional loess-smoothed line is plotted.

show.loess.ci

logical, if TRUE, a confidence region for the loess-smoothed line will be plotted. Default is FALSE. Only applies, if show.loess = TRUE (and for sjp.lmer, only applies if type = "fe.slope" or type = "fe.resid").

show.summary

Logical, if TRUE, a summary with model statistics is added to the plot.

show.scatter

Logical, if TRUE (default), adds a scatter plot of data points to the plot. Only applies for slope-type or predictions plots. For most plot types, dots are jittered to avoid overplotting, hence the points don't reflect exact values in the data.

point.alpha

Alpha value of point-geoms in the scatter plots. Only applies, if show.scatter = TRUE.

point.color

Color of of point-geoms in the scatter plots. Only applies, if show.scatter = TRUE.

jitter.ci

Logical, if TRUE and show.ci = TRUE and confidence bands are displayed as error bars, adds jittering to lines and error bars to avoid overlapping.

digits

Numeric, amount of digits after decimal point when rounding estimates and values.

vline.type

Linetype of the vertical "zero point" line. Default is 2 (dashed line).

vline.color

Color of the vertical "zero point" line. Default value is "grey70".

coord.flip

logical, if TRUE, the x and y axis are swapped.

y.offset

numeric, offset for text labels when their alignment is adjusted to the top/bottom of the geom (see hjust and vjust).

facet.grid

TRUE to arrange the lay out of of multiple plots in a grid of an integrated single plot. This argument calls facet_wrap or facet_grid to arrange plots. Use plot_grid to plot multiple plot-objects as an arranged grid with grid.arrange.

complete.dgns

Logical, if TRUE, additional tests are performed. Default is FALSE Only applies if type = "ma".

prnt.plot

logical, if TRUE (default), plots the results as graph. Use FALSE if you don't want to plot any graphs. In either case, the ggplot-object will be returned as value.

...

Other arguments passed down to further functions. Currently, following arguments are supported:

?effects::effect: Any arguments accepted by the effect resp. allEffects function, for type = "eff".
width: The width-argument for error bars.
alpha: The alpha-argument for confidence bands.
level: The level-argument confidence bands.

Value

Depending on the type, in most cases (insisibily) returns the ggplot-object with the complete plot (plot) as well as the data frame that was used for setting up the ggplot-object (df). For type = "ma", an updated model with removed outliers is returned.

Details

type = "lm": if fitted model only has one predictor, no forest plot is shown. Instead, a regression line with confidence interval (in blue) is plotted by default, and a loess-smoothed line without confidence interval (in red) can be added if argument show.loess = TRUE.
type = "std2": plots standardized beta values, however, standardization follows Gelman's (2008) suggestion, rescaling the estimates by dividing them by two standard deviations instead of just one. Resulting coefficients are then directly comparable for untransformed binary predictors. This standardization uses the standardize-function from the arm-package.
type = "slope": regression lines (slopes) with confidence intervals for each single predictor of the fitted model are plotted, i.e. all predictors of the fitted model are extracted and for each of them, the linear relationship is plotted against the response variable. Other predictors are omitted, so this plot type is intended to check the linear relationship between a predictor and the response.
type = "resid": is similar to the type = "slope" option, however, each predictor is plotted against the residuals (instead of response).
type = "pred": plots predicted values of the response, related to specific model predictors. This plot type calls predict(fit, newdata = model.frame, type = "response") and requires the vars argument to select specific terms that should be used for the x-axis and - optional - as grouping factor. Hence, vars must be a character vector with the names of one or two model predictors. See 'Examples'.
type = "eff": computes the marginal effects for all predictors, using the allEffects function. I.e. for each predictor, the predicted values towards the response are plotted, with all remaining co-variates set to the mean. Due to possible different scales of predictors, a faceted plot is printed (instead of plotting all lines in one plot). You can pass further arguments down to allEffects for flexible function call via the ...-argument.
type = "poly": plots the marginal effects of polynomial terms in fit, using the effect function, but only for a selected polynomial term, which is specified with poly.term. This function helps undertanding the effect of polynomial terms by plotting the curvilinear relationships of response and quadratic, cubic etc. terms. This function accepts following argument.
type = "ma": checks model assumptions. Please note that only three arguments are relevant: fit and complete.dgns. All other arguments are ignored.
type = "vif": Variance Inflation Factors (check for multicollinearity) are plotted. As a rule of thumb, values below 5 are considered as good and indicate no multicollinearity, values between 5 and 10 may be tolerable. Values greater than 10 are not acceptable and indicate multicollinearity between model's predictors.

References

Gelman A (2008) "Scaling regression inputs by dividing by two standard deviations." Statistics in Medicine 27: 2865<U+2013>2873. http://www.stat.columbia.edu/~gelman/research/published/standardizing7.pdf

Hyndman RJ, Athanasopoulos G (2013) "Forecasting: principles and practice." OTexts; accessed from https://www.otexts.org/fpp/5/4.

Examples

Run this code

# NOT RUN {
# --------------------------------------------------
# plotting estimates of linear models as forest plot
# --------------------------------------------------
# fit linear model
fit <- lm(airquality$Ozone ~ airquality$Wind + airquality$Temp + airquality$Solar.R)

# plot estimates with CI
sjp.lm(fit, grid.breaks = 2)

# plot estimates with CI
# and with narrower tick marks
# (because "grid.breaks" was not specified)
sjp.lm(fit)

# ---------------------------------------------------
# plotting regression line of linear model (done
# automatically if fitted model has only 1 predictor)
# ---------------------------------------------------
library(sjmisc)
data(efc)
# fit model
fit <- lm(neg_c_7 ~ quol_5, data=efc)
# plot regression line with label strings
sjp.lm(fit, resp.label = "Burden of care",
       axis.labels = "Quality of life", show.loess = TRUE)

# --------------------------------------------------
# plotting regression lines of each single predictor
# of a fitted model
# --------------------------------------------------
library(sjmisc)
data(efc)
# fit model
fit <- lm(tot_sc_e ~ c12hour + e17age + e42dep, data=efc)

# reression line and scatter plot
sjp.lm(fit, type = "slope")

# reression line w/o scatter plot
sjp.lm(fit, type = "slope", show.scatter = FALSE)

# --------------------------
# plotting model assumptions
# --------------------------
sjp.lm(fit, type = "ma")

# }
# NOT RUN {
# --------------------------
# grouping estimates
# --------------------------
library(sjmisc)
data(efc)
fit <- lm(barthtot ~ c160age + e17age + c12hour + e16sex + c161sex + c172code,
          data = efc)

# order estimates according to coefficient's order
sjp.lm(fit, group.estimates = c(1, 1, 2, 3, 3, 4),
       geom.colors = c("green", "red", "blue", "grey"), sort.est = FALSE)

fit <- lm(barthtot ~ c160age + c12hour + e17age+ c161sex + c172code + e16sex,
          data = efc)

# force order of estimates according to group assignment
sjp.lm(fit, group.estimates = c(1, 2, 1, 3, 4, 3),
       geom.colors = c("green", "red", "blue", "grey"), sort.est = TRUE)

# --------------------------
# predicted values for response
# --------------------------
library(sjmisc)
data(efc)
efc$education <- to_label(to_factor(efc$c172code))
efc$gender <- to_label(to_factor(efc$c161sex))
fit <- lm(barthtot ~ c160age + c12hour + e17age + gender + education,
          data = efc)

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

# with loess
sjp.lm(fit, type = "pred", vars = "e17age", show.loess = TRUE)

# grouped
sjp.lm(fit, type = "pred", vars = c("c12hour", "education"))

# grouped, non-facet
sjp.lm(fit, type = "pred", vars = c("c12hour", "education"),
       facet.grid = FALSE)

# two groupings
sjp.lm(fit, type = "pred", vars = c("c12hour", "gender", "education"))

# --------------------------
# plotting polynomial terms
# --------------------------
library(sjmisc)
data(efc)
# fit sample model
fit <- lm(tot_sc_e ~ c12hour + e17age + e42dep, data = efc)
# "e17age" does not seem to be linear correlated to response
# try to find appropiate polynomial. Grey line (loess smoothed)
# indicates best fit. Looks like x^3 has a good fit.
# (not checked for significance yet).
sjp.poly(fit, "e17age", 2:4, show.scatter = FALSE)
# fit new model
fit <- lm(tot_sc_e ~ c12hour + e42dep +
          e17age + I(e17age^2) + I(e17age^3),
          data = efc)
# plot marginal effects of polynomial term
sjp.lm(fit, type = "poly", poly.term = "e17age")

library(splines)
# fit new model with "splines"-package, "bs"
fit <- lm(tot_sc_e ~ c12hour + e42dep + bs(e17age, 3), data = efc)
# plot marginal effects of polynomial term, same call as above
sjp.lm(fit, type = "poly", poly.term = "e17age")
# }
# NOT RUN {
# }

Run the code above in your browser using DataLab