Computing speed slows down as iteration increases

[xanalee]

Hi~

I am currently doing GLMM estimation on 1000 simulated datasets by brms and evaluating the computational cost. However, I find the compute time becomes longer and longer as interation increases. The two figures below show the computational slowdown observed when I running the experiments locally and on ALICE HPC platform, respectively.

Here are the versions of R and brms I used:

• Local environment: R version 4.4.3, brms version 2.22.0
• ALICE HPC platform: R version 4.4.0, brms version 2.22.0

And I am pretty sure the simulated datasets are homogeneous since they are generated by Monte Carlo simulation. The super smooth curves make me guess there might be accumulation of internal variables or uncleaned caches. Could you please help me to figure out the issue?

I give my script here for reference.

library(brms)

path = ''  # edit as appropriate

n_cores = 4
n_threads = 2

sim_brms = function(data_ls, num_sub, num_time, num_sim, resp_dist){
  out_path = paste0(path, 'brms_', 'n', num_sub, '_t', num_time, '_sim', num_sim, '_', resp_dist, '/')
  if (!dir.exists(out_path)) {
    dir.create(out_path, recursive = TRUE, showWarnings = FALSE)
  }
  for (i in 1:num_sim){
    # Monitor
    cat(paste('brms_i: ', i, '\n'))
    
    data = data_ls[[i]]
    
    # Define helper functions
    check_convergence = function(model, n_term) {
      if (inherits(model, "try-error")) return("error")
      Rhat = rhat(model)
      if (any(abs(Rhat[1:n_term] - 1) >= 0.1)) return("official_fail")
      "good"
    }
    
    # Specify response distribution
    if (resp_dist == 'binomial'){
      fami = bernoulli(link = 'logit')
    }
    if (resp_dist == 'Poisson'){
      fami = poisson(link = 'log')
    }
    
    # Fit full model on full dataset
    compute_time = system.time({
      brms_model_full = try(brm(y_full ~ time + group + time:group + (1|id),
                                family = fami,
                                data = data,
                                threads = threading(n_threads),
                                cores = n_cores,
                                seed = 123
      ))})[['elapsed']]
    conv_status = check_convergence(brms_model_full, 5)
    out = c('compute_time' = compute_time, 'converge_rate' = conv_status)
    save(out, file = paste0(out_path, i, '.RData'))
    gc()
  }
}

n_sim = 1000
n_sub = 20
n_time = 2
re_dist = 'binomial'
cat(paste0('n', n_sub, '_t', n_time, '_', re_dist, '\n'))
load(paste0(path, 'data/n', n_sub, '_t', n_time, '_sim', n_sim, '_', re_dist,  '_rdi.RData'))
sim_brms(data_ls = simulate_data_ls, num_sub = n_sub, num_time = n_time, num_sim = n_sim, resp_dist = re_dist)

Thank you in advance and look forward for a reply.

[mirkosignorelli]

In order to be able to run the code in the example above, I am adding here the script used to generate the datasets in the code given above:

library(MASS)
path = '' # edit as appropriate

