Simulation_Code.R

#####################################
# R Code for Simulation
#####################################

# Batra, R.*, Johal, S.K.*, Chen, M., Ferrer, E. (in prep.). 
# Consequences of sampling frequency on the estimated dynamics of AR processes using continuous time models.

# Load Libraries ----------------------------------------------------------
remove(list = ls())
library(dplyr)
library(data.table)
library(OpenMx) # OpenMx was used to fit the CT and DT models. See Note 1 in the end for comments on other packages for estimation.
library(tidyverse)
library(ggh4x)

# Condition Lists ---------------------------------------------------------

# This is a data frame for all the conditions in the simulation:
# True_Fluctuation: This is the data generating frequency that is of interest for capturing the AR parameter. We use Daily and Weekly for our simulation.
# AR_Parameter: This is the strength of the autoregressive parameter at a particular data generating frequency. We use four different values of .05, .2, .5 and .8.
# CT_Autoeffect: We generate the data at hourly time interval. To generate hourly data, we need the autoeffect for the hourly interval that corresponds to a particular data generating frequency and AR parameter. See Data Generation section of the manuscript for more details.
# Sampling: We use two types of equality of sampling intervals. SameTimestamp is for equal intervals, where the data point is observed at the same time every day/week/month. RandomTimestamp is for random intervals, where the data point is observed at any point within a day/week/month.
# Sample_Size: This is the condition for the number of observations sampled for a single individual. Can be either 100, 300 or 1000. 1000 serves as the asymptotic level for the number of observations. 

ct_conditions = data.frame("True_Fluctuation" = rep(c("Daily", "Weekly"), each = 8, times = 3),
                           "AR_Parameter" = rep(c(0.05, 0.2, 0.5, 0.8), 12),
                           "CT_Autoeffect" = rep(c(rep(c(-0.1248222, -0.06706, -0.02888113, -0.00923), 2), 
                                                   rep(c(-0.01783174, -0.009579988, -0.004125876, -0.001328235), 2)), 3),
                           "Sampling" = rep(rep(c("SameTimestamp", "RandomTimestamp"), each = 4), 6),
                           "Sample_Size" = rep(c(100, 300, 1000), each = 16))

# We use equation 3 to generate hourly data, where the error variance uses the value of the autoeffect and the variance of process y (given under the section of Base Time Scale). var(y) is set outisde of the simulation.

var_asym <- 0.25 
ct_conditions$e_var = var_asym * (1-exp(2*ct_conditions$CT_Autoeffect*1))

######################################################################
# Generate Data ------------------------------------------------------
######################################################################

# true_process: (same as before) This is the data generating frequency that is of interest for capturing the AR parameter. We use Daily and Weekly for our simulation.
# ar_parameter: (same as before) This is the strength of the autoregressive parameter at a particular data generating frequency. We use four different values of .05, .2, .5 and .8.
# total_months: This argument helps determine how many hours of data to generate for the number of months because that is the highest time metric for sampling. So if we need X number of months observations, we need a 672*X amount of hourly data
# setseed_error: We set the seed for errors that are derived from a normal distribution for each iteration so the simulation can be replicated by others.

generate_data = function(true_process, ar_parameter, total_months, setseed_error){
  
  #### Generate Data ####
  
  # get the ct effect
  
  ct_auto = ct_conditions %>% filter(True_Fluctuation == true_process & AR_Parameter == ar_parameter) %>%
    distinct(CT_Autoeffect) %>%
    as.numeric()
  
  # get discrete time hourly effect using formula
  
  hour_dt = exp(ct_auto * 1)
  
  # generate vector of errors
  # (total_months*672 - 1) values instead of (total_months*672) bc the first value (y0)
  # has no error
  
  total_hours = total_months*672
  
  e_var = ct_conditions$e_var[ct_conditions$True_Fluctuation == true_process & 
                                ct_conditions$AR_Parameter == ar_parameter][1]
  
  set.seed(setseed_error) # 500 seeds for 500 replications
  
  error_vec = rnorm(total_hours - 1, mean = 0, sd = sqrt(e_var))
  
  # generate full dataset
  
  Y = c(0, rep(NA, total_hours - 1))
  
  for(i in 2:(total_hours - 1)) {
    Y[i] = hour_dt*Y[i-1] + error_vec[i-1]
  }
  
  return(Y) 
  
}

######################################################################
# Sampling Generated Data Functions ----------------------------------
######################################################################

