Nothing
#' Fit DLIM for simulation
#' @description Fit DLIM for simulation
#' @seealso \link[dlim]{dlim}
#' @seealso \link[dlim]{sim_data}
#' @seealso Type \code{vignette('dlimOverview')} for a detailed description.
#' @export
#' @import mgcv
#' @import dlnm
#' @param data output from \code{sim_data}
#' @param df_m degrees of freedom for modifiers
#' @param df_l degrees of freedom for lags
#' @param penalize True to penalize model
#' @param pen_fn if penalizing, can specify "ps" for penalized B-splines or "cr" for cubic regression splines with penalties on second derivatives
#' @param mod_args a list of additional arguments for the spline function (must be named by argument)
#' @param lag_args a list of additional arguments for the spline function (must be named by argument)
#' @param fit_dlm True to additionally fit dlm for comparison
#' @param model_type "linear" for a DLIM with linear interaction, "quadratic" for a DLIM with quadratic interaction, "standard" for a DLIM with splines
#' @param ... arguments to pass to model fitting function
#' @return This function returns an object of class "\code{sim_dlim}"
#' \item{cb}{DLIM cross-basis (class "\code{cross-basis}")}
#' \item{fit}{DLIM model fit (class "\code{lm}", "\code{glm}", "\code{gam}")}
#' \item{cb_dlm}{DLM cross-basis (class "\code{crossbasis}")}
#' \item{model_dlm}{DLM model fit (class "\code{lm}", "\code{glm}", "\code{gam}")}
#' \item{true_betas}{true linear effect of the exposure on the response for each individual and time point (class "\code{matrix}")}
#' \item{modifiers}{\code{modifiers} from \code{numeric}}
#' \item{data}{\code{data} (class "\code{list}")}
sim_dlim <- function(data, df_m, df_l, penalize=TRUE, pen_fn = "ps", mod_args=NULL, lag_args=NULL, fit_dlm=FALSE, model_type="standard",...){
#fit DLIM
model <- dlim(y = data$y,
x = data$x,
modifiers = data$modifiers,
z = data$Z,
df_m = df_m,
df_l = df_l,
penalize=penalize,
model_type=model_type,
mod_args = mod_args,
lag_args = lag_args,
...)
#fit DLM
if(fit_dlm){
#set up
modifiers <- matrix(data$modifiers, ncol=1)
if(!is.null(data$Z)){
design1 <- data.frame(intercept=rep(1,length(data$modifiers)), data$modifiers,data$Z)
}else{
design1 <- data.frame(intercept=rep(1,length(data$modifiers)), data$modifiers)
}
y <- data$y
#cross-basis
if(penalize){
cb_dlm <- crossbasis(x=data$x,argvar=list(fun="lin"),arglag = list(fun=pen_fn,df=df_l))
}else{
cb_dlm <- crossbasis(x=data$x,argvar=list(fun="lin"),arglag = list(fun="ns",df=df_l))
}
#fit model
Z <- model.matrix(~ 0+.,model.frame(~ ., design1, na.action=na.pass)) #handles factor covariates and missing values
if(penalize){
penalty <- cbPen(cb_dlm)
model_dlm <- do.call("gam",list(formula=0+y~cb_dlm+Z,paraPen = list(cb_dlm = penalty), method = "REML"))
}else{
model_dlm <- do.call("gam",list(formula=y~0+cb_dlm+Z))
}
}
if(fit_dlm){
results <- list("cb" = model$cb, "fit" = model$fit, "cb_dlm" = cb_dlm, "model_dlm" = model_dlm, "true_betas" = data$betas, "modifiers" = data$modifiers, data = data)
}else{
results <- list("cb" = model$cb, "fit" = model$fit, "true_betas" = data$betas, "modifiers" = data$modifiers, data = data)
}
attr(results, "model_type") <- model_type
class(results) <- "sim_dlim"
return(results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.