# prepare dummy variables for binary logistic regression
swiss$y1 <- ifelse(swiss$Fertility < median(swiss$Fertility), 0, 1)
swiss$y2 <- ifelse(swiss$Infant.Mortality < median(swiss$Infant.Mortality), 0, 1)
swiss$y3 <- ifelse(swiss$Agriculture < median(swiss$Agriculture), 0, 1)
# Now fit the models. Note that both models share the same predictors
# and only differ in their dependent variable (y1, y2 and y3)
fitOR1 <- glm(y1 ~ Education + Examination + Catholic, data = swiss,
family = binomial(link = "logit"))
fitOR2 <- glm(y2 ~ Education + Examination + Catholic, data = swiss,
family = binomial(link = "logit"))
fitOR3 <- glm(y3 ~ Education + Examination + Catholic, data = swiss,
family = binomial(link = "logit"))
## Not run:
# # open HTML-table in RStudio Viewer Pane or web browser
# sjt.glm(fitOR1, fitOR2,
# depvar.labels = c("Fertility", "Infant Mortality"),
# pred.labels = c("Education", "Examination", "Catholic"),
# ci.hyphen = " to ")
#
# # open HTML-table in RStudio Viewer Pane or web browser,
# # integrate CI in OR column
# sjt.glm(fitOR1, fitOR2, fitOR3,
# pred.labels = c("Education", "Examination", "Catholic"),
# separate.ci.col = FALSE)
#
# # open HTML-table in RStudio Viewer Pane or web browser,
# # indicating p-values as numbers and printing CI in a separate column
# sjt.glm(fitOR1, fitOR2, fitOR3,
# depvar.labels = c("Fertility", "Infant Mortality", "Agriculture"),
# pred.labels = c("Education", "Examination", "Catholic"))
#
# # --------------------------------------------
# # User defined style sheet
# # --------------------------------------------
# sjt.glm(fitOR1, fitOR2, fitOR3,
# depvar.labels = c("Fertility", "Infant Mortality", "Agriculture"),
# pred.labels = c("Education", "Examination", "Catholic"),
# show.header = TRUE,
# CSS = list(css.table = "border: 2px solid;",
# css.tdata = "border: 1px solid;",
# css.depvarhead = "color:#003399;"))
#
# # --------------------------------------------
# # Compare models with different link functions,
# # but same predictors and response
# # --------------------------------------------
# library(sjmisc)
# # load efc sample data
# data(efc)
# # dichtomozize service usage by "service usage yes/no"
# efc$services <- sjmisc::dicho(efc$tot_sc_e, 0, as.num = TRUE)
# # fit 3 models with different link-functions
# fit1 <- glm(services ~ neg_c_7 + c161sex + e42dep,
# data = efc, family = binomial(link = "logit"))
# fit2 <- glm(services ~ neg_c_7 + c161sex + e42dep,
# data = efc, family = binomial(link = "probit"))
# fit3 <- glm(services ~ neg_c_7 + c161sex + e42dep,
# data = efc, family = poisson(link = "log"))
#
# # compare models
# sjt.glm(fit1, fit2, fit3, show.aic = TRUE, show.family = TRUE)
#
# # --------------------------------------------
# # Change style of p-values and CI-appearance
# # --------------------------------------------
# # open HTML-table in RStudio Viewer Pane or web browser,
# # table indicating p-values as stars
# sjt.glm(fit1, fit2, fit3, p.numeric = FALSE,
# show.aic = TRUE, show.family = TRUE)
#
# # open HTML-table in RStudio Viewer Pane or web browser,
# # indicating p-values as stars and integrate CI in OR column
# sjt.glm(fit1, fit2, fit3, p.numeric = FALSE, separate.ci.col = FALSE,
# show.aic = TRUE, show.family = TRUE, show.r2 = TRUE)
#
# # ----------------------------------
# # automatic grouping of predictors
# # ----------------------------------
# library(sjmisc)
# # load efc sample data
# data(efc)
# # dichtomozize service usage by "service usage yes/no"
# efc$services <- sjmisc::dicho(efc$tot_sc_e, 0, as.num = TRUE)
# # make dependency categorical
# efc$e42dep <- to_factor(efc$e42dep)
# # fit model with "grouped" predictor
# fit <- glm(services ~ neg_c_7 + c161sex + e42dep, data = efc)
#
# # automatic grouping of categorical predictors
# sjt.glm(fit)
#
# # ----------------------------------
# # compare models with different predictors
# # ----------------------------------
# fit2 <- glm(services ~ neg_c_7 + c161sex + e42dep + c12hour, data = efc)
# fit3 <- glm(services ~ neg_c_7 + c161sex + e42dep + c12hour + c172code,
# data = efc)
#
# # print models with different predictors
# sjt.glm(fit, fit2, fit3)
#
# efc$c172code <- to_factor(efc$c172code)
# fit2 <- glm(services ~ neg_c_7 + c161sex + c12hour, data = efc)
# fit3 <- glm(services ~ neg_c_7 + c161sex + c172code, data = efc)
#
# # print models with different predictors
# sjt.glm(fit, fit2, fit3, group.pred = FALSE)## End(Not run)
Run the code above in your browser using DataLab