Plots

library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(ggplot2)
theme_set(theme_bw(base_size=16))
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
## 1      12      12     14        0        4      0      NA
## 2      20      11     10        4        0      0      NA
## 3      11      13      7        0        1      2       2
## 4      10      15     12        2        0      3       3
## 5      12      15     19        2        1      2       2
## 6      19       9      8        1        3      2       2

plot_grpfrq(efc$e42dep, efc$c172code, geom.colors = "gs")

# create binrary response
y <- ifelse(efc$neg_c_7 < median(na.omit(efc$neg_c_7)), 0, 1)

# create data frame for fitting model
df <- data.frame(
  y = to_factor(y),
  sex = to_factor(efc$c161sex),
  dep = to_factor(efc$e42dep),
  barthel = efc$barthtot,
  education = to_factor(efc$c172code)
)
head(df)

##   y sex dep barthel education
## 1 1   2   3      75         2
## 2 1   2   3      75         2
## 3 1   1   3      35         1
## 4 0   1   4       0         2
## 5 1   2   4      25         2
## 6 1   1   4      60         2

# set variable label for response
set_label(df$y) <- "High Negative Impact"

# fit model
fit <- glm(y ~., data = df, family = binomial(link = "logit"))

# plot marginal effects
plot_model(
  fit, 
  type = "pred", 
  terms = c("barthel", "sex","dep"), 
  colors = "bw",
  ci.lvl = NA
)+ theme(title = element_blank())

plot_model(fit, colors = "black")+theme_bw(base_size=20) + theme(title = element_blank())