I'm not sure why ggpredict
isn't adding confidence intervals. However, with some additional coding you can roll your own ggplot effects plot using the output of ggpredict
and effects
. (Ideally, you would generate the predictions and confidence bands directly from the model object, but it was quicker to just get them from ggpredict and effects.)
library(tidyverse)
library(patchwork)
# Function to combine output of effects and ggpredict into a single
# tidy data frame
ci.fnc = function(effects.obj, ggpredict.obj) {
# Capture the name of the effects term
col = sym(names(effects.obj$x))
# Generate a tidy data frame with lower and upper confidence intervals
cis = with(effects.obj, cbind(x, upper.prob, lower.prob)) %>%
pivot_longer(cols=-col) %>%
separate(name, into=c('type', 'response.level'), sep=".prob.X") %>%
mutate(response.level=gsub("\\.", " ", response.level)) %>%
spread(type, value)
# Rename the term column from "x" to the actual term name
ggpredict.obj = ggpredict.obj %>% rename(!!col := x)
# joint the predictions and confidence intervals into a single data frame
cis = left_join(cis, ggpredict.obj) %>% as.data.frame()
return(cis)
}
Now, generate the effects plots:
# Numeric x-axis variable
count = ci.fnc(pred1e, pred1) %>%
ggplot(aes(count, predicted)) +
geom_ribbon(aes(ymin=L, ymax=U), fill="blue", alpha=0.15) +
geom_line(aes(group=response.level), colour="blue", size=1) +
facet_grid(. ~ response.level) +
theme_bw() +
theme(plot.title=element_text(hjust=0.5),
panel.grid.minor.x=element_blank()) +
scale_y_continuous(limits=c(0,1), expand=c(0,0)) +
scale_x_continuous(breaks=0:10) +
labs(y="probability", x="count", title="Count effect plot")
# Categorical x-axis variable
case = ci.fnc(pred2e, pred2) %>%
ggplot(aes(case, predicted)) +
geom_line(aes(group=response.level), colour="grey70") +
geom_pointrange(aes(ymin=L, ymax=U), fatten=7, shape=21,
fill="blue", colour="red", stroke=0) +
facet_grid(. ~ response.level) +
theme_bw() +
theme(plot.title=element_text(hjust=0.5)) +
scale_y_continuous(limits=c(0,1), expand=c(0,0)) +
labs(y="probability", x="case", title="Case effect plot")
# Lay out both plots together
count / case
To create the ci.fnc
function, I looked at the structure of the objects returned by effects
and ggpredict
. For example, ggpredict
returns a data frame (run class(pred2)
or str(pred2)
), and we print it using as.data.frame
to override the default print method for ggpredict
objects.
as.data.frame(pred2)
x predicted response.level group
1 N 0.3614441 1 Better 1
2 N 0.2951495 2 Medium 1
3 N 0.3434064 3 Worse 1
4 Y 0.3448997 1 Better 1
5 Y 0.3290675 2 Medium 1
6 Y 0.3260327 3 Worse 1
effects
returns a list. Below is the output of str(pred2e)
. You can see that the x
element contains the values of term variable (the x-axis value in each plot. Then if you scroll to the bottom, you can see that the lower.prob
and upper.prob
elements contain the confidence limits. So, ci.fnc
is set up to extract these. You can also manually extract these values. For example:
pred2e$lower.prob
pred2e$upper.prob
pred2e$upper.prob
pred2e$x
with(pred2e, cbind(x, upper.prob, lower.prob))
str(pred2e)
List of 19
$ term : chr "case"
$ formula :Class 'formula' language rating ~ count + case
.. ..- attr(*, "variables")= language list(rating, count, case)
.. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:3] "rating" "count" "case"
.. .. .. ..$ : chr [1:2] "count" "case"
.. ..- attr(*, "term.labels")= chr [1:2] "count" "case"
.. ..- attr(*, "order")= int [1:2] 1 1
.. ..- attr(*, "intercept")= int 1
.. ..- attr(*, "response")= int 1
.. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
.. ..- attr(*, "predvars")= language list(rating, count, case)
.. ..- attr(*, "dataClasses")= Named chr [1:3] "factor" "numeric" "factor"
.. .. ..- attr(*, "names")= chr [1:3] "rating" "count" "case"
$ response : chr "rating"
$ y.levels : chr [1:3] "1 Better" "2 Medium" "3 Worse"
$ variables :List of 1
..$ case:List of 3
.. ..$ name : chr "case"
.. ..$ is.factor: logi TRUE
.. ..$ levels : chr [1:2] "N" "Y"
$ x :'data.frame': 2 obs. of 1 variable:
..$ case: Factor w/ 2 levels "N","Y": 1 2
$ model.matrix : num [1:2, 1:3] 1 1 3.56 3.56 0 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "1" "2"
.. ..$ : chr [1:3] "(Intercept)" "count" "caseY"
..- attr(*, "assign")= int [1:3] 0 1 2
..- attr(*, "contrasts")=List of 1
.. ..$ case: chr "contr.treatment"
$ data :'data.frame': 18 obs. of 3 variables:
..$ rating: Factor w/ 3 levels "1 Better","2 Medium",..: 1 1 1 2 2 2 3 3 3 1 ...
..$ count : num [1:18] 2 0 5 8 10 3 2 1 0 0 ...
..$ case : Factor w/ 2 levels "N","Y": 2 1 2 2 1 2 1 2 1 1 ...
..- attr(*, "terms")=Classes 'terms', 'formula' language rating ~ count + case + (count + case)
.. .. ..- attr(*, "variables")= language list(rating, count, case)
.. .. ..- attr(*, "factors")= int [1:3, 1:2] 0 1 0 0 0 1
.. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. ..$ : chr [1:3] "rating" "count" "case"
.. .. .. .. ..$ : chr [1:2] "count" "case"
.. .. ..- attr(*, "term.labels")= chr [1:2] "count" "case"
.. .. ..- attr(*, "order")= int [1:2] 1 1
.. .. ..- attr(*, "intercept")= int 1
.. .. ..- attr(*, "response")= int 1
.. .. ..- attr(*, ".Environment")=<environment: 0x7fddc4d92df8>
.. .. ..- attr(*, "predvars")= language list(rating, count, case)
.. .. ..- attr(*, "dataClasses")= Named chr [1:3] "factor" "numeric" "factor"
.. .. .. ..- attr(*, "names")= chr [1:3] "rating" "count" "case"
$ discrepancy : num 0
$ model : chr "multinom"
$ prob : num [1:2, 1:3] 0.361 0.345 0.295 0.329 0.343 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "prob.X1.Better" "prob.X2.Medium" "prob.X3.Worse"
$ logit : num [1:2, 1:3] -0.569 -0.642 -0.871 -0.712 -0.648 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "logit.X1.Better" "logit.X2.Medium" "logit.X3.Worse"
$ se.prob : num [1:2, 1:3] 0.17 0.165 0.167 0.168 0.17 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "se.prob.X1.Better" "se.prob.X2.Medium" "se.prob.X3.Worse"
$ se.logit : num [1:2, 1:3] 0.735 0.729 0.803 0.76 0.754 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "se.logit.X1.Better" "se.logit.X2.Medium" "se.logit.X3.Worse"
$ lower.logit : num [1:2, 1:3] -2.01 -2.07 -2.44 -2.2 -2.13 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "L.logit.X1.Better" "L.logit.X2.Medium" "L.logit.X3.Worse"
$ upper.logit : num [1:2, 1:3] 0.872 0.787 0.704 0.778 0.829 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "U.logit.X1.Better" "U.logit.X2.Medium" "U.logit.X3.Worse"
$ lower.prob : num [1:2, 1:3] 0.1181 0.112 0.0798 0.0995 0.1067 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "L.prob.X1.Better" "L.prob.X2.Medium" "L.prob.X3.Worse"
$ upper.prob : num [1:2, 1:3] 0.705 0.687 0.669 0.685 0.696 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:3] "U.prob.X1.Better" "U.prob.X2.Medium" "U.prob.X3.Worse"
$ confidence.level: num 0.95
- attr(*, "class")= chr "effpoly"