#' Utility function that will return numeric columns used in the regression model.
#'
#' @param glm_out the results of [stats::glm()].
#' @export
<- (glm_out) {
numeric_vars <- ((glm_out$, ))
numeric_vars <- (numeric_vars)[numeric_vars]
numeric_vars <- numeric_vars
((glm_out$,numeric_vars], FUN = (x) { ((x)) })) > 2]
(numeric_vars[numeric_vars %in% (glm_out$)-1]])
}
#' Visual plot to check the assumption that the predicted values are linearly
#' related to the predictors.
#'
#' @param glm_out the results of the logistic regression (i.e. the result from
#' \code{\link{glm}}).
#' @param vars character vector with the name of the independent variables to
#' include. If omitted, all numeric variables will be used.
#' @param n_groups the number of groups to use. By default groups will have a minimum of
#' six observations up to a maximum of 10 groups.
#' @param widths the relative widths of the two plots. If you only want the first
#' plot, use \code{widths = c(0, 1)}; if only the second plot, use
#' \code{widths = c(1, 0)}.
#' @param na.rm whether missing values should be removed.
#' @param x_lab label for the x-axis.
#' @param y_lab label for the y-axis.
#' @param x_limits limits for the x-axis.
#' @param coord_equal whether the x- and y-axes should be equal.
#' @param coord_flip whether the x- and y-axes should be flipped.
#' @param legend.position location of the legend.
#' @param varplot_nrow parameter passed to [ggplot2::facet_wrap()].
#' @param varplot_ncol parameter passed to [ggplot2::facet_wrap()].
#' @return a ggplot2 expression.
#' @export
#' @examples
#' study <- data.frame(
#' Hours=c(0.50,0.75,1.00,1.25,1.50,1.75,1.75,2.00,2.25,2.50,2.75,3.00,
#' 3.25,3.50,4.00,4.25,4.50,4.75,5.00,5.50),
#' Pass=c(0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,1,1,1)
#' )
#' lr.out <- glm(Pass ~ Hours, data = study,
#' family = binomial(link='logit'))
#' plot_linear_assumption_check(glm_out = lr.out, n_groups = 5)
<- (glm_out,
n_groups,
vars = (glm_out),
# boxTidwell = TRUE,
# boxTidwell_digits = 3,
# boxTidwell_position_x,
# boxTidwell_position_y = .85,
na.rm = ,
x_lab = 'Predicted Probabilities',
y_lab = (glm_out$)[1],
x_limits = (0, 1),
coord_equal = ,
coord_flip = ,
widths = ,
legend.position = 'bottom',
varplot_nrow = ,
varplot_ncol = ) {
boxTidwell <- # Not working reliability so leaving out for now.
(glm_out$$ != 'binomial') {
(('This function is intended for logistic regression, not ',
glm_out$$))
}
((vars) == 0) {
('No numeric variables found. Excluding variable plot.')
widths <- (1, 0)
}
((n_groups)) {
n_groups <- (((glm_out$) / 6), 10)
}
((x_limits)) {
x_limits <- ((, ))
}
probs <- (0, 1, 1 / n_groups)
predicted_prob <- (glm_out, = 'response')
y_var <- (glm_out$)[1]
Y <- glm_out$,y_var, = ]
p2 <-
((vars) > 0) {
X <- (glm_out$,vars, = ])
X_melt <- reshape2::melt(
(X, predicted_prob = predicted_prob),
id.var = 'predicted_prob'
)
p2 <- ggplot(X_melt, aes(x = predicted_prob, y = value, = variable)) +
geom_smooth(method = 'loess', = y ~ x) +
geom_point(alpha = 0.4) +
facet_wrap(~ variable,
scales = (coord_flip, "free_x", "free_y"),
= varplot_nrow,
= varplot_ncol) +
(x_limits) +
ylab() + xlab(x_lab)
}
<- (predicted_prob, probs = probs, na.rm = na.rm)
breaks <- (predicted_prob, breaks = , include.lowest = )
tab <- (
(Y, = breaks, mat = , = ),('group1', 'mean', 'se')],
x_pos = (predicted_prob, = breaks, mat = ),'median',=]
)
tab$ymin <- (tab,('mean', 'se')], 1, FUN = (x) { (x'mean'] - 1.96 * x'se'], 0) })
tab$ymax <- (tab,('mean', 'se')], 1, FUN = (x) { (x'mean'] + 1.96 * x'se'], 1) })
tab$color <- (tab$ymin <= tab$x_pos & tab$x_pos <= tab$ymax, 'p > 0.05', 'p <= 0.05')
p1 <- ggplot((predicted_prob = predicted_prob, Y = Y),
aes(x = predicted_prob, y = Y)) +
geom_vline(xintercept = , linetype = 2, alpha = 0.3) +
geom_jitter(alpha = 0.4, = 0, = 0.05) +
geom_segment( = tab, aes(x = x_pos, xend = x_pos, y = ymin, yend = ymax, color = color), linewidth = 4) +
geom_point( = tab, aes(x = x_pos, y = ), color = 'white', size = 2) +
geom_abline(slope = 1, intercept = 0, alpha = 0.3) +
scale_color_manual(values = ('#e41a1c', '#377eb8'),
limits = ('p \u2264 0.05', 'p > 0.05')) +
theme(legend.position = legend.position) +
xlab(x_lab) + ylab(y_lab) +
(x_limits) #+ ylim(c(0, 1))
((vars) > 0 & boxTidwell) {
numeric_vars <- (glm_out)
<- glm_out$,(y_var, numeric_vars)]
# out <- car::boxTidwell(as.formula(paste0(y_var, ' ~ ', paste0(numeric_vars, collapse = ' + '))),
# data = df) # Looking for a non-significant p-value
vars
out <- car::boxTidwell(y = glm_out$,y_var,=],
x1 = glm_out$,numeric_vars,=])
txt <- 'Box-Tidwell Test'
(i 1:(out$result)) {
txt <- (txt, '\n', (out$result)[i], ': ',
(out$result[i], digits = boxTidwell_digits))
}
p1 <- p1 + annotate("label",
x = ((boxTidwell_position_x),
((x_limits, predicted_prob), na.rm = ),
boxTidwell_position_x),
y = boxTidwell_position_y,
label = txt,
hjust = 1)
}
(coord_equal) {
p1 <- p1 + coord_equal()
}
(coord_flip) {
p1 <- p1 + coord_flip()
(!(p2)) {
p2 <- p2 + coord_flip()
}
}
(!(widths)) {
(widths[1] == 0) {
(p2)
} (widths[2] == 0) {
(p1)
}
} ((p2)) {
(p1)
}
# cowplot::plot_grid(
# p1 + theme_vs(),
# p2 + theme_vs(),
# nrow = 1)
egg::ggarrange(p1 + (),
p2 + (),
= 1,
widths = widths)
}
() { # Testing
(ggplot2)
('mtcars', envir = ())
mtcars_glm_out <- ( = vs ~ wt + disp, = , = (link='logit'))
(mtcars_glm_out)
(mtcars_glm_out)
y_var <- (mtcars_glm_out$)[1]
<- mtcars_glm_out$,(y_var, (mtcars_glm_out))]
# df[,y_var] <- df[,y_var] + 1
out <- car::boxTidwell(((y_var, ' ~ .')),
= ) # Looking for a non-significant p-value
boxTidwell((vs) ~ cyl + disp + hp, ~(am) + (gear, 2),
= )
# This example is from OpenIntro, 4th edition beginning on page 371
(openintro)
('resume', envir = ())
glm_out1 <- (received_callback ~ job_city + college_degree + years_experience +
honors + military + has_email_address + race + gender + job_ad_id,
= resume,
= (link = 'logit'))
(glm_out1)
(glm_out1,
x_limits = ,
coord_equal = ,
coord_flip = ,
legend.position = 'none',
varplot_nrow = 2,
widths = (1,2))
(glm_out1,
x_limits = ,
coord_equal = ,
widths = (1, 0))
(glm_out1,
x_limits = ,
coord_equal = ,
widths = (0, 1))
(glm_out1,
vars = ('years_experience'),
x_limits = ,
coord_equal = ,
widths = (0, 1))
# http://www.statisticalassociates.com/logistic10.htm
# The right-hand predictor side of the equation must be linear with the
# left-hand outcome side of the equation. You must test for linearity in the
# logit (in logistic regression the logit is the outcome side). This is
# commonly done with the Box-Tidwell Transformation (Test): Add to the
# logistic model interaction terms which are the crossproduct of each
# independent times its natural logarithm [(X)ln(X)]. If these terms are
# significant, then there is nonlinearity in the logit. This method is not
# sensitive to small nonlinearities.
y_var <- (glm_out1$)[1]
<- glm_out1$,(y_var, (glm_out1))]
out <- car::boxTidwell(((y_var, ' ~ .')),
= ) # Looking for a non-significant p-value
glm_out2 <- (received_callback ~ job_city + years_experience +
honors + military + has_email_address + race + gender,
= resume,
= (link = 'logit'))
(glm_out2)
(glm_out2,
x_limits = ,
coord_equal = ,
legend.position = 'none',
varplot_nrow = 2,
widths = (1,2))
admit <- ("https://stats.idre.ucla.edu/stat/data/binary.csv")
admit_glm_out <- (admit ~ gre + gpa + ,
= admit,
= (link = 'logit'))
(admit_glm_out, widths = (1, 2))
}
