-
Notifications
You must be signed in to change notification settings - Fork 92
/
plot_model.R
870 lines (808 loc) · 39.5 KB
/
plot_model.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
#' @title Plot regression models
#' @name plot_model
#'
#' @description
#' \code{plot_model()} creates plots from regression models, either
#' estimates (as so-called forest or dot whisker plots) or marginal effects.
#'
#' @param model A regression model object. Depending on the \code{type}, many
#' kinds of models are supported, e.g. from packages like \pkg{stats},
#' \pkg{lme4}, \pkg{nlme}, \pkg{rstanarm}, \pkg{survey}, \pkg{glmmTMB},
#' \pkg{MASS}, \pkg{brms} etc.
#' @param type Type of plot. There are three groups of plot-types: \cr \cr
#' \emph{Coefficients} (\href{https://strengejacke.github.io/sjPlot/articles/plot_model_estimates.html}{related vignette})
#' \describe{
#' \item{\code{type = "est"}}{Forest-plot of estimates. If the fitted model
#' only contains one predictor, slope-line is plotted.}
#' \item{\code{type = "re"}}{For mixed effects models, plots the random
#' effects.}
#' \item{\code{type = "std"}}{Forest-plot of standardized coefficients.}
#' \item{\code{type = "std2"}}{Forest-plot of standardized coefficients,
#' however, standardization is done by dividing by two SD (see 'Details').}
#' }
#' \emph{Marginal Effects} (\href{https://strengejacke.github.io/sjPlot/articles/plot_marginal_effects.html}{related vignette})
#' \describe{
#' \item{\code{type = "pred"}}{Predicted values (marginal effects) for
#' specific model terms. See \code{\link[ggeffects]{ggpredict}} for details.}
#' \item{\code{type = "eff"}}{Similar to \code{type = "pred"}, however,
#' discrete predictors are held constant at their proportions (not reference
#' level). See \code{\link[ggeffects]{ggeffect}} for details.}
#' \item{\code{type = "emm"}}{Similar to \code{type = "eff"}, see
#' \code{\link[ggeffects]{ggemmeans}} for details.}
#' \item{\code{type = "int"}}{Marginal effects of interaction terms in
#' \code{model}.}
#' }
#' \emph{Model diagnostics}
#' \describe{
#' \item{\code{type = "slope"}}{Slope of coefficients for each single
#' predictor, against the response (linear relationship between each model
#' term and response). See 'Details'.}
#' \item{\code{type = "resid"}}{Slope of coefficients for each single
#' predictor, against the residuals (linear relationship between each model
#' term and residuals). See 'Details'.}
#' \item{\code{type = "diag"}}{Check model assumptions. See 'Details'.}
#' }
#' \strong{Note:} For mixed models, the diagnostic plots like linear relationship
#' or check for Homoscedasticity, do \strong{not} take the uncertainty of
#' random effects into account, but is only based on the fixed effects part
#' of the model.
#' @param transform A character vector, naming a function that will be applied
#' on estimates and confidence intervals. By default, \code{transform} will
#' automatically use \code{"exp"} as transformation for applicable classes of
#' \code{model} (e.g. logistic or poisson regression). Estimates of linear
#' models remain untransformed. Use \code{NULL} if you want the raw,
#' non-transformed estimates.
#' @param terms Character vector with the names of those terms from \code{model}
#' that should be plotted. This argument depends on the plot-type:
#' \describe{
#' \item{\emph{Coefficients}}{Select terms that should be plotted. All other
#' term are removed from the output. Note that the term names must match
#' the names of the model's coefficients. For factors, this means that
#' the variable name is suffixed with the related factor level, and each
#' category counts as one term. E.g. \code{rm.terms = "t_name [2,3]"}
#' would remove the terms \code{"t_name2"} and \code{"t_name3"} (assuming
#' that the variable \code{t_name} is categorical and has at least
#' the factor levels \code{2} and \code{3}). Another example for the
#' \emph{iris}-dataset: \code{terms = "Species"} would not work, instead
#' you would write \code{terms = "Species [versicolor,virginica]"} to
#' remove these two levels, or \code{terms = "Speciesversicolor"} if you
#' just want to remove the level \emph{versicolor} from the plot.}
#' \item{\emph{Marginal Effects}}{Here \code{terms} indicates for which
#' terms marginal effects should be displayed. At least one term is
#' required to calculate effects, maximum length is three terms, where
#' the second and third term indicate the groups, i.e. predictions of
#' first term are grouped by the levels of the second (and third) term.
#' \code{terms} may also indicate higher order terms (e.g. interaction
#' terms). Indicating levels in square brackets allows for selecting only
#' specific groups. Term name and levels in brackets must be separated by
#' a whitespace character, e.g. \code{terms = c("age", "education [1,3]")}.
#' It is also possible to specify a range of numeric values for the
#' predictions with a colon, for instance \code{terms = c("education [1,3]",
#' "age [30:50]")}. Furthermore, it is possible to specify a function name.
#' Values for predictions will then be transformed, e.g.
#' \code{terms = "income [exp]"}. This is useful when model predictors were
#' transformed for fitting the model and should be back-transformed to the
#' original scale for predictions. Finally, numeric vectors for which no
#' specific values are given, a "pretty range" is calculated, to avoid
#' memory allocation problems for vectors with many unique values. If a
#' numeric vector is specified as second or third term (i.e. if this vector
#' represents a grouping structure), representative values (see
#' \code{\link[ggeffects]{values_at}}) are chosen. If all values for a
#' numeric vector should be used to compute predictions, you may use
#' e.g. terms = "age [all]". For more details, see
#' \code{\link[ggeffects]{ggpredict}}.}
#' }
#' @param sort.est Determines in which way estimates are sorted in the plot:
#' \itemize{
#' \item If \code{NULL} (default), no sorting is done and estimates are sorted in the same order as they appear in the model formula.
#' \item If \code{TRUE}, estimates are sorted in descending order, with highest estimate at the top.
#' \item If \code{sort.est = "sort.all"}, estimates are re-sorted for each coefficient (only applies if \code{type = "re"} and \code{grid = FALSE}), i.e. the estimates of the random effects for each predictor are sorted and plotted to an own plot.
#' \item If \code{type = "re"}, specify a predictor's / coefficient's name to sort estimates according to this random effect.
#' }
#' @param rm.terms Character vector with names that indicate which terms should
#' be removed from the plot. Counterpart to \code{terms}. \code{rm.terms =
#' "t_name"} would remove the term \emph{t_name}. Default is \code{NULL}, i.e.
#' all terms are used. For factors, levels that should be removed from the plot
#' need to be explicitely indicated in square brackets, and match the model's
#' coefficient names, e.g. \code{rm.terms = "t_name [2,3]"} would remove the terms
#' \code{"t_name2"} and \code{"t_name3"} (assuming that the variable \code{t_name}
#' was categorical and has at least the factor levels \code{2} and \code{3}).
#' Another example for the \emph{iris} dataset would be
#' \code{rm.terms = "Species [versicolor,virginica]"}. Note that the
#' \code{rm.terms}-argument does not apply to \emph{Marginal Effects} plots.
#' @param group.terms Numeric vector with group indices, to group coefficients.
#' Each group of coefficients gets its own color (see 'Examples').
#' @param order.terms Numeric vector, indicating in which order the coefficients
#' should be plotted. See examples in
#' \href{https://strengejacke.github.io/sjPlot/articles/plot_model_estimates.html}{this package-vignette}.
#' @param pred.type Character, only applies for \emph{Marginal Effects} plots
#' with mixed effects models. Indicates whether predicted values should be
#' conditioned on random effects (\code{pred.type = "re"}) or fixed effects
#' only (\code{pred.type = "fe"}, the default). For details, see documentation
#' of the \code{type}-argument in \code{\link[ggeffects]{ggpredict}}.
#' @param mdrt.values Indicates which values of the moderator variable should be
#' used when plotting interaction terms (i.e. \code{type = "int"}). \describe{
#' \item{\code{"minmax"}}{(default) minimum and maximum values (lower and
#' upper bounds) of the moderator are used to plot the interaction between
#' independent variable and moderator(s).} \item{\code{"meansd"}}{uses the
#' mean value of the moderator as well as one standard deviation below and
#' above mean value to plot the effect of the moderator on the independent
#' variable (following the convention suggested by Cohen and Cohen and
#' popularized by Aiken and West (1991), i.e. using the mean, the value one
#' standard deviation above, and the value one standard deviation below the
#' mean as values of the moderator, see
#' \href{https://www.theanalysisfactor.com/3-tips-interpreting-moderation/}{Grace-Martin
#' K: 3 Tips to Make Interpreting Moderation Effects Easier}).}
#' \item{\code{"zeromax"}}{is similar to the \code{"minmax"} option, however,
#' \code{0} is always used as minimum value for the moderator. This may be
#' useful for predictors that don't have an empirical zero-value, but absence
#' of moderation should be simulated by using 0 as minimum.}
#' \item{\code{"quart"}}{calculates and uses the quartiles (lower, median and
#' upper) of the moderator value.} \item{\code{"all"}}{uses all values of the
#' moderator variable.} }
#' @param ri.nr Numeric vector. If \code{type = "re"} and fitted model has more
#' than one random intercept, \code{ri.nr} indicates which random effects of
#' which random intercept (or: which list elements of
#' \code{\link[lme4]{ranef}}) will be plotted. Default is \code{NULL}, so all
#' random effects will be plotted.
#' @param title Character vector, used as plot title. By default,
#' \code{\link[sjlabelled]{response_labels}} is called to retrieve the label of
#' the dependent variable, which will be used as title. Use \code{title = ""}
#' to remove title.
#' @param 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. \strong{Note:} Some plot types
#' may not support this argument sufficiently. In such cases, use the returned
#' ggplot-object and add axis titles manually with
#' \code{\link[ggplot2]{labs}}. Use \code{axis.title = ""} to remove axis
#' titles.
#' @param axis.labels Character vector with labels for the model terms, used as
#' axis labels. By default, \code{\link[sjlabelled]{term_labels}} is
#' called to retrieve the labels of the coefficients, which will be used as
#' axis labels. Use \code{axis.labels = ""} or \code{auto.label = FALSE} to
#' use the variable names as labels instead. If \code{axis.labels} is a named
#' vector, axis labels (by default, the names of the model's coefficients)
#' will be matched with the names of \code{axis.label}. This ensures that
#' labels always match the related axis value, no matter in which way
#' axis labels are sorted.
#' @param 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. For
#' \emph{Marginal Effects} plots, \code{axis.lim} may also be a list of two
#' vectors of length 2, defining axis limits for both the x and y axis.
#' @param legend.title Character vector, used as legend title for plots that
#' have a legend.
#' @param grid.breaks Numeric value or vector; if \code{grid.breaks} is a
#' single value, sets the distance between breaks for the axis at every
#' \code{grid.breaks}'th position, where a major grid line is plotted. If
#' \code{grid.breaks} is a vector, values will be used to define the
#' axis positions of the major grid lines.
#' @param ci.lvl Numeric, the level of the confidence intervals (error bars).
#' Use \code{ci.lvl = NA} to remove error bars. For \code{stanreg}-models,
#' \code{ci.lvl} defines the (outer) probability for the \emph{credible interval}
#' that is plotted (see \code{\link[bayestestR]{ci}}). By
#' default, \code{stanreg}-models are printed with two intervals: the "inner"
#' interval, which defaults to the 50\%-CI; and the "outer" interval, which
#' defaults to the 89\%-CI. \code{ci.lvl} affects only the outer interval in
#' such cases. See \code{prob.inner} and \code{prob.outer} under the
#' \code{...}-argument for more details.
#' @param se Logical, if \code{TRUE}, the standard errors are
#' also printed. If robust standard errors are required, use arguments
#' \code{vcov.fun}, \code{vcov.type} and \code{vcov.args} (see
#' \code{\link[parameters]{standard_error_robust}} and
#' \href{https://easystats.github.io/parameters/articles/model_parameters_robust.html}{this vignette}
#' for details), or use argument \code{robust} as shortcut. \code{se} overrides
#' \code{ci.lvl}: if not \code{NULL}, arguments \code{ci.lvl} and \code{transform}
#' will be ignored. Currently, \code{se} only applies to \emph{Coefficients} plots.
#' @param show.intercept Logical, if \code{TRUE}, the intercept of the fitted
#' model is also plotted. Default is \code{FALSE}. If \code{transform =
#' "exp"}, please note that due to exponential transformation of estimates,
#' the intercept in some cases is non-finite and the plot can not be created.
#' @param show.values Logical, whether values should be plotted or not.
#' @param show.p Logical, adds asterisks that indicate the significance level of
#' estimates to the value labels.
#' @param show.data Logical, for \emph{Marginal Effects} plots, also plots the
#' raw data points.
#' @param show.legend For \emph{Marginal Effects} plots, shows or hides the
#' legend.
#' @param show.zeroinf Logical, if \code{TRUE}, shows the zero-inflation part of
#' hurdle- or zero-inflated models.
#' @param robust Logical, shortcut for arguments \code{vcov.fun} and \code{vcov.type}.
#' If \code{TRUE}, uses \code{vcov.fun = "vcovHC"} and \code{vcov.type = "HC3"} as
#' default, that is, \code{\link[sandwich]{vcovHC}} with default-type is called
#' (see \code{\link[parameters]{standard_error_robust}} and
#' \href{https://easystats.github.io/parameters/articles/model_parameters_robust.html}{this vignette}
#' for further details).
#' @param vcov.fun Character vector, indicating the name of the \code{vcov*()}-function
#' from the \pkg{sandwich} or \pkg{clubSandwich} package, e.g. \code{vcov.fun = "vcovCL"},
#' if robust standard errors are required.
#' @param vcov.type Character vector, specifying the estimation type for the
#' robust covariance matrix estimation (see \code{\link[sandwich:vcovHC]{vcovHC()}}
#' or \code{clubSandwich::vcovCR()} for details).
#' @param vcov.args List of named vectors, used as additional arguments that
#' are passed down to \code{vcov.fun}.
#' @param value.offset Numeric, offset for text labels to adjust their position
#' relative to the dots or lines.
#' @param dot.size Numeric, size of the dots that indicate the point estimates.
#' @param line.size Numeric, size of the lines that indicate the error bars.
#' @param colors May be a character vector of color values in hex-format, valid
#' color value names (see \code{demo("colors")}) or a name of a pre-defined
#' color palette. Following options are valid for the \code{colors} argument:
#' \itemize{
#' \item If not specified, a default color brewer palette will be used, which is suitable for the plot style.
#' \item If \code{"gs"}, a greyscale will be used.
#' \item If \code{"bw"}, and plot-type is a line-plot, the plot is black/white and uses different line types to distinguish groups (see \href{https://strengejacke.github.io/sjPlot/articles/blackwhitefigures.html}{this package-vignette}).
#' \item If \code{colors} is any valid color brewer palette name, the related palette will be used. Use \code{RColorBrewer::display.brewer.all()} to view all available palette names.
#' \item There are some pre-defined color palettes in this package, see \code{\link{sjPlot-themes}} for details.
#' \item Else specify own color values or names as vector (e.g. \code{colors = "#00ff00"} or \code{colors = c("firebrick", "blue")}).
#' }
#' @param grid Logical, if \code{TRUE}, multiple plots are plotted as grid
#' layout.
#' @param p.threshold Numeric vector of length 3, indicating the treshold for
#' annotating p-values with asterisks. Only applies if
#' \code{p.style = "asterisk"}.
#' @param wrap.title Numeric, determines how many chars of the plot title are
#' displayed in one line and when a line break is inserted.
#' @param 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.
#' @param case Desired target case. Labels will automatically converted into the
#' specified character case. See \code{snakecase::to_any_case()} for more
#' details on this argument. By default, if \code{case} is not specified,
#' it will be set to \code{"parsed"}, unless \code{prefix.labels} is not
#' \code{"none"}. If \code{prefix.labels} is either \code{"label"} (or
#' \code{"l"}) or \code{"varname"} (or \code{"v"}) and \code{case} is not
#' specified, it will be set to \code{NULL} - this is a more convenient
#' default when prefixing labels.
#' @param auto.label Logical, if \code{TRUE} (the default),
#' and \href{https://strengejacke.github.io/sjlabelled/articles/intro_sjlabelled.html}{data is labelled},
#' \code{\link[sjlabelled]{term_labels}} is called to retrieve the labels
#' of the coefficients, which will be used as predictor labels. If data is
#' not labelled, \href{https://easystats.github.io/parameters/reference/format_parameters.html}{format_parameters()}
#' is used to create pretty labels. If \code{auto.label = FALSE},
#' original variable names and value labels (factor levels) are used.
#' @param prefix.labels Indicates whether the value labels of categorical variables
#' should be prefixed, e.g. with the variable name or variable label. See
#' argument \code{prefix} in \code{\link[sjlabelled]{term_labels}} for
#' details.
#' @param jitter Numeric, between 0 and 1. If \code{show.data = TRUE}, you can
#' add a small amount of random variation to the location of each data point.
#' \code{jitter} then indicates the width, i.e. how much of a bin's width
#' will be occupied by the jittered values.
#' @param digits Numeric, amount of digits after decimal point when rounding
#' estimates or values.
#' @param p.adjust Character vector, if not \code{NULL}, indicates the method
#' to adjust p-values. See \code{\link[stats]{p.adjust}} for details.
#' @param value.size Numeric, indicates the size of value labels. Can be used
#' for all plot types where the argument \code{show.values} is applicable,
#' e.g. \code{value.size = 4}.
#' @param vline.color Color of the vertical "zero effect" line. Default color is
#' inherited from the current theme.
#' @param bpe For \strong{Stan}-models (fitted with the \pkg{rstanarm}- or
#' \pkg{brms}-package), the Bayesian point estimate is, by default, the median
#' of the posterior distribution. Use \code{bpe} to define other functions to
#' calculate the Bayesian point estimate. \code{bpe} needs to be a character
#' naming the specific function, which is passed to the \code{fun}-argument in
#' \code{\link[sjmisc]{typical_value}}. So, \code{bpe = "mean"} would
#' calculate the mean value of the posterior distribution.
#' @param bpe.style For \strong{Stan}-models (fitted with the \pkg{rstanarm}- or
#' \pkg{brms}-package), the Bayesian point estimate is indicated as a small,
#' vertical line by default. Use \code{bpe.style = "dot"} to plot a dot
#' instead of a line for the point estimate.
#' @param bpe.color Character vector, indicating the color of the Bayesian
#' point estimate. Setting \code{bpe.color = NULL} will inherit the color
#' from the mapped aesthetic to match it with the geom's color.
#' @param ci.style Character vector, defining whether inner and outer intervals
#' for Bayesion models are shown in boxplot-style (\code{"whisker"}) or in
#' bars with different alpha-levels (\code{"bar"}).
#' @param ... Other arguments, passed down to various functions. Here is a list
#' of supported arguments and their description in detail.
#' \describe{
#' \item{\code{prob.inner} and \code{prob.outer}}{For \strong{Stan}-models
#' (fitted with the \pkg{rstanarm}- or \pkg{brms}-package) and coefficients
#' plot-types, you can specify numeric values between 0 and 1 for
#' \code{prob.inner} and \code{prob.outer}, which will then be used as inner
#' and outer probabilities for the uncertainty intervals (HDI). By default,
#' the inner probability is 0.5 and the outer probability is 0.89 (unless
#' \code{ci.lvl} is specified - in this case, \code{ci.lvl} is used as outer
#' probability).
#' }
#' \item{\code{size.inner}}{For \strong{Stan}-models and \emph{Coefficients}
#' plot-types, you can specify the width of the bar for the inner
#' probabilities. Default is \code{0.1}. Setting \code{size.inner = 0}
#' removes the inner probability regions.
#' }
#' \item{\code{width}, \code{alpha}, and \code{scale}}{Passed
#' down to \code{geom_errorbar()} or \code{geom_density_ridges()}, for
#' forest or diagnostic plots.
#' }
#' \item{\code{width}, \code{alpha}, \code{dot.alpha}, \code{dodge} and \code{log.y}}{Passed
#' down to \code{\link[ggeffects]{plot.ggeffects}} for \emph{Marginal Effects}
#' plots.
#' }
#' \item{\code{show.loess}}{Logical, for diagnostic plot-types \code{"slope"}
#' and \code{"resid"}, adds (or hides) a loess-smoothed line to the plot.
#' }
#' \item{\emph{Marginal Effects} plot-types}{When plotting marginal effects,
#' arguments are also passed down to \code{\link[ggeffects]{ggpredict}},
#' \code{\link[ggeffects]{ggeffect}} or \code{\link[ggeffects]{plot.ggeffects}}.
#' }
#' \item{Case conversion of labels}{For case conversion of labels (see argument
#' \code{case}), arguments \code{sep_in} and \code{sep_out} will be passed
#' down to \code{snakecase::to_any_case()}. This only
#' applies to automatically retrieved term labels, \emph{not} if
#' term labels are provided by the \code{axis.labels}-argument.
#' }
#' }
#'
#' @return
#' Depending on the plot-type, \code{plot_model()} returns a
#' \code{ggplot}-object or a list of such objects. \code{get_model_data}
#' returns the associated data with the plot-object as tidy data frame, or
#' (depending on the plot-type) a list of such data frames.
#'
#' @details
#' \subsection{Different Plot Types}{
#' \describe{
#' \item{\code{type = "std"}}{Plots standardized estimates. See details below.}
#' \item{\code{type = "std2"}}{Plots standardized estimates, 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.
#' }
#' \item{\code{type = "pred"}}{Plots estimated marginal means (or marginal effects).
#' Simply wraps \code{\link[ggeffects]{ggpredict}}. See also
#' \href{https://strengejacke.github.io/sjPlot/articles/plot_marginal_effects.html}{this package-vignette}.
#' }
#' \item{\code{type = "eff"}}{Plots estimated marginal means (or marginal effects).
#' Simply wraps \code{\link[ggeffects]{ggeffect}}. See also
#' \href{https://strengejacke.github.io/sjPlot/articles/plot_marginal_effects.html}{this package-vignette}.
#' }
#' \item{\code{type = "int"}}{A shortcut for marginal effects plots, where
#' interaction terms are automatically detected and used as
#' \code{terms}-argument. Furthermore, if the moderator variable (the second
#' - and third - term in an interaction) is continuous, \code{type = "int"}
#' automatically chooses useful values based on the \code{mdrt.values}-argument,
#' which are passed to \code{terms}. Then, \code{\link[ggeffects]{ggpredict}}
#' is called. \code{type = "int"} plots the interaction term that appears
#' first in the formula along the x-axis, while the second (and possibly
#' third) variable in an interaction is used as grouping factor(s)
#' (moderating variable). Use \code{type = "pred"} or \code{type = "eff"}
#' and specify a certain order in the \code{terms}-argument to indicate
#' which variable(s) should be used as moderator. See also
#' \href{https://strengejacke.github.io/sjPlot/articles/plot_interactions.html}{this package-vignette}.
#' }
#' \item{\code{type = "slope"} and \code{type = "resid"}}{Simple diagnostic-plots,
#' where a linear model for each single predictor is plotted against the
#' response variable, or the model's residuals. Additionally, a loess-smoothed
#' line is added to the plot. The main purpose of these plots is to check whether
#' the relationship between outcome (or residuals) and a predictor is roughly
#' linear or not. Since the plots are based on a simple linear regression with
#' only one model predictor at the moment, the slopes (i.e. coefficients) may
#' differ from the coefficients of the complete model.
#' }
#' \item{\code{type = "diag"}}{For \strong{Stan-models}, plots the prior versus
#' posterior samples. For \strong{linear (mixed) models}, plots for
#' multicollinearity-check (Variance Inflation Factors), QQ-plots,
#' checks for normal distribution of residuals and homoscedasticity
#' (constant variance of residuals) are shown. For \strong{generalized
#' linear mixed models}, returns the QQ-plot for random effects.
#' }
#' }
#' }
#' \subsection{Standardized Estimates}{
#' Default standardization is done by completely refitting the model on the
#' standardized data. Hence, this approach is equal to standardizing the
#' variables before fitting the model, which is particularly recommended for
#' complex models that include interactions or transformations (e.g., polynomial
#' or spline terms). When \code{type = "std2"}, standardization of estimates
#' follows \href{http://www.stat.columbia.edu/~gelman/research/published/standardizing7.pdf}{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.
#' }
#'
#' @references
#' Gelman A (2008) "Scaling regression inputs by dividing by two
#' standard deviations." \emph{Statistics in Medicine 27: 2865-2873.}
#' \url{http://www.stat.columbia.edu/~gelman/research/published/standardizing7.pdf}
#' \cr \cr
#' Aiken and West (1991). Multiple Regression: Testing and Interpreting Interactions.
#'
#' @examples
#' # prepare data
#' library(sjmisc)
#' data(efc)
#' efc <- to_factor(efc, c161sex, e42dep, c172code)
#' m <- lm(neg_c_7 ~ pos_v_4 + c12hour + e42dep + c172code, data = efc)
#'
#' # simple forest plot
#' plot_model(m)
#'
#' # grouped coefficients
#' plot_model(m, group.terms = c(1, 2, 3, 3, 3, 4, 4))
#'
#' # keep only selected terms in the model: pos_v_4, the
#' # levels 3 and 4 of factor e42dep and levels 2 and 3 for c172code
#' plot_model(m, terms = c("pos_v_4", "e42dep [3,4]", "c172code [2,3]"))
#'
#' # multiple plots, as returned from "diagnostic"-plot type,
#' # can be arranged with 'plot_grid()'
#' \dontrun{
#' p <- plot_model(m, type = "diag")
#' plot_grid(p)}
#'
#' # plot random effects
#' if (require("lme4")) {
#' m <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
#' plot_model(m, type = "re")
#'
#' # plot marginal effects
#' plot_model(m, type = "pred", terms = "Days")
#' }
#' # plot interactions
#' \dontrun{
#' m <- glm(
#' tot_sc_e ~ c161sex + c172code * neg_c_7,
#' data = efc,
#' family = poisson()
#' )
#' # type = "int" automatically selects groups for continuous moderator
#' # variables - see argument 'mdrt.values'. The following function call is
#' # identical to:
#' # plot_model(m, type = "pred", terms = c("c172code", "neg_c_7 [7,28]"))
#' plot_model(m, type = "int")
#'
#' # switch moderator
#' plot_model(m, type = "pred", terms = c("neg_c_7", "c172code"))
#' # same as
#' # ggeffects::ggpredict(m, terms = c("neg_c_7", "c172code"))}
#'
#' # plot Stan-model
#' \dontrun{
#' if (require("rstanarm")) {
#' data(mtcars)
#' m <- stan_glm(mpg ~ wt + am + cyl + gear, data = mtcars, chains = 1)
#' plot_model(m, bpe.style = "dot")
#' }}
#'
#' @importFrom insight model_info find_predictors
#' @importFrom sjmisc word_wrap str_contains
#' @importFrom sjlabelled response_labels term_labels
#' @importFrom dplyr if_else n_distinct
#' @importFrom graphics plot
#' @importFrom ggeffects ggpredict ggeffect
#' @importFrom stats terms
#'
#' @export
plot_model <- function(model,
type = c("est", "re", "eff", "emm", "pred", "int", "std", "std2", "slope", "resid", "diag"),
transform,
terms = NULL,
sort.est = NULL,
rm.terms = NULL,
group.terms = NULL,
order.terms = NULL,
pred.type = c("fe", "re"),
mdrt.values = c("minmax", "meansd", "zeromax", "quart", "all"),
ri.nr = NULL,
title = NULL,
axis.title = NULL,
axis.labels = NULL,
legend.title = NULL,
wrap.title = 50,
wrap.labels = 25,
axis.lim = NULL,
grid.breaks = NULL,
ci.lvl = NULL,
se = NULL,
robust = FALSE,
vcov.fun = NULL,
vcov.type = c("HC3", "const", "HC", "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5"),
vcov.args = NULL,
colors = "Set1",
show.intercept = FALSE,
show.values = FALSE,
show.p = TRUE,
show.data = FALSE,
show.legend = TRUE,
show.zeroinf = TRUE,
value.offset = NULL,
value.size,
jitter = NULL,
digits = 2,
dot.size = NULL,
line.size = NULL,
vline.color = NULL,
p.threshold = c(0.05, 0.01, 0.001),
p.adjust = NULL,
grid,
case,
auto.label = TRUE,
prefix.labels = c("none", "varname", "label"),
bpe = "median",
bpe.style = "line",
bpe.color = "white",
ci.style = c("whisker", "bar"),
...
) {
type <- match.arg(type)
pred.type <- match.arg(pred.type)
mdrt.values <- match.arg(mdrt.values)
prefix.labels <- match.arg(prefix.labels)
vcov.type <- match.arg(vcov.type)
ci.style <- match.arg(ci.style)
# if we prefix labels, use different default for case conversion,
# else the separating white spaces after colon are removed.
if (missing(case)) {
if (prefix.labels == "none")
case <- "parsed"
else
case <- NULL
}
if (isTRUE(robust)) {
vcov.type <- "HC3"
vcov.fun <- "vcovHC"
}
# check se-argument
vcov.fun <- check_se_argument(se = vcov.fun, type = type)
# get info on model family
fam.info <- insight::model_info(model)
if (insight::is_multivariate(model))
fam.info <- fam.info[[1]]
# check whether estimates should be transformed or not
if (missing(transform)) {
if (fam.info$is_linear)
transform <- NULL
else
transform <- "exp"
}
# get titles and labels for axis ----
# this is not appropriate when plotting random effects,
# so retrieve labels only for other plot types
if (type %in% c("est", "std", "std2") && isTRUE(auto.label)) {
# get labels of dependent variables, and wrap them if too long
if (is.null(title)) title <- sjlabelled::response_labels(model, case = case, mv = fam.info$is_multivariate, ...)
title <- sjmisc::word_wrap(title, wrap = wrap.title)
# labels for axis with term names
if (is.null(axis.labels)) {
term_labels <- sjlabelled::term_labels(model, case = case, prefix = prefix.labels, ...)
if (.labelled_model_data(model) || is.stan(model)) axis.labels <- term_labels
}
axis.labels <- sjmisc::word_wrap(axis.labels, wrap = wrap.labels)
# title for axis with estimate values
if (is.null(axis.title)) axis.title <- sjmisc::word_wrap(estimate_axis_title(fit = model, axis.title = axis.title, type = type, transform = transform, include.zeroinf = TRUE), wrap = wrap.title)
axis.title <- sjmisc::word_wrap(axis.title, wrap = wrap.labels)
}
# check nr of estimates. if only one, plot slope
if (type == "est" &&
length(insight::find_predictors(model, component = "conditional", flatten = TRUE)) == 1 &&
length(insight::find_predictors(model, component = "instruments", flatten = TRUE)) == 0 &&
fam.info$is_linear && one_par(model)) type <- "slope"
# set some default options for stan-models, which are not
# available or appropriate for these
if (is.stan(model)) {
# no p-values
show.p <- FALSE
# no standardized coefficients
if (type %in% c("std", "std2", "slope")) type <- "est"
}
# set defaults for arguments, depending on model ----
if (is.null(ci.lvl)) ci.lvl <- dplyr::if_else(is.stan(model), .89, .95)
if (is.null(dot.size)) dot.size <- dplyr::if_else(is.stan(model), 1, 2.5)
if (is.null(line.size)) line.size <- dplyr::if_else(is.stan(model), .7, .7)
if (is.null(value.offset)) value.offset <- dplyr::if_else(is.stan(model), .25, .15)
# check if plot-type is applicable
if (type == "slope" && !fam.info$is_linear) {
type <- "est"
message("Plot-type \"slope\" only available for linear models. Using `type = \"est\"` now.")
}
if (type %in% c("est", "std", "std2") || (is.stan(model) && type == "re")) {
# plot estimates ----
p <- plot_type_est(
type = type,
ci.lvl = ci.lvl,
se = se,
tf = transform,
model = model,
terms = terms,
group.terms = group.terms,
rm.terms = rm.terms,
sort.est = sort.est,
title = title,
axis.title = axis.title,
axis.labels = axis.labels,
axis.lim = axis.lim,
grid.breaks = grid.breaks,
show.intercept = show.intercept,
show.values = show.values,
show.p = show.p,
value.offset = value.offset,
digits = digits,
geom.colors = colors,
geom.size = dot.size,
line.size = line.size,
order.terms = order.terms,
vline.color = vline.color,
value.size = value.size,
bpe = bpe,
bpe.style = bpe.style,
bpe.color = bpe.color,
facets = grid,
show.zeroinf = show.zeroinf,
p.threshold = p.threshold,
vcov.fun = vcov.fun,
vcov.type = vcov.type,
vcov.args = vcov.args,
ci.style = ci.style,
p_adjust = p.adjust,
...
)
} else if (type == "re") {
# plot random effects ----
p <- plot_type_ranef(
model = model,
ri.nr = ri.nr,
ci.lvl = ci.lvl,
se = se,
tf = transform,
sort.est = sort.est,
title = title,
axis.labels = axis.labels,
axis.lim = axis.lim,
grid.breaks = grid.breaks,
show.values = show.values,
value.offset = value.offset,
digits = digits,
facets = grid,
geom.colors = colors,
geom.size = dot.size,
line.size = line.size,
vline.color = vline.color,
value.size = value.size,
bpe.color = bpe.color,
ci.style = ci.style,
...
)
} else if (type %in% c("pred", "eff", "emm")) {
# plot marginal effects ----
p <- plot_type_eff(
type = type,
model = model,
terms = terms,
ci.lvl = ci.lvl,
pred.type = pred.type,
facets = grid,
show.data = show.data,
jitter = jitter,
geom.colors = colors,
axis.title = axis.title,
title = title,
legend.title = legend.title,
axis.lim = axis.lim,
case = case,
show.legend = show.legend,
dot.size = dot.size,
line.size = line.size,
...
)
} else if (type == "int") {
# plot interaction terms ----
p <- plot_type_int(
model = model,
mdrt.values = mdrt.values,
ci.lvl = ci.lvl,
pred.type = pred.type,
facets = grid,
show.data = show.data,
jitter = jitter,
geom.colors = colors,
axis.title = axis.title,
title = title,
legend.title = legend.title,
axis.lim = axis.lim,
case = case,
show.legend = show.legend,
dot.size = dot.size,
line.size = line.size,
...
)
} else if (type %in% c("slope", "resid")) {
# plot slopes of estimates ----
p <- plot_type_slope(
model = model,
terms = terms,
rm.terms = rm.terms,
ci.lvl = ci.lvl,
colors = colors,
title = title,
show.data = show.data,
jitter = jitter,
facets = grid,
axis.title = axis.title,
case = case,
useResiduals = type == "resid",
...
)
} else if (type == "diag") {
# plot diagnostic plots ----
if (is.stan(model)) {
p <- plot_diag_stan(
model = model,
geom.colors = colors,
axis.lim = axis.lim,
facets = grid,
axis.labels = axis.labels,
...
)
} else if (fam.info$is_linear) {
p <- plot_diag_linear(
model = model,
geom.colors = colors,
dot.size = dot.size,
line.size = line.size,
...
)
} else {
p <- plot_diag_glm(
model = model,
geom.colors = colors,
dot.size = dot.size,
line.size = line.size,
...
)
}
}
p
}
#' @importFrom purrr map
#' @rdname plot_model
#' @export
get_model_data <- function(model,
type = c("est", "re", "eff", "pred", "int", "std", "std2", "slope", "resid", "diag"),
transform,
terms = NULL,
sort.est = NULL,
rm.terms = NULL,
group.terms = NULL,
order.terms = NULL,
pred.type = c("fe", "re"),
ri.nr = NULL,
ci.lvl = NULL,
colors = "Set1",
grid,
case = "parsed",
digits = 2,
...) {
p <- plot_model(
model = model,
type = type,
transform = transform,
terms = terms,
sort.est = sort.est,
rm.terms = rm.terms,
group.terms = group.terms,
order.terms = order.terms,
pred.type = pred.type,
ri.nr = ri.nr,
ci.lvl = ci.lvl,
colors = colors,
grid = grid,
case = case,
digits = digits,
auto.label = FALSE,
...
)
if (inherits(p, "list"))
purrr::map(p, ~ .x$data)
else
p$data
}
#' @importFrom insight has_intercept
one_par <- function(model) {
tryCatch(
{
length(stats::coef(model)) < 2 & !insight::has_intercept(model)
},
error = function(x) { FALSE }
)
}