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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(devtools)
#install_github("sciarraseb/nonlinSims", force=T)

#library(devtools)
library(easypackages)
packages <- c('devtools', 'nonlinSims', 'parallel', 'tidyverse', "OpenMx", "data.table", 'progress', 'tictoc')
libraries(packages)


libraries(packages)
load_all()

Experiment 1 (Does equal spacing perform best across all patterns of s-shaped change?)

To test whether equal spacing should be used when one does not know the pattern of s-shape change, Experiment 1 will manipulate the pattern of change, measurement spacing, and the number of measurements. Measurement spacing will be manipulated such that measurements are either separated by equal-length intervals, intervals that increase in length over time (time-increasing spacing), intervals that decrease in length over time (time-decreasing spacing), or intervals that separate measurements at the beginning, middle, and end of the measurement window (middle-and-extreme spacing). Number of measurements will be manipulated to either be 5, 7, 9, or 11 measurements. The patterns of s-shape change will be manipulated by modulating the midpoint parameters, thus shifting the inflection point. Change will be assumed to occur over a 360-day period. Sample size will be fixed at 225.

time_period <- 360

#fixed effects
sd_scale <- 1
common_effect_size <- 0.32
theta_fixed <- 3
alpha_fixed <- theta_fixed + common_effect_size
beta_fixed <- 180
gamma_fixed <- 20

#random effects 
sd_theta <- 0.05
sd_alpha <- 0.05
sd_beta <- 10
sd_gamma <- 4
sd_error <- 0.05
num_measurements <- 5
#Set up of population-level parameters
pop_params_4l <- generate_four_param_pop_curve(theta_fixed =  theta_fixed, alpha_fixed = alpha_fixed, 
                                               beta_fixed = beta_fixed, gamma_fixed = gamma_fixed, 
                                               sd_theta = sd_theta, sd_alpha = sd_alpha, 
                                               sd_beta = sd_beta, sd_gamma = sd_gamma, sd_error = sd_error) 

cov_matrix <- generate_four_param_cov_matrix(num_time_points = num_measurements, pop_param_list = pop_params_4l)


num_iterations <- 1e3
seed <- 27

# Experiment 1 (number measurements,  spacing, midpoint) ----------------------------------------------------------------------------------------
factor_list_exp_1 <- list('num_measurements' = seq(from = 5, to = 11, by = 2), 
                          'time_structuredness' = c('time_structured'),
                          'spacing' = c('equal', 'time_inc', 'time_dec', 'mid_ext'), 
                          'midpoint' = c(80, 180, 280), 
                          'sample_size' = 225)


tic()
exp_1_data <- run_exp_simulation(factor_list = factor_list_exp_1, num_iterations = 1, pop_params = pop_params_4l, 
                                 num_cores = detectCores()-1, seed = seed, definition = NA)
toc()


write_csv(x = exp_1_data, file = '~/Desktop/exp_1_data.csv')

# Experiment 2 (number measurements, spacing,  sample size) -------------------------------------------------------------------------------------

factor_list_exp_2 <- list('num_measurements' = seq(from = 5, to = 11, by = 2), 
                          'time_structuredness' = c('time_structured'),
                          'spacing' = c('equal', 'time_inc', 'time_dec', 'mid_ext'),
                          'midpoint' = 180, 
                          'sample_size' = c(30, 50, 100, 200, 500, 1000))

tic()
exp_2_data <- run_exp_simulation(factor_list = factor_list_exp_2, num_iterations = num_iterations, pop_params = pop_params_4l, 
                                 num_cores = detectCores(), seed = seed)
toc()

write_csv(x = exp_2_data, file = 'Desktop/exp_2_data.csv')

# Experiment 3 (number measurements, sample size, time structuredness) -------------------------------------------------------------------------

factor_list_exp_3 <- list('num_measurements' = seq(from = 5, to = 11, by = 2), 
                          'time_structuredness' = c('time_structured', 'fast_response', 'slow_response'),
                          'spacing' = c('equal'), 
                          'midpoint' = 180, 
                          'sample_size' = c(30, 50, 100, 200, 500, 1000))