# dataset: This is the time series generated by generate_data
# samplingfreq: We sample the generated data at three different frequencies, that is daily, weekly, or monthly
# startindex: This argument tells us where to start sampling from in the process (we have it set to 1, so starting the sampling from the first hour)
# observations: This is the condition for the number of observations sampled for a single individual. Can be either 100, 300 or 1000. For example, if we want to sample 100 daily observations, we will need 100*24 hours of data from which to sample the data from.
# setseed_sampling: When we are sampling randomly, we also set the seed for sampling as it helps in future replications picking the same group of observations. 

# Sampling Randomly

sample_randomly = function(dataset, samplingfreq, startindex, observations, setseed_sampling){
  
  if(samplingfreq == "Daily"){
    
    total_hours = observations*24
    
    # check that we have enough data in our dataset to get the number of observations
    
    if((startindex + total_hours - 1) > length(dataset)){
      print("Not enough data")
      stop()
    }else{
      data_subset = dataset[startindex:(startindex + total_hours - 1)] 
    }
    
    # now create groupings of 24 hours
    # timestamp is so we can keep the exact time elapsed between
    # measurement occasions; however, divided by 24 because we want
    # to keep the timestamp in the metric of the sampling frame
    # not the hourly metric
    
    data_chunked = data.frame("Id" = rep(1, total_hours), "Measurement_Occasion" = rep(1:observations, each = 24),
                              "Value" = data_subset,
                              "Timestamp" = 1:total_hours, "Base_Time" = (1:total_hours)/24)
    
    # now from each "day", choose one observation at random, and
    # return its timsestamp
    set.seed(setseed_sampling)
    observations = data_chunked %>%
      group_by(Measurement_Occasion) %>%
      slice_sample(n = 1)
    
  }else if(samplingfreq == "Weekly"){
    
    total_hours = observations*168
    
    # check that we have enough data in our dataset to get the number of observations
    
    if((startindex + total_hours - 1) > length(dataset)){
      print("Not enough data")
      stop()
    }else{
      data_subset = dataset[startindex:(startindex + total_hours - 1)] 
    }
    
    # now create groupings of 168 hours
    # timestamp is so we can keep the exact time elapsed between
    # measurement occasions
    
    data_chunked = data.frame("Id" = rep(1, total_hours), "Measurement_Occasion" = rep(1:observations, each = 168),
                              "Value" = data_subset,
                              "Timestamp" = 1:total_hours, "Base_Time" = (1:total_hours)/168)
    
    # now from each "week", choose one observation at random, and
    # return its timestamp
    set.seed(setseed_sampling)
    observations = data_chunked %>%
      group_by(Measurement_Occasion) %>%
      slice_sample(n = 1)
    
  }else if(samplingfreq == "Monthly"){
    
    total_hours = observations*672
    
    # check that we have enough data in our dataset to get the number of observations
    
    if((startindex + total_hours - 1) > length(dataset)){
      print("Not enough data")
      stop()
    }else{
      data_subset = dataset[startindex:(startindex + total_hours - 1)] 
    }
    
    data_chunked = data.frame("Id" = rep(1, total_hours), "Measurement_Occasion" = rep(1:observations, each = 672),
                              "Value" = dataset,
                              "Timestamp" = 1:total_hours, "Base_Time" = (1:total_hours)/672)
    
    # now from each "month", choose one observation at random, and
    # return its index
    set.seed(setseed_sampling)
    observations = data_chunked %>%
      group_by(Measurement_Occasion) %>%
      slice_sample(n = 1)
  }
  
  observations <- as.data.frame(observations)
  return(observations)
  
}

# Sampling at the Same Timestamp

