R/predict.sim_dlim.R

In dlim: Distributed Lag Interaction Model

#' Simulated DLIM Predictions
#' @description This function estimates cumulative and non-cumulative lag/modifier coefficients from a model in which the response is regressed on a cross-basis generated by the \code{cross_basis()} function.
#' @seealso \link[dlim]{predict.dlim}
#' @seealso Type \code{vignette('dlimOverview')} for a detailed description.
#' @export
#' @import dlnm 
#' @param object an object of class "\code{dlim}" 
#' @param newdata vector of modifiers for inference (class "\code{numeric}")
#' @param type Type of prediction. "response" for predicted responses, "DLF" for the estimated distributed lag functions, "CE" for cumulative effects (class "\code{character}")
#' @param ... additional arguments affecting the predictions produced
#' @return This function returns a list of 4 or 7 elements:
#' \item{est_dlim}{\code{est_dlim} element from \code{predict.dlim} (class "\code{list}")}
#' \item{cb}{cross-bais from \code{object} (class "\code{cross-basis}")}
#' \item{fit}{\code{fit} from \code{object} (class "\code{lm}", "\code{glm}", "\code{gam}")}
#' \item{true_betas}{\code{true_betas} from \code{object} (class "\code{matrix}")}
#' \item{cb_dlm}{\code{cb_dlm} from \code{object} (class "\code{crosspred}")}
#' \item{model_dlm}{\code{model_dlm} from \code{object} (class "\code{lm}", "\code{glm}", "\code{gam}")}
#' \item{est_dlm}{cumulative and/or point-wise estimates, standard errors, and confidence intervals for the DLM (class "\code{list}")}


predict.sim_dlim <- function(object, newdata=NULL, type=c("DLF","CE", "response"), ...){

  if(!inherits(object, "sim_dlim")){
    stop("Object not of class sim_dlim")
  }


  #DLIM
  class(object) <- "dlim"
  dlim_pred <- predict(object = object, newdata = newdata, type=type)

  result <- c()
  result$est_dlim <- dlim_pred$est_dlim
  result$cb <- object$cb
  result$fit <- object$fit
  result$true_betas <- object$true_betas

  #DLM

  if(!is.null(object$cb_dlm)){
    cb_dlm <- object$cb_dlm
    model_dlm <- object$model_dlm
    cb <- object$cb
    m <- length(object$modifiers)

    #get the crosspred object using dlnm package
    dlm_crosspred <- crosspred(cb_dlm,model_dlm,at=rep(1,object$cb$L+1),cen = FALSE)    #beta estimates for each lag (are different for each lag)
    beta_by_lag <- dlm_crosspred$matfit
    #cumulative beta estimate
    cumul_betas <- dlm_crosspred$allfit
    #estimated non-cumulative betas (m X lag matrix, each row is the same)
    betas <- matrix(rep(beta_by_lag,m),nrow=m,byrow=TRUE)
    #cumulative beta estimates (length #observations/modifiers, will all be the same value)
    betas_cumul <- rep(cumul_betas,m)

    est_dlm <- list()
    est_dlm$betas <- betas
    est_dlm$betas_cumul <- betas_cumul

    #Cumulative Confidence intervals
    cumul_LB <- rep(dlm_crosspred$alllow,m)
    cumul_UB <- rep(dlm_crosspred$allhigh,m)
    cumul_SE <- rep(dlm_crosspred$allse,m)

    est_dlm$cumul_LB <- cumul_LB
    est_dlm$cumul_UB <- cumul_UB
    est_dlm$cumul_SE <- cumul_SE

    #DLF Confidence intervals
    LB <- matrix(rep(dlm_crosspred$matlow,m),nrow=m,byrow=TRUE)
    UB <- matrix(rep(dlm_crosspred$mathigh,m),nrow=m,byrow=TRUE)
    SE <- matrix(rep(dlm_crosspred$matse,m),nrow=m,byrow=TRUE)

    est_dlm$LB <- LB
    est_dlm$UB <- UB
    est_dlm$SE <- SE

    #save DLM results
    result$cb_dlm <- object$cb_dlm
    result$model_dlm <- object$model_dlm
    result$est_dlm <- est_dlm

  }

  return(result)
}