tic()
exp_3_data <- run_exp_simulation(factor_list = factor_list_exp_3, num_iterations = num_iterations, pop_params = pop_params_4l, 
                                 num_cores = detectCores(), seed = seed)
toc()                

write_csv(x = exp_3_data, file = '~/Desktop/exp_3_data.csv')


factor_list_exp_def <- list('num_measurements' = seq(from = 5, to = 11, by = 2), 
                          'time_structuredness' = c('slow_response'),
                          'spacing' = c('equal'), 
                          'midpoint' = 180, 
                          'sample_size' = c(30, 50, 100, 200, 500, 1000))
tic()
exp_3_def_data <- run_exp_simulation(factor_list = factor_list_exp_def, num_iterations = 1, pop_params = pop_params_4l, 
                                 num_cores = detectCores() - 1, seed = seed, definition = T)
toc()                




write_csv(x = exp_3_data, file = '~/Desktop/exp_3_data.csv')

Code tests

#Fixed_effect starting values  
#In MPlus framework, w represents starting values (pg. 520 of manual; defaults for growth models)
##fixed-effect values for  are obtained using self-starting function and not fitting any nonlinear regression model to the data (fitting a 
#model at this stage can yield impossible values). Also note that starting value for diff_fixed is calculated by subtracting mean ending value from mean
#beginning value
##In Mplus, automatic starting values for the growth factor means and variances are generated based on individual regressions of the outcome variable on time. 
##Unfortunately, following this approach in this specific situation with nonlinear models can hijacks the starting value procedure (see above comment). 
fixed_effect_starts <- logistic_self_start(data = data)

#Random-effect starting values 
##diff_rand is obtained by obtaining variance of change scores from each person (i.e., var((score_final_time_point - score_final_time_point))
diff_rand_start <- compute_diff_rand_start(data)

rand_starts <- compute_rand_effect_starts(data = data,time_rand_effects = c('beta_rand', 'gamma_rand'))
fixed_starts <- compute_fixed_effect_starts(data = data, rand_effect_starts = rand_starts)


#w  + srb
##r = unif(-0,5, 0.5) 
##s = 5 (determines strength of random perturbation)
##b = base scale variable = 2max{sqrt(var), 1}, where var = starting value for variance of parameter
###b for diff_fixed = 2 (almost all of the time)
###b for beta_fixed = sqrt(var) almost all of the time 
###b for gamma_fixed = sqrt(var) almost all of the time
###b for random effects and residual variance = 0 (i.e., no perturbation in generation process of random variables or residual variance)

#Residual variance starting value (Hipp & Bauer, 2006)
##20-80% of observed variance (averaged across each time point)
res_var <- compute_residual_start_value(data_wide = data_wide)

Appendix

Apppendix A: Computation of spacing for time-increasing and time-decreasing conditions

Several attempts were made to compute interval lengths such that they were of integer lengths. Unfortunately, it was not possible to find a set of intervals that followed the equation $i_{length} = kc + b$ (where $c$ represents the constant length and $b$ represents the base length) such that all interval lengths were of integer values. The code below shows that, for all base lengths of 1--30 days (assuming a time period of 360 days), no set of intervals is only integes.

compute_time_increasing_intervals <- function(time_period, num_measurements, base_time_length) {

  #compute length of constant by first calculating how many days remain after subtracting base_time_length for each interval.
  ##num_measurements-1 = number of intervals
  remaining_num_days <- time_period - (num_measurements-1)*base_time_length
  ##The number of constants = num_measurements - 1
  constant_length <- remaining_num_days/sum(seq(0,(num_measurements-2)))

  interval_lengths <- seq(0,(num_measurements-2))*constant_length + base_time_length
  return(interval_lengths)
}

testx <- function() {
  for (base_time_length in seq(1:30)) {
    print(base_time_length)
    print("=================")
    for (num_measurements in c(5, 7, 9, 11)) {
      xo <- get_stuff(n = n, x = x)
      xo <- c(n, xo, sum(xo))
      print(sprintf("%1.3f",xo))
    }
    print("=================")
  }
}

Appendix B: Other functions

sciarraseb/nonlinSims documentation built on Jan. 30, 2023, 8:17 p.m.