diff --git a/R/robincar-contrast.R b/R/robincar-contrast.R index 1d6ed6f..803e9ec 100644 --- a/R/robincar-contrast.R +++ b/R/robincar-contrast.R @@ -5,12 +5,30 @@ diff <- function(est){ return(cont) } +odds_ratio <- function(est){ + base <- est[1] / (1-est[1]) + cont <- est[-1] / (1-est[-1]) / base + return(cont) +} + +log_odds_ratio <- function(est){ + base <- est[1] / (1-est[1]) + cont <- log(est[-1] / (1-est[-1]) / base) + return(cont) +} + ratio <- function(est){ base <- est[1] cont <- est[-1] / base return(cont) } +log_ratio <- function(est){ + base <- est[1] + cont <- log(est[-1] / base) + return(cont) +} + jacobian <- function (settings, est) { UseMethod("jacobian", settings) } @@ -30,6 +48,30 @@ jacobian.RATIO <- function(settings, est){ return(jacob) } +#' @exportS3Method +jacobian.LOGRATIO <- function(settings, est){ + col1 <- -1 / est[1] + mat1 <- diag(1/est[-1], length(est)-1) + jacob <- cbind(col1, mat1) + return(jacob) +} + +#' @exportS3Method +jacobian.ODDSRATIO <- function(settings, est){ + col1 <- -1 / est[1] ** 2 * est[-1] / (1-est[-1]) + mat1 <- diag((1-est[-1]) / est[-1] / (1-est[1]) ** 2, length(est)-1) + jacob <- cbind(col1, mat1) + return(jacob) +} + +#' @exportS3Method +jacobian.LOGODDSRATIO <- function(settings, est){ + col1 <- -1 / (est[1] * (1 - est[1])) + mat1 <- diag(1 / (est[-1] * (1 - est[-1])), length(est)-1) + jacob <- cbind(col1, mat1) + return(jacob) +} + #' @exportS3Method jacobian.CUSTOM <- function(settings, est){ if(!is.null(settings$dh)){ @@ -59,6 +101,18 @@ contrast.settings <- function(k, contrast_h, contrast_dh=NULL){ h <- ratio type <- "RATIO" name <- function(tl) paste0("treat ", tl[-1], " / ", tl[1]) + } else if(contrast_h == "log_ratio"){ + h <- log_ratio + type <- "LOGRATIO" + name <- function(tl) paste0("log(treat ", tl[-1], " / ", tl[1], ")") + } else if(contrast_h == "odds_ratio"){ + h <- odds_ratio + type <- "ODDSRATIO" + name <- function(tl) paste0("odds(treat ", tl[-1], ") / odds(", tl[1], ")") + } else if(contrast_h == "log_odds_ratio"){ + h <- log_odds_ratio + type <- "LOGODDSRATIO" + name <- function(tl) paste0("log(odds(treat ", tl[-1], ")) - log(odds(", tl[1], "))") } else { stop("Unrecognized contrast function name.") } @@ -80,6 +134,22 @@ contrast <- function(settings, treat, est, varcov){ # Get labels lab <- settings$name(treat) + # If using an odds ratio contrast, make sure + # give a warning if the estimates are not between 1 and 0. + if(class(settings)[2] %in% c("ODDSRATIO", "LOGODDSRATIO")){ + if(any(est < 0) || any(est > 1)){ + warning("Estimates are not between 0 and 1: are you sure you want an odds ratio?") + } + } + + # If using a ratio or odds ratio, give a warning + # that the performance is better with logs. + if(class(settings)[2] %in% c("ODDSRATIO", "RATIO")){ + warning("Using a ratio or odds ratio is not recommended. + Consider using log_ratio or log_odds_ratio for better + performance of the variance estimator.") + } + # Get transformed estimate c_est <- settings$h(est) diff --git a/tests/testthat/test-contrast.R b/tests/testthat/test-contrast.R index 7d86c8a..70bb580 100644 --- a/tests/testthat/test-contrast.R +++ b/tests/testthat/test-contrast.R @@ -26,3 +26,93 @@ test_that("Test contrast function ordering", { expect_equal(as.vector(fit.anova$contrast$result$estimate[1]), anova) }) + +test_that("Test warnings for ratio and odds ratio functions", { + + data <- RobinCar:::data_sim %>% + # recode A from 0 to "a", 1 to "b", 2 to "c" + mutate(A=factor(A, levels=0:2, labels=c("a", "b", "c"))) %>% + mutate(y_bin=ifelse(y > mean(y), 1, 0)) # create binary outcome + + warning_txt <- "Using a ratio or odds ratio is not recommended. + Consider using log_ratio or log_odds_ratio for better + performance of the variance estimator." + + # Test for linear model with ratio + expect_warning( + fit <- robincar_glm( + df=data, + response_col="y_bin", + treat_col="A", + formula= y_bin ~ A, + contrast_h="odds_ratio" + ), + warning_txt + ) + + # Test for linear model with ratio + expect_warning( + fit <- robincar_glm( + df=data, + response_col="y", + treat_col="A", + formula= y ~ A, + contrast_h="ratio" + ), + warning_txt + ) + +}) + + + +test_that("Test warnings for contrast functions that are not probabilities", { + # We want these to give warnings when someone + # specifies that they want a log odds ratio. + + data <- RobinCar:::data_sim %>% + # recode A from 0 to "a", 1 to "b", 2 to "c" + mutate(A=factor(A, levels=0:2, labels=c("a", "b", "c"))) %>% + mutate(y_bin=ifelse(y > mean(y), 1, 0)) # create binary outcome + + # Test for linear model with ratio + expect_warning( + fit <- robincar_glm( + df=data, + response_col="y", + treat_col="A", + formula= y ~ A, + contrast_h="log_odds_ratio" + ), + "Estimates are not between 0 and 1: are you sure you want an odds ratio?" + ) + +}) + +test_that("Compare log odds ratio SEs to GLM", { + # We want these to give warnings when someone + # specifies that they want a log odds ratio. + + data <- RobinCar:::data_sim %>% + # recode A from 0 to "a", 1 to "b", 2 to "c" + mutate(A=factor(A, levels=0:2, labels=c("a", "b", "c"))) %>% + mutate(y_bin=ifelse(y > mean(y), 1, 0)) # create binary outcome + + # Test for linear model with ratio + fit <- robincar_glm( + df=data, + response_col="y_bin", + treat_col="A", + formula= y_bin ~ A, + contrast_h="log_odds_ratio" + ) + + or <- with(data, glm(y_bin ~ A, family=binomial(link="logit"))) + + # Check that log odds ratios match + expect_equal( + unname(or$coefficients[2:3]), + unname(fit$contrast$result$estimate), tolerance=1e-6) + +}) +