mlVAR0 <- function(
data, # Data frame
vars, # Vector of variables to include in analysis
idvar, # String indicating the subject id variable name
lags = 1, # Vector indicating the lags to include
dayvar, # string indicating the measurement id variable name (if missing, every measurement is set to one day)
beepvar, # String indicating beep per day (is missing, is added)
periodvar, # string indicating the period of measurement.
treatmentvar, # character vector indicating treatment
covariates, # character indicating covariates independent of measurement.
timevar,
maxTimeDiff, # If not missing. maximum time difference.
control = list(optimizer = "bobyqa"), # "bobyqa" or "Nelder_Mead"
# windowSize, # Assign to use moving window estimation. If missing, does not use moving window estimation
verbose = TRUE, # Include progress bar?
# standardize = c("inSubject","none","general"),
orthogonal, # by default set to TRUE when nvar > 6 and FALSE otherwise
estimator = c("lmer","lmmlasso"),
method = c("default","stepwise","movingWindow"),
laginteractions = c("none","mains","interactions"), # Include interactions with lag?
critFun = BIC,
lambda = 0,
center = c("inSubject","general","none")
)
{
warning("Function is deprecated. Use mlVAR instead.")
# standardize = c("inSubject","none","general")
center <- match.arg(center)
if (is(data,"mlVARsim0")){
vars <- data$vars
idvar <- data$idvar
data <- data$Data
}
method <- match.arg(method)
estimator <- match.arg(estimator)
laginteractions <- match.arg(laginteractions)
if (missing(orthogonal)){
if (length(vars) > 6 & method != "movingWindow"){
if (verbose) message("More than 6 nodes and method != 'movingWindow', correlations between random effects are set to zero (orthogonal = TRUE)")
orthogonal <- TRUE
} else {
orthogonal <- FALSE
}
}
# Check input:
stopifnot(!missing(vars))
stopifnot(!missing(idvar))
# inout list (to include in output):
input <- list(vars = vars, lags = lags)
# Add day id if missing:
if (missing(idvar))
{
idvar <- "ID"
data[[idvar]] <- 1
} else input$idvar <- idvar
# Add day var if missing:
if (missing(dayvar))
{
dayvar <- "DAY"
data[[dayvar]] <- 1
} else input$dayvar <- dayvar
# Add beep var if missing:
if (missing(beepvar))
{
beepvar <- "BEEP"
data[[beepvar]] <- ave(data[[idvar]],data[[idvar]],data[[dayvar]],FUN = seq_along)
} else input$beepvar <- beepvar
# Add period var if missing:
if (missing(periodvar))
{
periodvar <- "PERIOD"
data[[periodvar]] <- 1
} else input$periodvar <- periodvar
if (!missing(treatmentvar)) input$treatmentvar <- treatmentvar
if (!missing(covariates)) input$covariates <- covariates
# Remove NA period, day or beeps:
data <- data[!is.na(data[[idvar]]) & !is.na(data[[periodvar]]) & !is.na(data[[dayvar]]) & !is.na(data[[beepvar]]), ]
# Create augmented lagged data:
augData <- plyr::ddply(data, c(idvar,dayvar,periodvar), function(x){
# Check for duplicate beeps:
if (any(duplicated(x[[beepvar]])))
{
stop("Duplicated beepnumbers found. ")
}
# Order by beep:
x <- x[order(x[[beepvar]]),]
# Augment missing:
beepseq <- seq(1, max(x[[beepvar]]))
if (!identical(x[[beepvar]], beepseq))
{
dummy <- data.frame(ID = unique(x[[idvar]]), PERIOD= unique(x[[periodvar]]), DAY = unique(x[[dayvar]]), BEEP = beepseq)
names(dummy) <- c(idvar,periodvar,dayvar,beepvar)
x <- plyr::join(dummy,x, by = names(dummy))
}
# Lag variables:
for (l in lags)
{
if (l > nrow(x)) stop("Lag is larger than number of measurements")
lagDF <- x[-(nrow(x)-(1:l)+1),vars]
lagDF <- lagDF[c(rep(NA,l),seq_len(nrow(lagDF))),]
names(lagDF) <- paste0("L",l,"_",names(lagDF))
x <- cbind(x,lagDF)
}
return(x)
})
### MULTILEVEL ANALYSIS ###
# All lagged variables:
AllLagVars <- c(sapply(lags, function(l) paste0("L",l,"_",vars)))
# Lagged interactions:
# Lag interactions:
if (laginteractions != "none" || !missing(maxTimeDiff)){
if (missing(timevar)) stop("'timevar' needed to include time difference")
call <- substitute(augData %>% group_by(.data[['idvar']]) %>% mutate(LAGDIFF = as.numeric(c(NA,difftime(x[-1],x[-length(x)],units="hours")))),
list(x = as.name(timevar)))
augData <- eval(call)
}
if (laginteractions != "none"){
AllLagVars <- c(AllLagVars,"LAGDIFF")
if (laginteractions == "interactions"){
for (it in seq_along(vars)){
augData[[paste0("lag_x_",vars[it])]] <- augData[[vars[it]]] * augData$LAGDIFF
}
AllLagVars <- c( AllLagVars ,paste0("lag_x_",vars))
}
}
# Center predictors:
# standardize <- match.arg(standardize)
Center <- function(x) return(x - mean(x,na.rm=TRUE))
if (center =="inSubject"){
for(i in unique(data[[idvar]])){
for (var in AllLagVars){
augData[augData[[idvar]]==i,var] <- Center(augData[augData[[idvar]]==i,var])
}
}
} else if (center == "general"){
augData[,AllLagVars] <- sapply(augData[,AllLagVars] ,Center)
}
if ( !missing(maxTimeDiff)){
augData <- augData %>% filter(.data[['LAGDIFF']] < maxTimeDiff)
}
### RUN NODEWISE ANALYSES ###
NodeWise_Results <- list()
Results <- list()
formulas <- list()
FixEf <- list()
FixEf_SE <- list()
# Run models:
if (verbose){
pb <- txtProgressBar(min=0,max=length(vars),style=3)
}
for (j in seq_along(vars)){
whichAuto <- grepl(paste0("^L1_",vars[j],"$"), AllLagVars)
if (method == "stepwise"){
NodeWise_Results[[j]] <- Stepwise(
aData = augData,
response = vars[j],
idvar = idvar, # String indicating the subject id variable name
autoLaggedVars = AllLagVars[whichAuto],
laggedVars = AllLagVars[!whichAuto], # Vector indicating the lags to include
periodvar = periodvar, # string indicating the period of measurement.
treatmentvar = treatmentvar, # character vector indicating treatment
covariates = covariates, # character indicating covariates independent of measurement.
control = control, # "bobyqa" or "Nelder_Mead"
timevar = timevar,
critFun=critFun,
progress=verbose,
estimator = estimator,
lambda = lambda,
orthogonal = orthogonal
)
} else if (method == "default") {
NodeWise_Results[[j]] <- NodeWise(
aData = augData,
response = vars[j],
idvar = idvar, # String indicating the subject id variable name
autoLaggedVars = AllLagVars[whichAuto],
laggedVars = AllLagVars[!whichAuto], # Vector indicating the lags to include
periodvar = periodvar, # string indicating the period of measurement.
treatmentvar = treatmentvar, # character vector indicating treatment
covariates = covariates, # character indicating covariates independent of measurement.
control = control, # "bobyqa" or "Nelder_Mead"
timevar = timevar,
estimator = estimator,
lambda = lambda,
orthogonal = orthogonal
)
} else if (method=="movingWindow") {
NodeWise_Results[[j]] <- movingWindow(
aData = augData,
response = vars[j],
idvar = idvar, # String indicating the subject id variable name
autoLaggedVars = AllLagVars[whichAuto],
laggedVars = AllLagVars[!whichAuto], # Vector indicating the lags to include
periodvar = periodvar, # string indicating the period of measurement.
treatmentvar = treatmentvar, # character vector indicating treatment
covariates = covariates, # character indicating covariates independent of measurement.
control = control, # "bobyqa" or "Nelder_Mead"
timevar = timevar,
estimator = estimator,
lambda = lambda,
orthogonal = orthogonal
)
} else stop("Method not implemented")
if (verbose){
setTxtProgressBar(pb, j)
}
}
if (verbose){
close(pb)
}
# Fixed effects for each lag:
# ### MOVE TO MOVING WINDOW FUNCTION ####
#
# # Construct the window matrix:
# Neffect <- length(AllLagVars) - 1
#
#
# if (!missing(windowSize)){
# # Number of possible effects to include is amount of lagged - 1 (autocor is always included):
#
#
# # Window matrix is Neffect * windowSize
# # Randomize:
# if (Neffect>2) Rand <- sample(seq_len(Neffect)) else Rand <- 1
# WindowMatrix <- matrix(Rand[(0:(windowSize-1) + rep(0:(Neffect-1),each=windowSize)) %% Neffect + 1],
# Neffect, windowSize, byrow=TRUE)
#
# } else WindowMatrix <- matrix(1:Neffect,1)
#
# cur <- 1
# for (j in seq_along(vars)){
#
# for (w in seq_len(nrow(WindowMatrix))){
# PredCur <- Pred
#
# # Always include the autoregression
# whichAuto <- grep(paste0("^L1_",vars[j]), AllLagVars)
# RandsCur <- c(AllLagVars[whichAuto],AllLagVars[-j][WindowMatrix[w,]])
# PredCur <- paste0(PredCur," + (",paste(RandsCur,collapse="+")," |",idvar,")")
#
# ff <- as.formula(paste(vars[j],"~",PredCur))
# Results[[cur]] <- lme4::lmer(ff,data=augData, control = do.call('lmerControl',control),REML=FALSE)
# formulas[[cur]] <- ff
#
# # Extract fixed effects:
# feCur <- fixef(Results[[cur]])
# FixEf[[cur]] <- feCur[names(feCur) %in% c("(Intercept)",RandsCur)]
#
# feSECur <- se.fixef(Results[[cur]])
# FixEf_SE[[cur]] <- feSECur[names(feSECur) %in% c("(Intercept)",RandsCur)]
#
# cur <- cur + 1
# }
#
# if (progress){
# setTxtProgressBar(pb, j)
# }
# }
#
# if (progress){
# close(pb)
# }
#
if (estimator=="lmmlasso"){
out <- list(
fixedEffects = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","Coef")))
# se.fixedEffects = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","se.Coef"))),
# randomEffects = ranPerID,
# randomEffectsVariance = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","Variance"))),
# pvals = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","pvals"))),
# pseudologlik = logLik,
# df = df,
# BIC = BIC,
# input = input,
# lmerResults = lapply(NodeWise_Results,"[[","lmerResult"),
# lmerFormulas = lapply(NodeWise_Results,"[[","formula")
)
class(out) <- "mlVAR"
return(out)
}
# Extract info:
if (method == "movingWindow"){
logLik <- NA
df <- NA
BIC <- NA
} else {
logLik <- sum(unlist(lapply(lapply(NodeWise_Results,"[[","Result"),logLik)))
df <- sum(unlist(lapply(lapply(lapply(NodeWise_Results,"[[","Result"),logLik),attr,'df')))
BIC <- sum(unlist(lapply(lapply(NodeWise_Results,"[[","Result"),BIC)))
}
ranlist <- lapply(NodeWise_Results,"[[","ranPerID")
ranPerID <- lapply(lapply(seq_along(ranlist[[1]]),function(i)lapply(ranlist,function(x)as.data.frame(x[[i]]))),function(xx)as.data.frame(bind_rows(xx)))
# Output:
out <- list(
fixedEffects = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","Coef"))),
se.fixedEffects = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","se.Coef"))),
randomEffects = ranPerID,
randomEffectsVariance = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","Variance"))),
pvals = as.data.frame(bind_rows(lapply(NodeWise_Results,"[[","pvals"))),
# pseudologlik = logLik,
# df = df,
# BIC = BIC,
input = input,
lmerResults = lapply(NodeWise_Results,"[[","lmerResult"),
lmerFormulas = lapply(NodeWise_Results,"[[","formula"))
class(out) <- "mlVAR0"
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.