sampling_same_timestamp <- function(dataset, samplingfreq, timestamp, observations, start_index){
  
  # Each hourly data starts at 12 am so the hourly data is for starting of the hour
  # so first observation is 12 am, 2nd observation is 1 am,..., 10th obs is 9 am
  # timestamp is the hour when we start observing the process, which is the same for each day, week, month
  # start_index is the index where the window of our observation period starts
  # Base_Time is the time 
  
  if(samplingfreq == "Daily"){
    
    total_hours = observations*24
    
    # check that we have enough data in our dataset to get the number of observations
    
    if((start_index + total_hours - 1) > length(dataset)){
      print("Not enough data")
      stop()
    }else{
      new_dataset <- dataset[c(start_index: (start_index + total_hours - 1))] 
    }
    
    sample_index <- seq(from = timestamp, to = total_hours, by = 24)
    sample_data <- data.frame("Id" = rep(1, observations), "Measurement_Occasion" = c(1:observations), "Value" = new_dataset[sample_index],
                              "Timestamp" = sample_index, "Base_Time" = c(1:observations))
  }else if(samplingfreq == "Weekly"){
    
    # We have 7 days in a week here, so 7*24 = 168 will be the increment for each value at
    # the exact same time/day each month
    
    total_hours = observations*168
    
    # check that we have enough data in our dataset to get the number of observations
    
    if((start_index + total_hours - 1) > length(dataset)){
      print("Not enough data")
      stop()
    }else{
      new_dataset <- dataset[c(start_index: (start_index + total_hours - 1))] 
    }
    
    sample_index <- seq(from = timestamp, to = total_hours, by = 168)
    sample_data <- data.frame("Id" = rep(1, observations), "Measurement_Occasion" = c(1:observations), "Value" = new_dataset[sample_index],
                              "Timestamp" = sample_index, "Base_Time" = c(1:observations))
    
  }else if(samplingfreq == "Monthly"){
    
    # We have 28 days in a month here, so 7*24*4 = 672 will be the increment for each value at
    # the exact same time/day each month
    
    total_hours = observations*672
    
    # check that we have enough data in our dataset to get the number of observations
    
    if((start_index + total_hours - 1) > length(dataset)){
      print("Not enough data")
      stop()
    }else{
      new_dataset <- dataset[c(start_index: (start_index + total_hours - 1))] 
    }
    
    sample_index <- seq(from = timestamp, to = total_hours, by = 672)
    sample_data <- data.frame("Id" = rep(1, observations), "Measurement_Occasion" = c(1:observations), "Value" = dataset[sample_index],
                              "Timestamp" = sample_index, "Base_Time" = c(1:observations))
  }
  
  return(sample_data)
}

######################################################################
# Estimation of the Model (OpenMx) -----------------------------------
######################################################################

# CT-AR Model

