Added log ratio and log odds ratio contrasts

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)
 

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)
+
+})
+