R/path_specific_onesimulation.R

In causalPAF: Causal Effect for Population Attributable Fractions (PAF)

#' @title Path specific population attributable fraction for a mediating pathway.
#' @description Path specific population attributable fraction for a mediating pathway.
#' This metric has several desirable properties. This is a kind of sequential PAF on pathways,
#' but now corresponding to eliminating the mediating pathway first.
#' \deqn{PAF_{A->M->Y} = (P(Y=1)-P(Y_{A,M_{0}}=1)/P(Y=1) }
#' \if{html}{\figure{CausalDAG.jpg} {options: width=100 alt="Causal Bayesian network DAG"} }
#' \if{latex}{\figure{CausalDAG.jpg}{options: width=1in}}
#' @param data_frame Data frame containing the data. The data frame has cases in rows and variables in columns.
#' @param exposure The exposure name in the form of character string e.g. "phys" for physical exercise.
#' @param mediator The mediator name in the form of character string e.g. "whr" for waist hip ratio.
#' @param response The outcome name in the form of character string e.g. "case" for a stroke case.
#' @param mediator_model A list containing each of the fitted mediator regression models e.g.
#' \code{mediator_model=list(model_list[[6]],model_list[[7]],model_list[[8]])}.
#' @param response_model is a regression for the outcome on all mediators together with all parents and confounders of the mediators in a Markov causal Bayesian network DAG. A regression model fitted for the response in a causal Bayesian network excluding ``children'' of the mediators in the causal Bayesian network. See example in tutorial.This model can be listed either as (1) an empty list ( response_model_mediators = list() ) or (2) the user can specify their own customised causal regression model(s) to use. When it is listed as an empty list the 'causalPAF' package will fit the response_model_mediators regression model automatically based on the causal DAG supplied by the user in in_outArg. Alternatively, the user can specify the exact model(s) that the user wishes to use, these model(s) must be in list format (list() where length(response_model_mediators) == length(mediator) ), the same length as the parameter, mediator, with the user customised model for each mediator listed in the same order as in the parameter, mediator, and if there is only one model, it must be listed each time within the list() so that length(response_model_mediators) == length(mediator).
#' @param response_model_2 A regression model fitted for the response in a causal Bayesian network excluding ``children'' of the exposure in the causal Bayesian network. This regression model will not adjust for mediators (exclude mediators) of the exposure in the regression model so that the total effect of the exposure on the response can be modelled. This model can be listed either as (1) an empty list ( response_model_exposure = list() ) or (2) the user can specify their own customised causal regression model to use. If specified as an empty list, list(), then the 'causalPAF' package will define and fit the model automatically based on the causal DAG defined by the in_outArg parameter. Alternatively, the user can specify the exact model that the user wishes to use, this model must be in list format (list() where length(response_model_exposure) == 1 ), of length 1, assuming only one exposure of interest (other exposures can be risk factors) and the model must be defined within a list() since the package assumes a list() format is supplied. See example in tutorial. E.G. If physical exercise ("exer") in the example given in the diagram is the exposure. Then the regression would include all parents of "exer" (i.e. sex, region, educ, age) as well as risk factors at the same level of the causal Bayesian network (i.e. stress, smoke, diet, alcoh).
#' @param weights  A numeric n x \eqn{1} vector where n is the number of patients in the case control data frame.
#' Different weighting approaches can be applied as per the literature, Pathway-specific population attributable
#' fractions (PS-PAF) O'Connell and Ferguson (\eqn{2022}), Revisiting sequential population attributable fractions Ferguson, O'Connell and O'Donnell (\eqn{2020}).
#' For more information on weighting, a tutorial paper will be published and linked here when it is published.
#' For example, O'Connell and Ferguson (\eqn{2022}), for a case-control study design, when prevalence pi is known, and
#' the sampled disease cases and controls are randomly selected from their respective populations. We assume for
#' simplicity that the case to control matching ration is \eqn{1} to r, for some r greater than or equal to \eqn{1}.
#' Under assumptions outlined in O'Connell and Ferguson (\eqn{2022}), the components of the PAF can be found as the
#' corresponding empirical expectations and distributions in the re-weighted dataset where cases are assigned the weights,
#' \eqn{1}, and controls are assigned the weights of (((\eqn{1} divided by pi) minus \eqn{1}) all divided by r). Effectively, then we can think of the re-weighted population as a random
#' sample. A tutorial paper will be linked here when published for more information which will show how to apply the weighting for different study designs.
#' @export
#' @import stats
# #' @importFrom stats formula predict filter lag
#' @keywords internal
#' @return \item{path_specific_PAF }{path specific PAF}

path_specific_onesimulation <- function(data_frame, exposure, mediator, response, mediator_model, response_model, response_model_2, weights ){

  mediator_outcomes <- c()
  for(i in 1:length(mediator_model)) mediator_outcomes[i] <-  as.character(formula(mediator_model[[i]]))[2]
  which.model <- grep(paste('^',mediator,'$',sep=''),mediator_outcomes,perl=TRUE)
  data_frame_pathspecific=data_frame
  data_frame_pathspecific[,grep(paste('^',exposure,'$',sep=''),colnames(data_frame_pathspecific),perl=TRUE)] <- levels(data_frame_pathspecific[,grep(paste('^',exposure,'$',sep=''),colnames(data_frame),perl=TRUE)])[1]
  for(i in which.model){

    ### simulate mediators given exposure at reference (except for mediator in question)
    thecol <- grep(paste('^',mediator_outcomes[i],'$',sep=''),colnames(data_frame_pathspecific),perl=TRUE)
    data_frame_pathspecific <- do_sim(mediator_model[[i]],data_frame_pathspecific)

  }
  #  make sure exposure is set at natural value
  data_frame_pathspecific[,grep(paste('^',exposure,'$',sep=''),colnames(data_frame_pathspecific),perl=TRUE)] <- data_frame[,grep(paste('^',exposure,'$',sep=''),colnames(data_frame_pathspecific),perl=TRUE)]
  predicted_response <- predict(response_model,newdata=data_frame_pathspecific,type="response")

  path_specific_PAF <- sum(weights*(predict(response_model,type="response")-predicted_response))/sum(weights*predict(response_model,type="response"))

  return(path_specific_PAF)
}