estimate_model_openmx <- function(dataset, cond_no, iter_no, samplingfreq){
  
  # Empty List for Model Specification
  opmxL = list()

  opmxL$bmat <- mxMatrix(name = "B", "Zero", 1, 1) # No Covariates in the model

  opmxL$cmat <- mxMatrix(name = "C", "Full", nrow = 1, ncol = 1, free = FALSE,
                       values = c(1), 
                       dimnames = list(c("Value"), c("Y"))) # No latent variables/Factor Loading is 1

  opmxL$dmat <- mxMatrix("Zero", 1, 1, name = "D") 

  opmxL$qmat <- mxMatrix("Full", 1, 1, free = TRUE, 
                       values = c(1),
                       name = "Q", 
                       labels = c("error_var")) # Dynamic residual variance of the process

  opmxL$rmat <- mxMatrix("Zero", 1, 1,
                       name = "R") 

  opmxL$xmat <- mxMatrix(name = "x0", "Full", 1, 1, free = FALSE,
                       values = c(0),
                       labels = c("Y_Mean"))

  opmxL$pmat <- mxMatrix(name = "P0", "Full", 1, 1, free = FALSE,
                       values = c(0.25),
                       labels = c("var_Y"))

  opmxL$umat <- mxMatrix("Zero", 1, 1, name = "u") 

  opmxL$tmat <- mxMatrix('Full', 1, 1, name='time', labels = "data.Base_Time")

  check = 0
  
  for (k in 1:10) {
    set.seed(k*123)
    betstart <- -runif(1, 0 ,1)
  
    opmxL$amat_ct <- mxMatrix(name = "A", "Full", 1, 1, free = TRUE,
                          values = c(betstart),
                          labels = c("beta_y")) # dynamic model
    
    # Model Estimation
    ar_ct_model = mxModel(model = "AR(1) Model",
                          opmxL,
                          mxFitFunctionML(),
                          mxData(dataset, type = "raw"),
                          mxExpectationStateSpaceContinuousTime(A = "A", B = "B", C = "C", 
                                                            D = "D", Q = "Q", R = "R",
                                                            x0 = "x0", P0 = "P0",
                                                            u = "u", t = "time"))

    fit_ar_ct = try(mxRun(ar_ct_model))
    # fit_ar_ct = mxTryHard(ar_ct_model)
    
    if(class(fit_ar_ct)[1] == "MxModel"){
      if(all(!is.na(fit_ar_ct@output$standardErrors)) & (fit_ar_ct@output$status[[1]] == 0)){
      check = 1
      fit_mod_summary <- summary(fit_ar_ct)
      
      # Calculation of Wald CIs
      n = dim(dataset)[1]
      lb_error_var = fit_mod_summary$parameters$Estimate[1] - qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[1]
      up_error_var = fit_mod_summary$parameters$Estimate[1] + qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[1]
      
      lb_beta = fit_mod_summary$parameters$Estimate[2] - qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[2]
      ub_beta = fit_mod_summary$parameters$Estimate[2] + qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[2]
      
      fit_mod_summary$parameters$lbound = c(lb_error_var, lb_beta)
      fit_mod_summary$parameters$ubound = c(up_error_var, ub_beta)
      
      # Extracting the Estimates from the Model
      fit_mod_data <- data.frame("parameter" = c("error_var_ct", "beta_ct", "LogLik_ct", "AIC_ct", "BIC_ct"), 
                                 rbind(fit_mod_summary$parameters[, c(5, 6, 7, 8)], 
                                       c(fit_mod_summary$fit, rep(NA, 3)), 
                                       c(fit_mod_summary$informationCriteria[1, 2], rep(NA, 3)), 
                                       c(fit_mod_summary$informationCriteria[2, 2], rep(NA, 3))), 
                                 "Method" = "CT", "Check_of_Convergence" = check,
                                 row.names = NULL)
      
      break
    }else{
      fit_mod_summary <- summary(fit_ar_ct)
      
      # No calculation of CIs needed here because these fits are not used for analysis
      
      fit_mod_data <- data.frame("parameter" = c("error_var_ct", "beta_ct", "LogLik_ct", "AIC_ct", "BIC_ct"), 
                                 rbind(fit_mod_summary$parameters[, c(5, 6, 7, 8)], c(fit_mod_summary$fit, rep(NA, 3)), 
                                       c(fit_mod_summary$informationCriteria[1, 2], rep(NA, 3)), 
                                       c(fit_mod_summary$informationCriteria[2, 2], rep(NA, 3))), 
                                 "Method" = "CT", "Check_of_Convergence" = check,
                                 row.names = NULL)
    }
  }else{
    fit_mod_data <- data.frame("parameter" = c("error_var_ct", "beta_ct", "LogLik_ct", "AIC_ct", "BIC_ct"), 
                               matrix(rep(NA, 20), nrow = 5), 
                               "Method" = "CT", "Check_of_Convergence" = check,
                               row.names = NULL)
    colnames(fit_mod_data) <- c("parameter", "Estimate", "Std.Error", "lbound", "ubound", "Method", 
                                "Check_of_Convergence")
  }}
  
  assign(paste("cook_Cond", cond_no, "_iter_", iter_no, "_", samplingfreq, "CT", sep = ""), fit_ar_ct)

  save(list = paste("cook_Cond", cond_no, "_iter_", iter_no, "_", samplingfreq, "CT", sep = ""),
       file = paste("/Users/rbat/Box/Fluctuations in CT/shared CT Group Project/RData_Files_(OpenMx)/Cook_Cond",
                    cond_no, "_iter_", iter_no, "_", samplingfreq, "CT.RData", sep = ""))
  
  return(fit_mod_data)
}

# DT-AR Model

