library(ggeffects)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.2
data(efc)
head(efc)
## c12hour e15relat e16sex e17age e42dep c82cop1 c83cop2 c84cop3 c85cop4 c86cop5
## 1 16 2 2 83 3 3 2 2 2 1
## 2 148 2 2 88 3 3 3 3 3 4
## 3 70 1 2 82 3 2 2 1 4 1
## 4 168 1 2 67 4 4 1 3 1 1
## 5 168 2 2 84 4 3 2 1 2 2
## 6 16 2 2 85 4 2 2 3 3 3
## c87cop6 c88cop7 c89cop8 c90cop9 c160age c161sex c172code c175empl barthtot
## 1 1 2 3 3 56 2 2 1 75
## 2 1 3 2 2 54 2 2 1 75
## 3 1 1 4 3 80 1 1 0 35
## 4 1 1 2 4 69 1 2 0 0
## 5 2 1 4 4 47 2 2 0 25
## 6 2 2 1 1 56 1 2 1 60
## neg_c_7 pos_v_4 quol_5 resttotn tot_sc_e n4pstu nur_pst grp negc7d
## 1 12 12 14 0 4 0 NA child 1
## 2 20 11 10 4 0 0 NA child 1
## 3 11 13 7 0 1 2 2 spouse/partner 0
## 4 10 15 12 2 0 3 3 spouse/partner 0
## 5 12 15 19 2 1 2 2 child 1
## 6 19 9 8 1 3 2 2 child 1
library(skimr)
##
## Attaching package: 'skimr'
## The following object is masked from 'package:stats':
##
## filter
skim(efc)
## Skim summary statistics
## n obs: 908
## n variables: 28
## group variables:
##
## ── Variable type:factor ────────────────────────────────────────────────────────
## variable missing complete n n_unique top_counts
## grp 7 901 908 8 chi: 473, spo: 171, oth: 92, dau: 85
## negc7d 16 892 908 2 0: 496, 1: 396, NA: 16
## ordered
## FALSE
## FALSE
##
## ── Variable type:numeric ───────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50 p75 p100 hist
## barthtot 25 883 908 64.55 29.54 0 45 70 90 100 ▂▂▁▃▂▅▃▇
## c12hour 6 902 908 42.4 50.81 4 10 20 42.75 168 ▇▃▁▁▁▁▁▂
## c160age 7 901 908 53.46 13.35 18 44 54 63 89 ▁▃▅▇▇▆▂▁
## c161sex 7 901 908 1.76 0.43 1 2 2 2 2 ▂▁▁▁▁▁▁▇
## c172code 66 842 908 1.97 0.63 1 2 2 2 3 ▃▁▁▇▁▁▁▂
## c175empl 6 902 908 0.43 0.49 0 0 0 1 1 ▇▁▁▁▁▁▁▆
## c82cop1 7 901 908 3.12 0.58 1 3 3 3 4 ▁▁▁▁▁▇▁▃
## c83cop2 6 902 908 2.02 0.72 1 2 2 2 4 ▃▁▇▁▁▂▁▁
## c84cop3 6 902 908 1.63 0.87 1 1 1 2 4 ▇▁▃▁▁▁▁▁
## c85cop4 10 898 908 1.77 0.87 1 1 2 2 4 ▇▁▇▁▁▂▁▁
## c86cop5 6 902 908 1.39 0.67 1 1 1 2 4 ▇▁▃▁▁▁▁▁
## c87cop6 8 900 908 1.29 0.64 1 1 1 1 4 ▇▁▂▁▁▁▁▁
## c88cop7 8 900 908 1.92 0.91 1 1 2 2 4 ▇▁▇▁▁▂▁▂
## c89cop8 7 901 908 2.16 1.04 1 1 2 3 4 ▇▁▆▁▁▆▁▃
## c90cop9 20 888 908 2.93 0.96 1 2 3 4 4 ▂▁▆▁▁▇▁▇
## e15relat 7 901 908 2.85 2.08 1 2 2 4 8 ▃▇▁▂▁▁▁▂
## e16sex 7 901 908 1.67 0.47 1 1 2 2 2 ▃▁▁▁▁▁▁▇
## e17age 17 891 908 79.12 8.09 65 73 79 85 103 ▆▆▇▇▅▅▁▁
## e42dep 7 901 908 2.94 0.94 1 2 3 4 4 ▂▁▆▁▁▇▁▇
## n4pstu 9 899 908 1.02 1.09 0 0 1 2 4 ▇▃▁▃▁▂▁▁
## neg_c_7 16 892 908 11.85 3.86 7 9 11 14 28 ▇▇▃▃▂▁▁▁
## nur_pst 419 489 908 1.86 0.75 1 1 2 2 3 ▇▁▁▇▁▁▁▅
## pos_v_4 27 881 908 12.48 2.24 5 11 13 14 16 ▁▁▂▂▃▇▃▅
## quol_5 11 897 908 14.37 5.31 0 11 15 18 25 ▂▂▃▅▇▇▇▂
## resttotn 0 908 908 0.33 0.68 0 0 0 0 4 ▇▂▁▁▁▁▁▁
## tot_sc_e 0 908 908 1.01 1.28 0 0 1 1.25 9 ▇▂▁▁▁▁▁▁
fit <- lm(barthtot ~ c12hour + neg_c_7 + c161sex +
c172code, data = efc)
ggpredict(fit, terms = "c12hour")
## # Predicted values of Total score BARTHEL INDEX
## # x = average number of hours of care per week
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 75.44 | [73.26, 77.63]
## 20 | 70.38 | [68.56, 72.19]
## 45 | 64.05 | [62.39, 65.70]
## 65 | 58.98 | [57.16, 60.80]
## 85 | 53.91 | [51.71, 56.11]
## 105 | 48.85 | [46.15, 51.55]
## 125 | 43.78 | [40.52, 47.05]
## 170 | 32.38 | [27.73, 37.03]
##
## Adjusted for:
## * neg_c_7 = 11.84
## * c161sex = 1.76
## * c172code = 1.97
mydf <- ggpredict(fit, terms = c("c12hour", "c172code"))
ggplot(mydf, aes(x, predicted, colour = group)) + geom_line()+
theme_bw(base_size=16)
mydf <- ggpredict(fit, terms = c("c12hour", "c172code", "c161sex"))
mydf
## # Predicted values of Total score BARTHEL INDEX
## # x = average number of hours of care per week
##
## # c172code = low level of education
## # c161sex = [1] Male
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 73.95 | [69.35, 78.55]
## 45 | 62.56 | [58.23, 66.88]
## 85 | 52.42 | [47.90, 56.95]
## 170 | 30.89 | [24.85, 36.94]
##
## # c172code = intermediate level of education
## # c161sex = [1] Male
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 74.67 | [71.06, 78.29]
## 45 | 63.27 | [59.88, 66.66]
## 85 | 53.14 | [49.40, 56.89]
## 170 | 31.61 | [25.98, 37.24]
##
## # c172code = high level of education
## # c161sex = [1] Male
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 75.39 | [71.04, 79.74]
## 45 | 63.99 | [59.73, 68.26]
## 85 | 53.86 | [49.23, 58.49]
## 170 | 32.33 | [25.95, 38.71]
##
## # c172code = low level of education
## # c161sex = [2] Female
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 75.00 | [71.41, 78.59]
## 45 | 63.60 | [60.46, 66.74]
## 85 | 53.46 | [50.13, 56.80]
## 170 | 31.93 | [26.83, 37.04]
##
## # c172code = intermediate level of education
## # c161sex = [2] Female
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 75.71 | [73.31, 78.12]
## 45 | 64.32 | [62.42, 66.21]
## 85 | 54.18 | [51.81, 56.55]
## 170 | 32.65 | [27.94, 37.36]
##
## # c172code = high level of education
## # c161sex = [2] Female
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 76.43 | [72.89, 79.98]
## 45 | 65.03 | [61.68, 68.39]
## 85 | 54.90 | [51.16, 58.65]
## 170 | 33.37 | [27.70, 39.05]
##
## Adjusted for:
## * neg_c_7 = 11.84
ggplot(mydf, aes(x, predicted, colour = group)) +
geom_line() + theme_bw(base_size=18)+
facet_wrap(~facet)
mydf <- ggpredict(fit, terms = c("c12hour", "c172code", "c161sex", "neg_c_7"))
plot(mydf)
####
library(lme4)
## Loading required package: Matrix
data(sleepstudy)
skim(sleepstudy)
## Skim summary statistics
## n obs: 180
## n variables: 3
## group variables:
##
## ── Variable type:factor ────────────────────────────────────────────────────────
## variable missing complete n n_unique top_counts
## Subject 0 180 180 18 308: 10, 309: 10, 310: 10, 330: 10
## ordered
## FALSE
##
## ── Variable type:numeric ───────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50 p75 p100
## Days 0 180 180 4.5 2.88 0 2 4.5 7 9
## Reaction 0 180 180 298.51 56.33 194.33 255.38 288.65 336.75 466.35
## hist
## ▇▃▃▃▃▃▃▇
## ▂▇▇▆▆▃▁▁
m <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
pr <- ggpredict(m, "Days")
pr
## # Predicted values of Reaction
## # x = Days
##
## x | Predicted | 95% CI
## --------------------------------
## 0 | 251.41 | [238.03, 264.78]
## 1 | 261.87 | [248.57, 275.17]
## 2 | 272.34 | [258.44, 286.24]
## 3 | 282.81 | [267.71, 297.91]
## 5 | 303.74 | [284.96, 322.52]
## 6 | 314.21 | [293.17, 335.25]
## 7 | 324.68 | [301.21, 348.14]
## 9 | 345.61 | [316.94, 374.29]
##
## Adjusted for:
## * Subject = 0 (population-level)