R/create_definition_model.R

In sciarraseb/nonlinSims: Generates longitudinal data sets according to a logistic pattern.

#' Analyzes data with nonlinear latent growth curve model.
#'
#' Data are modelled with a structured latent growth curve model.
#' @md
#' @param data_wide wide version of data 
#' @param model_name name of model 
#' @export
create_definition_model <- function(data_wide, model_name) {
  
  #initial checks 
  tryCatch(expr = model_name, error = function(e) {message("Error: model_name is not a character vector")})
  
  manifest_vars <- extract_manifest_var_names(data_wide = data_wide)
  latent_vars <- c('theta', 'alpha', 'beta', 'gamma')
  obs_score_days <- paste('data.', extract_obs_score_days(data = data_wide), sep = '')

  model <- mxModel(model = model_name,
                   type = 'RAM', independent = T,
                   mxData(observed = data_wide, type = 'raw'),
                   
                   manifestVars = manifest_vars,
                   latentVars = latent_vars,
                   
                   #Residual variances; by using one label, they are assumed to all be equal (homogeneity of variance)
                   mxPath(from = manifest_vars,
                          arrows=2, free=TRUE,  labels='epsilon', values = 1, lbound = 0),
                   
                   #Set manifest means to observed means; do not enable following line of code or else convergence drops dramatically
                   #mxPath(from = 'one', to = manifest_vars, free = FALSE, arrows = 1, values = manifest_means),
                   
                   #Latent variable covariances and variances
                   mxPath(from = latent_vars,
                          connect='unique.pairs', arrows=2,
                          #aa(diff_rand), ab(cov_diff_beta), ac(cov_diff_gamma), bb(beta_rand), bc(var_beta_gamma), cc(gamma_rand)
                          free = c(TRUE,FALSE, FALSE, FALSE, 
                                   TRUE, FALSE, FALSE, 
                                   TRUE, FALSE, 
                                   TRUE), 
                          values=c(1, NA, NA, NA, 
                                   1, NA, NA, 
                                   1, NA,
                                   1),
                          labels=c('theta_rand', 'NA(cov_theta_alpha)', 'NA(cov_theta_beta)', 'NA(cov_theta_gamma)',
                                   'alpha_rand','NA(cov_alpha_beta)', 'NA(cov_alpha_gamma)', 
                                   'beta_rand', 'NA(cov_beta_gamma)', 
                                   'gamma_rand'), 
                          lbound = c(1e-3, NA, NA, NA, 
                                     1e-3, NA, NA, 
                                     1, NA,
                                     1), 
                          ubound = c(2, NA, NA, NA, 
                                     2, NA, NA, 
                                     90^2, NA, 
                                     45^2)),
                   
                   #Latent variable means (linear parameters). Note that the nonlinear parameters of beta and gamma do not have estimated means
                   mxPath(from = 'one', to = c('theta', 'alpha'), free = c(TRUE, TRUE), arrows = 1,
                          labels = c('theta_fixed', 'alpha_fixed'), lbound = 0, ubound = 7, 
                          values = c(1, 1)),
                   
                   #Functional constraints
                   mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE, 
                            labels = 'theta_fixed', name = 't', values = 1, lbound = 0,  ubound = 7), 
                   mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE, 
                            labels = 'alpha_fixed', name = 'a', values = 1, lbound = 0,  ubound = 7), 
                   mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE, 
                            labels = 'beta_fixed', name = 'b', values = 1, lbound = 1, ubound = 360),
                   mxMatrix(type = 'Full', nrow = length(manifest_vars), ncol = 1, free = TRUE, 
                            labels = 'gamma_fixed', name = 'g', values = 1, lbound = 1, ubound = 360), 
                   
                   #allows you to reduce error 
                   mxMatrix(type = 'Full', nrow = length(obs_score_days), ncol = 1, free = FALSE, 
                            labels = obs_score_days, name = 'time'),
                   
                   #Algebra specifying first partial derivatives; 
                   mxAlgebra(expression = 1 - 1/(1 + exp((b - time)/g)), name="Tl"),
                   mxAlgebra(expression = 1/(1 + exp((b - time)/g)), name = 'Al'), 
                   mxAlgebra(expression = -((a - t) * (exp((b - time)/g) * (1/g))/(1 + exp((b - time)/g))^2), name = 'Bl'),
                   mxAlgebra(expression =  (a - t) * (exp((b - time)/g) * ((b - time)/g^2))/(1 + exp((b -time)/g))^2, name = 'Gl'),
                   
                   #Factor loadings; all fixed and, importantly, constrained to change according to their partial derivatives (i.e., nonlinear functions) 
                   mxPath(from = 'theta', to = manifest_vars, arrows=1, free=FALSE,  
                          labels = sprintf(fmt = 'Tl[%d,1]', 1:length(manifest_vars))),
                   mxPath(from = 'alpha', to = manifest_vars, arrows=1, free=FALSE,  
                          labels = sprintf(fmt = 'Al[%d,1]', 1:length(manifest_vars))), 
                   mxPath(from='beta', to = manifest_vars, arrows=1,  free=FALSE,
                          labels =  sprintf(fmt = 'Bl[%d,1]', 1:length(manifest_vars))), 
                   mxPath(from='gamma', to = manifest_vars, arrows=1,  free=FALSE,
                          labels =  sprintf(fmt = 'Gl[%d,1]', 1:length(manifest_vars))), 
                   
                   mxFitFunctionML(vector = FALSE)
                   
  )
  
  return(model)
}

extract_obs_score_days <- function(data_wide) {
  
  obs_score_names <- names(data_wide)[str_detect(string = names(data_wide), pattern = 'measurement')]
  # manifest_names <- str_extract(string = names(data_wide[ , 2:ncol(data_wide)]), 
  #                               pattern = "^t\\d+_\\d+")
  return(obs_score_names)
}