gen_data = function(num_sub, num_time, num_sim, resp_dist, with_rd_slope) {
  
  n_obs = num_sub * num_time
  
  # Subject ID
  subject_id = factor(rep(1:num_sub, each = num_time))
  
  # Time variable
  time_variable = rep((0:(num_time-1))/(num_time-1), times = num_sub)
  
  simulate_data_ls = list()
  for (i in 1:num_sim) {
    # Group variable
    group = rbinom(num_sub, size = 1, prob = 0.4)
    group_variable = group[subject_id]
    
    # Fixed effects
    beta0 = 0.1
    beta1 = 0.3
    beta2 = -0.2
    beta3 = 0.1
    # z of full dataset
    z_full = beta0 +  beta1 * time_variable + beta2 * group_variable + beta3 * time_variable * group_variable
    # z of dataset with beta3 = 0
    z_0 = beta0 +  beta1 * time_variable + beta2 * group_variable
    # z of dataset with beta2 = beta3 = 0
    z_00 = beta0 +  beta1 * time_variable
    
    # Random effects
    if (!with_rd_slope) {
      # Random intercepts for each subject
      tau = 0.8
      random_effects = rnorm(num_sub, mean = 0, sd = tau)
      random_intercept = random_effects[subject_id]
      random_slope = rep(0, num_sub*num_time)
    } else {
      # Random intercepts and slopes for each subject
      tau0 = 0.8
      tau1 = 0.4
      rho01 = 0.2
      sigma_mat = matrix(c(tau0^2, rho01*tau0*tau1, rho01*tau0*tau1, tau1^2), nrow = 2, ncol = 2, byrow = T)
      random_effects = mvrnorm(num_sub, mu = c(0, 0), Sigma = sigma_mat)
      random_intercept = random_effects[subject_id, 1]
      random_slope = random_effects[subject_id, 2]
    }
    
    # z of full dataset
    z_full = z_full + random_intercept + random_slope * time_variable
    # z of dataset with beta3 = 0
    z_0 = z_0 + random_intercept + random_slope * time_variable
    # z of dataset with beta2 = beta3 = 0
    z_00 = z_00 + random_intercept + random_slope * time_variable
    
    if (resp_dist == 'binomial') {
      # pi of full dataset
      pi_full = 1/(1+exp(-z_full))
      # pi of dataset with beta3 = 0
      pi_0 = 1/(1+exp(-z_0))
      # pi of dataset with beta2 = beta3 = 0
      pi_00 = 1/(1+exp(-z_00))
      
      # Dataset
      simulate_data_ls[[i]] = data.frame(
        y_full = rbinom(n_obs, size = 1, prob = pi_full),
        y_0 = rbinom(n_obs, size = 1, prob = pi_0),
        y_00 = rbinom(n_obs, size = 1, prob = pi_00),
        id = subject_id,
        time = time_variable,
        group = group_variable,
        random_intercept = random_intercept,
        random_slope = random_slope
      )
    }
    
    if (resp_dist == 'Poisson'){
      # mu of full dataset
      mu_full = exp(z_full)
      # mu of dataset with beta3 = 0
      mu_0 = exp(z_0)
      # mu of dataset with beta2 = beta3 = 0
      mu_00 = exp(z_00)
      
      # Dataset
      simulate_data_ls[[i]] = data.frame(
        y_full = rpois(n_obs, lambda = mu_full),
        y_0 = rpois(n_obs, lambda = mu_0),
        y_00 = rpois(n_obs, lambda = mu_00),
        id = subject_id,
        time = time_variable,
        group = group_variable,
        random_intercept = random_intercept,
        random_slope = random_slope
      )
    }
  }
  save(simulate_data_ls, 
       file = paste0(path, 'data/', 'n', num_sub, '_t', num_time, '_sim', 
                     num_sim, '_', resp_dist, ifelse(with_rd_slope, '_rdis', '_rdi'), '.RData'))
}

set.seed(126) # for reproducibility

n_sim = 1000
n_sub_v = c(20, 50, 200)
n_time_v = c(2, 5)
re_dist_v = c('binomial', 'Poisson')
with_rd_slope_v = c(FALSE, TRUE)
for (with_rd_slope in with_rd_slope_v){
  for (n_sub in n_sub_v){
    for (n_time in n_time_v){
      for (re_dist in re_dist_v){
        gen_data(num_sub = n_sub, num_time = n_time, num_sim = n_sim, resp_dist = re_dist, with_rd_slope = with_rd_slope)
      }
    }
  }
}

[paul-buerkner]

Thank you! I am sorry, I didn't reply so far. I will look at this issue in a few weeks when I have more time at my hands. In the meantime, if someone else has time to investigate this issue, I would be grateful.