Nothing
#' performance_equiv function
#'
#' This function takes two logistic regression models \eqn{M_A, M_B},
#' test data, significance level \eqn{\alpha} and acceptable score
#' degradation \eqn{\delta_B}. It checks whether the models perform
#' equivalently on the test set and returns various figures.
#' @param model_a logistic regression model \eqn{M_A}
#' @param model_b logistic regression model \eqn{M_B}
#' @param test_data testing dataset
#' @param dv_index column number of the dependent variable
#' @param delta_B acceptable score degradation (defaults to 1.1)
#' @param alpha significance level \eqn{\alpha} (defaults to 0.05)
#' @return \describe{
#' \item{\code{equivalence}}{Are models \eqn{M_A,M_B} producing equivalent
#' Brier scores for the given test data? (boolean)}
#' \item{\code{brier_score_ac}}{\eqn{M_A} Brier score on the testing data}
#' \item{\code{brier_score_bc}}{\eqn{M_B} Brier score on the testing data}
#' \item{\code{diff_sd_l}}{SD of the lower Brier difference \eqn{BS^A-\delta_B^2BS^B}}
#' \item{\code{diff_sd_u}}{SD of the upper Brier difference \eqn{BS^A-\delta_B^{-2}BS^B}}
#' \item{\code{test_stat_l}}{\eqn{t_L} equivalence boundary for the test}
#' \item{\code{test_stat_u}}{\eqn{t_U} equivalence boundary for the test}
#' \item{\code{crit_val}}{a level-\eqn{\alpha} critical value for the test}
#' \item{\code{delta_B}}{Calculated equivalence parameter}
#' \item{\code{p_value_l}}{P-value for \eqn{t_L}}
#' \item{\code{p_value_u}}{P-value for \eqn{t_U}}
#' }
#' @export
#' @importFrom stats predict.glm var qt pt
performance_equiv <- function(model_a, model_b, test_data,
dv_index, delta_B = 1.1, alpha = 0.05) {
m <- nrow(test_data)
test_y <- test_data[, dv_index]
eps_b <- delta_B^2
pi_ac <- predict.glm(model_a, test_data, type = "response")
pi_bc <- predict.glm(model_b, test_data, type = "response")
bs_ac <- brier_score(test_y, pi_ac)
bs_bc <- brier_score(test_y, pi_bc)
b_ac <- (pi_ac - test_y)^2
b_bc <- (pi_bc - test_y)^2
d_l <- b_bc - (b_ac / eps_b)
d_u <- b_bc - (b_ac * eps_b)
test_stat_l <- sqrt(m) * mean(d_l) / sqrt(var(d_l))
test_stat_u <- sqrt(m) * mean(d_u) / sqrt(var(d_u))
equivalence_threshold <- qt(1 - 0.5 * alpha, df = (m - 1))
bse_l <- (test_stat_l > equivalence_threshold)
bse_u <- (test_stat_u < -equivalence_threshold)
return(list(
equivalence = (bse_l && bse_u),
brier_score_ac = bs_ac,
brier_score_bc = bs_bc,
diff_sd_l = sd(d_l),
diff_sd_u = sd(d_u),
test_stat_l = test_stat_l,
test_stat_u = test_stat_u,
crit_val = equivalence_threshold,
epsilon_B = eps_b,
p_value_l = pt(test_stat_l, m - 1),
p_value_u = pt(test_stat_u, m - 1)
)
)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.