R/main.R

In jhorzek/addModeratorToCtsem: Adds moderators to ctsemOMX objects

#' addModerator
#'
#' Adds moderators to a ctsem object. The general structure is y = (a1 + a2*moderatorValues[,1] + a3*moderatorValue[,2] + ...)*predictor
#'
#' @param model model of type ctsemFit or MxModel
#' @param moderatedVariable label of the variable which is to be moderated (MUST be identical to the label of the parameter in the ctsemFit or MxModel object)
#' @param variableInMatrix name of the mxMatrix where the moderatedVariable is located
#' @param moderatorValues matrix with nrow = sample size and ncol = number of moderators. Contains for each person the values of each moderator. The colnames of this matrix will be used to label the moderators.
#' @returns an MxModel with moderator on the specific variable. The main effect is labeled moderatedVariable_MAIN, the interaction effects are in moderatedVariable_INTERACTION
#' @examples
#' library(ctsemOMX)
#' library(addModeratorToCtsem)
#'
#' #### Example: ####
#' # Model: dx(t) = (a + a1*mod1 + a2*mod2)*x(t)dt + g*W(t)dt
#' n <- 100 # sample size
#' moderators <- matrix(rnorm(2*n,0,.05), ncol = 2) # values of the moderator variables
#'
#' a <- -.2 # main effect
#' aInd <- -.2 - 1*moderators[,1] - .5*moderators[,2] # individual drift values (a+a1*mod1+a2*mod2)
#' hist(aInd)
#'
#' # simulate data for a first order stochastic differential equation
#' data <- c()
#' for(i in 1:n){
#'   generatingModel<-ctModel(Tpoints=100,n.latent=1,n.TDpred=0,n.TIpred=0,n.manifest=1,
#'                            MANIFESTVAR=diag(0,1),
#'                            LAMBDA=diag(1,1),
#'                            DRIFT=matrix(c(aInd[i]),nrow=1),
#'                            DIFFUSION=matrix(c(1),1),
#'                            T0MEANS=matrix(0,ncol=1,nrow=1),
#'                            T0VAR=diag(1,1))
#'   dat <- ctGenerate(generatingModel,n.subjects=1,burnin=10)
#'   dat[,"id"] <- i
#'   data <-rbind(data, dat)
#' }
#'
#' # set up model with ctsem
#' model<-ctModel(Tpoints=100,n.latent=1,n.TDpred=0,n.TIpred=0,n.manifest=1,
#'                MANIFESTVAR=diag(0,1),
#'                LAMBDA=diag(1,1),
#'                DRIFT=matrix(c("a"),nrow=1),
#'                DIFFUSION=matrix("g",1),
#'                T0MEANS=matrix(0,ncol=1,nrow=1),
#'                T0VAR=diag(1,1))
#'
#' # fit model. IMPORTANT: objective must be set to "Kalman"
#' fitModel <- ctFit(dat = data,
#'                   ctmodelobj = model,
#'                   fit = T,
#'                   objective = "Kalman")
#'
#' # create moderated model
#' moderatedModel <- addModerator(model = fitModel,
#'                                moderatedVariable = "a",
#'                                variableInMatrix = "DRIFT",
#'                                moderatorValues = moderators)
#'
#' # fit moderated model
#' fit.moderatedModel <- mxRun(moderatedModel)
#'
#' # the main effect is always labeled moderatedVariable_MAIN
#' fit.moderatedModel$a_MAIN
#' # the interaction effects are always labeled moderatedVariable_INTERACTION
#' fit.moderatedModel$a_INTERACTION
#' @export

addModerator <- function(model, moderatedVariable, variableInMatrix, moderatorValues){
  if(class(model) == "ctsemFit"){
    model <- model$mxobj
  }else if(class(model) != "MxModel"){
    stop("model has to be of type ctsemFit or MxModel")
  }
  if(anyNA(moderatorValues)) stop("NA in moderatorValues detected. This function assumes that all moderator values are observed and without error.")

  if(!variableInMatrix %in% names(model)) stop(paste0("Could not find ", variableInMatrix, " in names(model)."))
  if(!moderatedVariable %in% model[[variableInMatrix]]$labels) stop(paste0("Could not find ", moderatedVariable, " in model[[variableInMatrix]]$labels."))

  if(!"id" %in% colnames(model$data$observed)){
    stop("The data set must have a variable called id which identifies the individuals. If you are using ctsemOMX, make sure to pass objective = 'Kalman' to ctFit.")
  }

  if(is.vector(moderatorValues)) moderatorValues <- matrix(moderatorValues, ncol = 1)

  if(nrow(moderatorValues) != length(unique(model$data$observed[,"id"]))){stop("The moderatorValues must be of length N, where N is the sample size.")}

  if(is.null(colnames(moderatorValues))){
    moderatorNames <- paste0("moderator", 1:ncol(moderatorValues))
    warning(paste0("The moderatorValues matrix does not have column names. The moderators will be labeled: ", moderatorNames, collapse = ","))
  }else{
    moderatorNames <- colnames(moderatorValues)
  }

  nModerators <- ncol(moderatorValues)

  # Set up model
  newModel <- mxModel(model)
  selectElement <- model[[variableInMatrix]]$labels == moderatedVariable

  newModel[[variableInMatrix]]$labels[selectElement] <- paste0("moderationAlgebra_", moderatedVariable, "[1,1]")
  newModel[[variableInMatrix]]$free[selectElement] <- F

  mainEffect <- mxMatrix(values = newModel[[variableInMatrix]]$values[selectElement],
                         labels = paste0(moderatedVariable,"_main"),
                         free = T, nrow = 1, ncol = 1,
                         name = paste0(moderatedVariable,"_MAIN"))

  interactionEffect <- mxMatrix(values = rep(0, nModerators),
                                labels = matrix(paste0(moderatedVariable,"_X_", moderatorNames), ncol = 1),
                                free = T,
                                nrow = nModerators, ncol = 1,
                                name = paste0(moderatedVariable,"_INTERACTION"))
  mxModerators <- mxMatrix(values = moderatorValues, nrow = nrow(moderatorValues),
                           ncol = nModerators,
                           free = FALSE,
                           name = paste0("moderatorsFor", moderatedVariable))

  moderationAlgebra <- mxAlgebraFromString(algString = paste0(paste0(moderatedVariable,"_main"),
                                                              "+",
                                                              paste0("moderatorsFor", moderatedVariable), "[data.id,]",
                                                              " %*% ",
                                                              paste0(moderatedVariable,"_INTERACTION")
  ),
  name = paste0("moderationAlgebra_", moderatedVariable)
  )

  newModel <- mxModel(newModel,
                      mainEffect,
                      interactionEffect,
                      mxModerators,
                      moderationAlgebra)
  return(newModel)
}