estimate_model_OpenMx_dt = function(dataset, cond_no, iter_no, samplingfreq){
  
  # Empty List for Model Specification
  opmxL = list()
  
  opmxL$bmat <- mxMatrix(name = "B", "Zero", 1, 1) # No Covariates in the model
  
  opmxL$cmat <- mxMatrix(name = "C", "Full", nrow = 1, ncol = 1, free = FALSE,
                         values = c(1), 
                         dimnames = list(c("Value"), c("Y"))) # No latent variables/Factor Loading is 1
  
  opmxL$dmat <- mxMatrix("Zero", 1, 1, name = "D") 
  
  opmxL$qmat <- mxMatrix("Full", 1, 1, free = TRUE, 
                         values = c(1),
                         name = "Q", 
                         labels = c("error_var")) # Dynamic residual variance of the process
  
  opmxL$rmat <- mxMatrix("Zero", 1, 1,
                         name = "R") 
  
  opmxL$xmat <- mxMatrix(name = "x0", "Full", 1, 1, free = FALSE,
                         values = c(0),
                         labels = c("Y_Mean"))
  
  opmxL$pmat <- mxMatrix(name = "P0", "Full", 1, 1, free = FALSE,
                         values = c(0.25),
                         labels = c("var_Y"))
  
  opmxL$umat <- mxMatrix("Zero", 1, 1, name = "u") 

  
  check = 0
  
  for (k in 1:10) {
    set.seed(k*123)
    betstart <- runif(1, 0 ,1)
    
    opmxL$amat_dt <- mxMatrix(name = "A", "Full", 1, 1, free = TRUE,
                              values = c(betstart),
                              labels = c("beta_y")) # dynamic model
    
    # Model Estimation
    ar_dt_model = mxModel(model = "AR(1) Model",
                          opmxL,
                          mxFitFunctionML(),
                          mxData(dataset, type = "raw"),
                          mxExpectationStateSpace(A = "A", B = "B", C = "C", 
                                                                D = "D", Q = "Q", R = "R",
                                                                x0 = "x0", P0 = "P0",
                                                                u = "u"))
    
    fit_ar_dt = try(mxRun(ar_dt_model))
    # fit_ar_ct = mxTryHard(ar_ct_model)
    
    if(class(fit_ar_dt)[1] == "MxModel"){
      if(all(!is.na(fit_ar_dt@output$standardErrors)) & (fit_ar_dt@output$status[[1]] == 0)){
        check = 1
        fit_mod_summary <- summary(fit_ar_dt)
        
        # Calculation of Wald CIs
        n = dim(dataset)[1]
        lb_error_var = fit_mod_summary$parameters$Estimate[1] - qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[1]
        up_error_var = fit_mod_summary$parameters$Estimate[1] + qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[1]
        
        lb_beta = fit_mod_summary$parameters$Estimate[2] - qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[2]
        ub_beta = fit_mod_summary$parameters$Estimate[2] + qt(0.025, n, lower.tail = FALSE)*fit_mod_summary$parameters$Std.Error[2]
        
        fit_mod_summary$parameters$lbound = c(lb_error_var, lb_beta)
        fit_mod_summary$parameters$ubound = c(up_error_var, ub_beta)
        
        # Extracting the Estimates from the Model
        fit_mod_data <- data.frame("parameter" = c("error_var_dt", "beta_dt", "LogLik_dt", "AIC_dt", "BIC_dt"), 
                                   rbind(fit_mod_summary$parameters[, c(5, 6, 7, 8)], c(fit_mod_summary$fit, rep(NA, 3)), 
                                         c(fit_mod_summary$informationCriteria[1, 2], rep(NA, 3)), 
                                         c(fit_mod_summary$informationCriteria[2, 2], rep(NA, 3))), 
                                   "Method" = "DT", "Check_of_Convergence" = check,
                                   row.names = NULL)
        
        break
      }else{
        fit_mod_summary <- summary(fit_ar_dt)
        
        fit_mod_data <- data.frame("parameter" = c("error_var_dt", "beta_dt", "LogLik_dt", "AIC_dt", "BIC_dt"), 
                                   rbind(fit_mod_summary$parameters[, c(5, 6, 7, 8)], c(fit_mod_summary$fit, rep(NA, 3)), 
                                         c(fit_mod_summary$informationCriteria[1, 2], rep(NA, 3)), 
                                         c(fit_mod_summary$informationCriteria[2, 2], rep(NA, 3))), 
                                   "Method" = "DT", "Check_of_Convergence" = check,
                                   row.names = NULL)
      }
    }else{
      fit_mod_data <- data.frame("parameter" = c("error_var_dt", "beta_dt", "LogLik_dt", "AIC_dt", "BIC_dt"), 
                                 matrix(rep(NA, 20), nrow = 5), 
                                 "Method" = "DT", "Check_of_Convergence" = check,
                                 row.names = NULL)
      colnames(fit_mod_data) <- c("parameter", "Estimate", "Std.Error", "lbound", "ubound", "Method", 
                                  "Check_of_Convergence")
    }}
  
  assign(paste("cook_Cond", cond_no, "_iter_", iter_no, "_", samplingfreq, "DT", sep = ""), ar_dt_model)

  save(list = paste("cook_Cond", cond_no, "_iter_", iter_no, "_", samplingfreq, "DT", sep = ""),
       file = paste("/Users/rbat/Box/Fluctuations in CT/shared CT Group Project/RData_Files_(OpenMx)/Cook_Cond",
                    cond_no, "_iter_", iter_no, "_", samplingfreq, "DT.RData", sep = ""))
  
  return(fit_mod_data)
  
}


