R/AFivglm.R

In AF: Model-Based Estimation of Confounder-Adjusted Attributable Fractions

############## AF function for a ivglm object #####################
#' @title Attributable fraction function based on Instrumental Variables (IV) regression as an \code{\link[ivtools]{ivglm}} object in the \code{ivtools} package.
#' @description \code{AFivglm} estimates the model-based adjusted attributable fraction from a Instrumental Variable regression from a \code{\link[ivtools]{ivglm}} object. The IV regression can be estimated by either G-estimation or Two Stage estimation for a binary exposure and outcome.
#' @param object a fitted Instrumental Variable regression of class "\code{\link[ivtools]{ivglm}}".
#' @param data an optional data frame, list or environment (or object coercible by \code{as.data.frame} to a data frame) containing the variables in the model. If not found in data, the variables are taken from environment (\code{formula}), typically the environment from which the function is called.
#' @return \item{AF.est}{estimated attributable fraction.}
#' @return \item{AF.var}{estimated variance of \code{AF.est}. The variance is obtained by combining the delta methods with the sandwich formula.}
#' @details \code{AFivglm} estimates the attributable fraction for an IV regression
#' under the hypothetical scenario where a binary exposure \code{X} is eliminated from the population.
#' The estimate can be adjusted for IV-outcome confounders \code{L} in the \code{\link[ivtools]{ivglm}} function. 
#' Let the AF function be defined as
#' \deqn{AF=1-\frac{Pr(Y_0=1)}{Pr(Y=1)}}{AF = 1 - \frac{Pr(Y0=1)}{Pr(Y=1)}}
#' where \eqn{Pr(Y_0=1)}{Pr(Y0=1)} denotes the counterfactual outcome prevalence had everyone been unexposed and \eqn{Pr(Y=1)}{Pr(Y=1)} denotes the factual outcome prevalence.
#' If the instrument \code{Z} is valid, conditional on covariates \code{L}, i.e. fulfills the IV assumptions 1) the IV should have a (preferably strong) association with
#'the exposure, 2) the effect of the IV on the outcome should only go through the exposure and 3) the IV-outcome association should be unconfounded
#' (Imbens and Angrist, 1994) then \eqn{Pr(Y_0=1)}{Pr(Y0=1)} can be estimated.
#' @author Elisabeth Dahlqwist, Arvid \enc{Sjölander}{Sjolander}
#' @seealso \code{\link[ivtools]{ivglm}} used for fitting the causal risk ratio or odds ratio using the G-estimator or Two stage estimator.
#' @references Dahlqwist E., Kutalik Z., \enc{Sjölander}{Sjolander}, A. (2019). Using Instrumental Variables to estimate the attributable fraction. \emph{Manuscript}.
#' @examples
#' # Example 1
#' set.seed(2)
#' n <- 5000
#' ## parameter a0 determines the outcome prevalence
#' a0 <- -4
#' psi.true <- 1
#' l <- rbinom(n, 1, 0.5)
#' u <- rbinom(n, 1, 0.5)
#' z <- rbinom(n, 1, plogis(a0))
#' x <- rbinom(n, 1, plogis(a0+3*z+ u))
#' y <- rbinom(n, 1, exp(a0+psi.true*x+u))
#' d <- data.frame(z,u,x,y,l)
#' ## Outcome prevalence 
#' mean(d$y)
#' 
#' ####### G-estimation
#' ## log CRR
#' fitz.l <- glm(z~1, family=binomial, data=d)
#' gest_log <- ivglm(estmethod="g", X="x", Y="y",
#'                   fitZ.L=fitz.l, data=d, link="log")
#' AFgestlog <- AFivglm(gest_log, data=d)
#' summary(AFgestlog)
#' 
#' ## log COR
#' ## Associational model, saturated
#' fit_y <- glm(y~x+z+x*z, family="binomial", data=d)
#' ## Estimations of COR and AF
#' gest_logit <- ivglm(estmethod="g", X="x", Y="y",
#'                     fitZ.L=fitz.l, fitY.LZX=fit_y,
#'                     data=d, link="logit")
#' AFgestlogit <- AFivglm(gest_logit, data = d)
#' summary(AFgestlogit)
#' 
#' ####### TS estimation
#' ## log CRR
#' # First stage
#' fitx <- glm(x ~ z, family=binomial, data=d)
#' # Second stage
#' fity <- glm(y ~ x, family=poisson, data=d)
#' ## Estimations of CRR and AF
#' TSlog <- ivglm(estmethod="ts", X="x", Y="y",
#'                fitY.LX=fity, fitX.LZ=fitx, data=d, link="log")
#' AFtslog <- AFivglm(TSlog, data=d)
#' summary(AFtslog)
#' 
#' ## log COR
#' # First stage
#' fitx_logit <- glm(x ~ z, family=binomial, data=d)
#' # Second stage
#' fity_logit <- glm(y ~ x, family=binomial, data=d)
#' ## Estimations of COR and AF
#' TSlogit <- ivglm(estmethod="ts", X="x", Y="y",
#'                  fitY.LX=fity_logit, fitX.LZ=fitx_logit,
#'                   data=d, link="logit")
#' AFtslogit <- AFivglm(TSlogit, data=d)
#' summary(AFtslogit)
#'
#' ## Example 2: IV-outcome confounding by L
#' ####### G-estimation
#' ## log CRR
#' fitz.l <- glm(z~l, family=binomial, data=d)
#' gest_log <- ivglm(estmethod="g", X="x", Y="y",
#'                   fitZ.L=fitz.l, data=d, link="log")
#' AFgestlog <- AFivglm(gest_log, data=d)
#' summary(AFgestlog)
#' 
#' ## log COR
#' ## Associational model
#' fit_y <- glm(y~x+z+l+x*z+x*l+z*l, family="binomial", data=d)
#' ## Estimations of COR and AF
#' gest_logit <- ivglm(estmethod="g", X="x", Y="y",
#'                     fitZ.L=fitz.l, fitY.LZX=fit_y,
#'                     data=d, link="logit")
#' AFgestlogit <- AFivglm(gest_logit, data = d)
#' summary(AFgestlogit)
#' 
#' ####### TS estimation
#' ## log CRR
#' # First stage
#' fitx <- glm(x ~ z+l, family=binomial, data=d)
#' # Second stage
#' fity <- glm(y ~ x+l, family=poisson, data=d)
#' ## Estimations of CRR and AF
#' TSlog <- ivglm(estmethod="ts", X="x", Y="y",
#'                fitY.LX=fity, fitX.LZ=fitx, data=d,
#'                link="log")
#' AFtslog <- AFivglm(TSlog, data=d)
#' summary(AFtslog)
#' 
#' ## log COR
#' # First stage
#' fitx_logit <- glm(x ~ z+l, family=binomial, data=d)
#' # Second stage
#' fity_logit <- glm(y ~ x+l, family=binomial, data=d)
#' ## Estimations of COR and AF
#' TSlogit <- ivglm(estmethod="ts", X="x", Y="y",
#'                  fitY.LX=fity_logit, fitX.LZ=fitx_logit,
#'                  data=d, link="logit")
#' AFtslogit <- AFivglm(TSlogit, data=d)
#' summary(AFtslogit)
#' @importFrom stats plogis
#' @import ivtools
#' @export
AFivglm<- function(object, data){
  call <- match.call()
  method <- object$input$estmethod 
  # Warning if the object is not a ivglm object
  if(!(as.character(object$call[1]) == "ivglm"))
    stop("The object is not an ivglm object", call. = FALSE)
  #######################
    if(method == "g"){
      result <- AFgest(object = object, data = data)
    }
    if(method == "ts"){
      result <- AFtsest(object = object, data = data)
    }
    #### Output
    out <- result
  
  return(out)
}