R/epi.interaction.R

In epiR: Tools for the Analysis of Epidemiological Data

epi.interaction <- function(model, coef, param = c("product", "dummy"), conf.level = 0.95){

   N. <- 1 - ((1 - conf.level)/2)
   z <- qnorm(N., mean = 0, sd = 1)
  
   if (class(model)[1] != "glm" & class(model)[2] != "lm" & class(model)[1] != "clogit" & class(model)[1] != "coxph" & class(model)[1] != "geeglm" & class(model)[1] != "glmerMod")
     stop("Error: model must be either a glm, clogit, coxph, geeglm or glmerMod object")      

   if(class(model)[1] == "glm" & class(model)[2] == "lm"){
     theta1 <- as.numeric(model$coefficients[coef[1]])
     theta2 <- as.numeric(model$coefficients[coef[2]])
     theta3 <- as.numeric(model$coefficients[coef[3]])

     theta1.se <- summary(model)$coefficients[coef[1],2]
     theta2.se <- summary(model)$coefficients[coef[2],2]
     theta3.se <- summary(model)$coefficients[coef[3],2]
   }
   
   if(class(model)[1] == "clogit" | class(model)[1] == "coxph"){
     theta1 <- as.numeric(model$coefficients[coef[1]])
     theta2 <- as.numeric(model$coefficients[coef[2]])
     theta3 <- as.numeric(model$coefficients[coef[3]])
     
     theta1.se <- summary(model)$coefficients[coef[1],3]
     theta2.se <- summary(model)$coefficients[coef[2],3]
     theta3.se <- summary(model)$coefficients[coef[3],3]
   }

   if(class(model)[1] == "geeglm" & class(model)[2] == "gee"){
     theta1 <- as.numeric(model$coefficients[coef[1]])
     theta2 <- as.numeric(model$coefficients[coef[2]])
     theta3 <- as.numeric(model$coefficients[coef[3]])
     
     theta1.se <- summary(model)$coefficients[coef[1],2]
     theta2.se <- summary(model)$coefficients[coef[2],2]
     theta3.se <- summary(model)$coefficients[coef[3],2]
   }
   
   if(class(model)[1] == "glmerMod"){
     theta1 <- as.numeric(summary(model)$coefficients[coef[1]])
     theta2 <- as.numeric(summary(model)$coefficients[coef[2]])
     theta3 <- as.numeric(summary(model)$coefficients[coef[3]])
     
     theta1.se <- as.numeric(summary(model)$coefficients[coef[1],2])
     theta2.se <- as.numeric(summary(model)$coefficients[coef[2],2])
     theta3.se <- as.numeric(summary(model)$coefficients[coef[3],2])
   }
   
   if(theta1 < 0 | theta2 < 0)
     warning("At least one of the two regression coefficients is less than zero (i.e., OR < 1). Estimates of RERI and AP will be invalid. Estimate of SI valid.")
      
   if(param == "product"){
     
     # RERI:
     cov.mat <- vcov(model)
     
     h1 <- exp(theta1 + theta2 + theta3) - exp(theta1)
     h2 <- exp(theta1 + theta2 + theta3) - exp(theta2)
     h3 <- exp(theta1 + theta2 + theta3)
     
     reri.var <- (h1^2 * theta1.se^2) + 
       (h2^2 * theta2.se^2) + 
       (h3^2 * theta3.se^2) + 
       (2 * h1 * h2 * cov.mat[coef[1],coef[2]]) + 
       (2 * h1 * h3 * cov.mat[coef[1],coef[3]]) + 
       (2 * h2 * h3 * cov.mat[coef[2],coef[3]])
     reri.se <- sqrt(reri.var)
     
     reri.p <- exp(theta1 + theta2 + theta3) - exp(theta1) - exp(theta2) + 1
     reri.l <- reri.p - (z * reri.se)
     reri.u <- reri.p + (z * reri.se)
     reri <- data.frame(est = reri.p, lower = reri.l, upper = reri.u)
     
     
     # Multiplicative interaction:
     mult.p <- as.numeric(exp(theta3))
     
     mult.ci <- suppressMessages(confint(object = model, parm = coef[3]))
     mult.l <- as.numeric(exp(mult.ci[1]))
     mult.u <- as.numeric(exp(mult.ci[2]))

     multiplicative <- data.frame(est = mult.p, lower = mult.l, upper = mult.u)
     
     
     # APAB:
     cov.mat <- vcov(model)
     
     h1 <- ((exp(theta1 + theta2 + theta3) - exp(theta1)) / (exp(theta1 + theta2 + theta3))) - ((exp(theta1 + theta2 + theta3) - exp(theta1) - exp(theta2) + 1) / (exp(theta1 + theta2 + theta3)))
     h2 <- ((exp(theta1 + theta2 + theta3) - exp(theta2)) / (exp(theta1 + theta2 + theta3))) - ((exp(theta1 + theta2 + theta3) - exp(theta1) - exp(theta2) + 1) / (exp(theta1 + theta2 + theta3)))
     h3 <- 1 -((exp(theta1 + theta2 + theta3) - exp(theta1) - exp(theta2) + 1) / exp(theta1 + theta2 + theta3))
     
     apab.var <- (h1^2 * theta1.se^2) + 
       (h2^2 * theta2.se^2) + 
       (h3^2 * theta3.se^2) + 
       (2 * h1 * h2 * cov.mat[coef[1],coef[2]]) + 
       (2 * h1 * h3 * cov.mat[coef[1],coef[3]]) + 
       (2 * h2 * h3 * cov.mat[coef[2],coef[3]])
     apab.se <- sqrt(apab.var)
     
     apab.p <- (exp(theta1 + theta2 + theta3) - exp(theta1) - exp(theta2) + 1) / exp(theta1 + theta2 + theta3)
     
     apab.l <- apab.p - (z * apab.se)
     apab.u <- apab.p + (z * apab.se)
     apab <- data.frame(est = apab.p, lower = apab.l, upper = apab.u)
      
      
     # S:
     s.p <- (exp(theta1 + theta2 + theta3) - 1) / (exp(theta1) + exp(theta2) - 2)
     cov.mat <- vcov(model)
     
     # If model type is glm, cph or geeglm and point estimate of S is negative terminate analysis. Advise user to use a linear odds model:
     
     if(class(model)[1] == "glm" & class(model)[2] == "lm" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     if(class(model)[1] == "clogit" & class(model)[2] == "coxph" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     if(class(model)[1] == "geeglm" & class(model)[2] == "gee" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     if(class(model)[1] == "glmerMod" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     h1 <- ((exp(theta1 + theta2 + theta3)) / (exp(theta1 + theta2 + theta3) - 1)) - (exp(theta1) / (exp(theta1) + exp(theta2) - 2))
     h2 <- ((exp(theta1 + theta2 + theta3)) / (exp(theta1 + theta2 + theta3) - 1)) - (exp(theta2) / (exp(theta1) + exp(theta2) - 2))
     h3 <- exp(theta1 + theta2 + theta3) / (exp(theta1 + theta2 + theta3) - 1)
     
     lns.var <- h1^2 * theta1.se^2 + 
       h2^2 * theta2.se^2 + 
       h3^2 * theta3.se^2 + 
       (2 * h1 * h2 * cov.mat[coef[2],coef[1]]) + 
       (2 * h1 * h3 * cov.mat[coef[3],coef[1]]) + 
       (2 * h2 * h3 * cov.mat[coef[3],coef[2]])
     
     lns.se <- sqrt(lns.var)
     
     lns.p <- log(s.p)
     lns.l <- lns.p - (z * lns.se)
     lns.u <- lns.p + (z * lns.se)
     
     s.l <- exp(lns.l)
     s.u <- exp(lns.u)
     s <- data.frame(est = s.p, lower = s.l, upper = s.u)
     
     rval <- list(reri = reri, apab = apab, s = s, multiplicative = multiplicative)
       
     }
          
   if(param == "dummy"){
       
     # RERI:
     cov.mat <- vcov(model)
     
     h1 <- -exp(theta1)
     h2 <- -exp(theta2)
     h3 <-  exp(theta3)
     
     reri.var <- (h1^2 * (cov.mat[coef[1],coef[1]])) + 
       (h2^2 * (cov.mat[coef[2],coef[2]])) + 
       (h3^2 * (cov.mat[coef[3],coef[3]])) + 
       (2 * h1 * h2 * cov.mat[coef[1],coef[2]]) + 
       (2 * h1 * h3 * cov.mat[coef[1],coef[3]]) + 
       (2 * h2 * h3 * cov.mat[coef[2],coef[3]])
     reri.se <- sqrt(reri.var)
     
     reri.p <- exp(theta3) - exp(theta1) - exp(theta2) + 1
     reri.l <- reri.p - (z * reri.se)
     reri.u <- reri.p + (z * reri.se)
     
     reri <- data.frame(est = reri.p, lower = reri.l, upper = reri.u)
     
     
     # Multiplicative interaction:
     mult.p <- as.numeric(exp(theta3))
     
     mult.ci <- suppressMessages(confint(object = model, parm = coef[3]))
     mult.l <- as.numeric(exp(mult.ci[1]))
     mult.u <- as.numeric(exp(mult.ci[2]))
     
     multiplicative <- data.frame(est = mult.p, lower = mult.l, upper = mult.u)

     
     # APAB:
     cov.mat <- vcov(model)
     
     h1 <- -exp(theta1 - theta3)
     h2 <- -exp(theta2 - theta3)
     h3 <- (exp(theta1) + exp(theta2) - 1) / exp(theta3)
     
     apab.var <- (h1^2 * (cov.mat[coef[1],coef[1]])) + 
       (h2^2 * (cov.mat[coef[2],coef[2]])) + 
       (h3^2 * (cov.mat[coef[3],coef[3]])) + 
       (2 * h1 * h2 * cov.mat[coef[1],coef[2]]) + 
       (2 * h1 * h3 * cov.mat[coef[1],coef[3]]) + 
       (2 * h2 * h3 * cov.mat[coef[2],coef[3]])
     apab.se <- sqrt(apab.var)
     
     # apab.p <- exp(-theta3) - exp(theta1 - theta3) - exp(theta2 - theta3) + 1
     # Equation 4 (Skrondal 2003):
     apab.p <- (exp(theta3) - exp(theta1) - exp(theta2) + 1) / exp(theta3)
     apab.l <- apab.p - (z * apab.se)
     apab.u <- apab.p + (z * apab.se)
     apab <- data.frame(est = apab.p, lower = apab.l, upper = apab.u)
     
     
     # S:
     s.p <- (exp(theta3) - 1) / (exp(theta1) + exp(theta2) - 2)
     cov.mat <- vcov(model)
     
     # If model type is glm or cph and point estimate of S is negative terminate analysis.
     # Advise user to use a linear odds model:
     if(class(model)[1] == "glm" & class(model)[2] == "lm" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     if(class(model)[1] == "clogit" & class(model)[2] == "coxph" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     if(class(model)[1] == "geeglm" & class(model)[2] == "gee" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     if(class(model)[1] == "glmerMod" & s.p < 0){
       warning(paste("Point estimate of synergy index (S) is less than zero (", round(s.p, digits = 2), ").\n  Confidence intervals cannot be calculated using the delta method. Consider re-parameterising as linear odds model.", sep = ""))
     }
     
     # Use delta method (Hosmer and Lemeshow 1992) if model type is glm, clogit or cph:
     h1 <- -exp(theta1) / (exp(theta1) + exp(theta2) - 2)
     h2 <- -exp(theta2) / (exp(theta1) + exp(theta2) - 2)
     h3 <- exp(theta3) / (exp(theta3) - 1)
     
     lns.var <- h1^2 * theta1.se^2 + 
       h2^2 * theta2.se^2 + 
       h3^2 * theta3.se^2 + 
       (2 * h1 * h2 * cov.mat[coef[2],coef[1]]) + 
       (2 * h1 * h3 * cov.mat[coef[3],coef[1]]) + 
       (2 * h2 * h3 * cov.mat[coef[3],coef[2]])
     
     lns.se <- sqrt(lns.var)
     
     lns.p <- log(s.p)
     lns.l <- lns.p - (z * lns.se)
     lns.u <- lns.p + (z * lns.se)
     
     s.l <- exp(lns.l)
     s.u <- exp(lns.u)
     s <- data.frame(est = s.p, lower = s.l, upper = s.u)
     
     rval <- list(reri = reri, apab = apab, s = s, multiplicative = multiplicative)
     }

   return(rval)
}