# Higher-level Simulation Function -----------------------------------------------------

# true_process is the time metric that the true process fluctuates at
# ar_parameter is strength of DT AR Parameter on the true metric
# sampling is whether sampling occurs randomly or at the same time
# observations tells us how many observations we are sampling (e.g., 100 days/weeks/months)
# however, we need to make sure we generate enough data to observe that many months, so we also use
# observations to guide how many months worth of data we generate, 
# set_seed_error and set_seed_y are used to set the seed for data generation


fluctuations_ct_simulation = function(true_process, ar_parameter, sampling, observations,
                                      set_seed_error, setseed_sampling, cond_no, iter_no){
  
  #### Generate Data ####
  
  # total_months determined by how many observations we want to sample, because we need to make sure we have enough
  Y = generate_data(true_process = true_process, ar_parameter = ar_parameter, total_months = observations,
                    setseed_error = set_seed_error)
  
  #### Sampling ####
  
  if(sampling == "RandomTimestamp"){

    daily_sample = sample_randomly(Y, "Daily", 1, observations, setseed_sampling)

    weekly_sample = sample_randomly(Y, "Weekly", 1, observations, setseed_sampling )
    
    monthly_sample = sample_randomly(Y, "Monthly", 1, observations, setseed_sampling)
    
  }else if(sampling == "SameTimestamp"){
    
    daily_sample = sampling_same_timestamp(Y, "Daily", 10, observations, 1)
    
    weekly_sample = sampling_same_timestamp(Y, "Weekly", 10, observations, 1)
    
    monthly_sample = sampling_same_timestamp(Y, "Monthly", 10, observations, 1)
    
  }
  
  #### Fitting Models ####
  
  fit_model_daily = cbind("Sample_Freq" = "Daily",
                          rbind(estimate_model_openmx(daily_sample, cond_no = cond_no, iter_no = iter_no, "Daily"),
                                estimate_model_OpenMx_dt(daily_sample, cond_no = cond_no, iter_no = iter_no, "Daily")))
  
  fit_modely_weekly = cbind("Sample_Freq" = "Weekly", 
                            rbind(estimate_model_openmx(weekly_sample, cond_no = cond_no, iter_no = iter_no, "Weekly"),
                                  estimate_model_OpenMx_dt(weekly_sample, cond_no = cond_no, iter_no = iter_no, "Weekly")))
  
  fit_model_monthly = cbind("Sample_Freq" = "Monthly", 
                            rbind(estimate_model_openmx(monthly_sample, cond_no = cond_no, iter_no = iter_no, "Monthly"),
                                  estimate_model_OpenMx_dt(monthly_sample, cond_no = cond_no, iter_no = iter_no, "Monthly")))
  
  #### Return Statement ####
  
  return(rbind(fit_model_daily, fit_modely_weekly, fit_model_monthly))
  
}

# data frame for the seeds
df_seeds <- data.frame("error" = c(150:649), "sample" = c(1550:2049))

# Higher level simulation

for (j in c(1:48)) {
  true_process = ct_conditions[j, "True_Fluctuation"]
  ar_parameter = ct_conditions[j, "AR_Parameter"]
  sampling = ct_conditions[j, "Sampling"]
  observations = ct_conditions[j, "Sample_Size"]
  data_1 <- data.frame(NULL)
  for(i in 1:500){
    set_seed_error = df_seeds[i, "error"]
    setseed_sampling = df_seeds[i, "sample"]
    iter_data <- fluctuations_ct_simulation(true_process, ar_parameter, sampling, observations, set_seed_error, 
                                            setseed_sampling, j, i)
    data_1 = rbind(data_1, cbind(iter_data, "Iteration" = i)) 
  }
  data_1 = cbind("True_Fluctuation" = true_process, "AR_Parameter" = ar_parameter, "Sampling" = sampling,
                 "Sample_Size" = observations, data_1)
  write.csv(data_1, 
            file = paste("/Users/rbat/Box/Fluctuations in CT/shared CT Group Project/Results_OpenMx/data_cond", 
                         j, ".csv", sep = ""), row.names = FALSE)
}

# Note 1: We also fit the models using Dynr package in R (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8297742/) but we discovered that for this specific group of conditions, the estimation was highly biased when the AR parameter is less than .3. This is because dynr uses a numerical solver and for this specific conditions of a CT-AR model used here, it faces issues with estimation.