R/msm.r

In CICI: Causal Inference with Continuous (Multiple Time Point) Interventions

msm <- function(
    X, formula, family = gaussian, se = NULL, cilevel = 0.95, abar= NULL
) {
  # --------------------------- Check input arguments --------------------------
  if ("simulated.data" %in% names(X)) {
    if (is.null(X[["simulated.data"]])) {
      stop("Set 'ret=TRUE' in gformula().")
    }
  } else {
    stop("Set 'X' to a 'gformula' object resulting from 'gformula', not 'sgf'.")
  }
  
  # Survival setting: transform matrices to dataframes
  if (is.matrix(X$simulated.data[[1]][[1]])) {
    X$simulated.data <- lapply(X$simulated.data, function(blevel) {
      lapply(blevel, function(m) as.data.frame(m))
    })
  }
  
  if (is.null(se)) {
    if (length(X$simulated.data) > 1) {
      se <- "bootstrap"
    } else {
      se <- "glm"
      warning("GLM CIs for MSM parameters suffer from undercoverage. Consider
      setting 'B > 0' in gformula() to use bootstrap CIs.")
    }
  } else {
    if (!(se %in% c("glm", "bootstrap"))) {
      stop("Set 'se' to either 'glm' or 'bootstrap' standard errors.")
    } else if (se == "bootstrap" && length(X$simulated.data) == 1) {
      stop("Set 'B > 0' in gformula() to use se = 'bootstrap'.")
    }
  }
  
  if (!(is.numeric(cilevel) && length(cilevel) == 1 &&
        cilevel > 0 && cilevel < 1)) {
    stop("'cilevel' must be a single numeric value between 0 and 1.")
  }
  
  if (is.character(formula)) {
    formula <- as.formula(formula)
  } else if (is.language(formula) && !inherits(formula, "formula")) {
    formula <- eval(formula)
  } else if (!inherits(formula, "formula")) {
    stop("'formula' must be a formula, quoted expression, or character string.")
  }
  vars <- all.vars(formula)
  varsNA <- setdiff(vars, names(X$simulated.data[[1]][[1]]))
  if (length(varsNA) > 0) {
    stop(sprintf(
      "The following variables used in 'formula' are not present in the data:
       %s.\n", paste(varsNA, collapse = ", ")
    ))
  }
  
  if(is.vector(X$setup$abar)){if(length(X$setup$abar)==1){stop("An MSM needs at least two interventions")}}
  if(is.matrix(X$setup$abar)){if(nrow(X$setup$abar)==1){stop("An MSM needs at least two interventions")}}
  if(is.list(X$setup$abar) & (is.null(abar)==FALSE)){stop("Using 'abar' with individual interventions not allowed")}
  if(is.matrix(X$setup$abar) & is.vector(abar)){stop("Chosen 'abar' is of wrong type.")}
  if(is.vector(X$setup$abar) & is.matrix(abar)){stop("Chosen 'abar' is of wrong type.")}
  
  # ----------------------------- Evaluate the MSM -----------------------------
  # Stack the simulated data
  if (X$setup$i.type != "natural") {
    df <- do.call(rbind, X$simulated.data[[1]])
  } else {
    stop("MSMs for specified 'abar' (natural intervention) type not possible.")
  }
  
  # use subset if specified
  if(is.null(abar)==FALSE){
    if(is.vector(abar)){abar<-matrix(rep(abar,length(X$setup$Anodes)),ncol=length(X$setup$Anodes))}
    colnames(abar) <- X$setup$Anodes
    key1 <- interaction(df[,X$setup$Anodes], drop = TRUE)
    key2 <- interaction(as.data.frame(abar), drop = TRUE)
    is_in <- key1 %in% key2
    if(any(key2%in%key1==FALSE)){stop(" (Some) Selected 'abar' do not appear in data")}
    df <- df[is_in,]
  }
  
  # Fit the MSM
  fit <- glm(formula, data = df, family = family)
  coefs <- coef(fit)
  
  # Calculate standard errors and confidence intervals
  if (se == "glm") {
    cibounds <- suppressMessages(confint(fit, level = cilevel, trace=F))
  } else if (se == "bootstrap") {
    B <- length(X$simulated.data)
    p <- length(coef(fit))       
    boots <- matrix(NA_real_, nrow = B-1, ncol = p)
    for (b in 2:B) {
      # Stack the simulated data for the current bootstrap sample
      df <- do.call(rbind, X$simulated.data[[b]])
      # Fit the MSM
      bootfit <- glm(formula, data = df, family = family)
      # Extract MSM parameters
      boots[b-1, ] <- coef(bootfit)
    }
    # Bootstrap CI
    cibounds <- t(apply(
      boots, 2, quantile,
      probs = c((1 - cilevel) / 2, 1 - (1 - cilevel) / 2)
    ))
  }
  
  if(length(X$simulated.data)>1 & se=="bootstrap"){vcov<-cov(boots)}else{vcov<-NULL}
  
  # ------------------------------- Return Object ------------------------------
  returnObj <- list(
    "MSM" = fit,
    "coefs" = coefs,
    "CIlow" = cibounds[, 1],
    "CIup" = cibounds[, 2],
    "formula" = formula,
    "se" = se,
    "vcov" = vcov
  )
  
  class(returnObj) <- "msmResult"
  return(returnObj)
}