Статистический тест, приводящий к NaNs и предупреждающий о сходимости

Question

Я немного новичок в R, поэтому примите это во внимание, посмотрев на код ниже. Итак, в принципе, у меня есть следующая функция. В основном он генерирует данные, а затем запускает статистический тест в самом конце, используя clmm2.

library(magrittr)
library(dplyr)  
library(data.table)
library(splitstackshape)
library(ordinal)
library(broom.mixed)
library(purrr)
library(ggplot2)

sim_experiment <- function ( n_sample = 1000,
                             n_sample20 = n_sample*20,
                             n_sample2 = n_sample*2,
                             prop_disp = 0.25,
                             prop_fem = 0.50,
                             disp_probability = 0.85,
                             nondisp_probability = 0.90,
                             fem_probability = 0.87,
                             mal_probability = 0.90) {

  id=factor(rep(1:n_sample, each=2))
  treat <- factor(rep(c("a","s"),times=2))
  datasocial <- data.frame(id,treat)
  datasocial$dispersal <- sample(c("n", "d"), length(levels(datasocial$id)), prob = c(1-prop_disp, prop_disp), T)[as.numeric(datasocial$id)]
  datasocial$sex <- sample(c("f","m"), length(levels(datasocial$id)), prob = c(prop_fem,1-prop_fem), T)[as.numeric(datasocial$id)]
  dataraw <- expandRows(datasocial, count=10, count.is.col=FALSE)
  for (i in 1:n_sample20) {
    dataraw$dispersal_prob[[i]] <- if (dataraw$dispersal[[i]]=="n") {
      nondisp_probability
    } else {disp_probability}
  }
  for (i in 1:n_sample20) {
    dataraw$sex_prob[[i]] <- if (dataraw$sex[[i]]=="f") {
      fem_probability
    } else {mal_probability}
  }
  dataraw$sample_asocial_prob1 = rep(c(1-0.22,1-0.23,1-0.25,1-0.26,1-0.27,1-0.29,1-0.30,1-0.31,1-0.32,1-0.33,0,0,0,0,0,0,0,0,0,0),times=n_sample)
  dataraw$sample_asocial_prob2 = rep(c(0.22,0.23,0.25,0.26,0.27,0.29,0.30,0.31,0.32,0.33,0,0,0,0,0,0,0,0,0,0),times=n_sample)
  for (i in 1:n_sample20) {
    dataraw$sample_social_prob1[[i]] <- if (dataraw$treat[[i]]=="a") {
      0
    } else {1 - dataraw$dispersal_prob[[i]]*dataraw$sex_prob[[i]]}
  }
  for (i in 1:n_sample20) {
    dataraw$sample_social_prob2[[i]] <- if (dataraw$treat[[i]]=="a") {
      0
    } else {dataraw$dispersal_prob[[i]]*dataraw$sex_prob[[i]]}
  }
  for (i in 1:n_sample20) {
    dataraw$risk_sensitivity[[i]] <- if (dataraw$treat[[i]]=="a") {
      sample(c(0, 1), prob = c(dataraw$sample_asocial_prob1[[i]],dataraw$sample_asocial_prob2[[i]]), T)}
    else {
      sample(c(0, 1), prob = c(dataraw$sample_social_prob1[[i]],dataraw$sample_social_prob2[[i]]), T)
    }}
  result <- aggregate(dataraw$risk_sensitivity,by=list(dataraw$id,dataraw$treat),sum)
  result <- result[order(result$Group.1) , ]
  datasocial$risk_sensitivity <- result$x 
  for (i in 1:n_sample2) {
    datasocial$risk_sensitivity[[i]] <- if (datasocial$treat[[i]]=="s" & datasocial$risk_sensitivity[[i]]==10) {
      9
    } else {datasocial$risk_sensitivity[[i]]}
  }
  for (i in 1:n_sample2) {
    datasocial$risk_sensitivity[[i]] <- if (datasocial$treat[[i]]=="s" & datasocial$risk_sensitivity[[i]]==0) {
      1
    } else {datasocial$risk_sensitivity[[i]]}
  }
  datasocial$risk_sensitivity <- as.factor(datasocial$risk_sensitivity)
  test <- clmm2(risk_sensitivity ~ treat + sex + dispersal + sex*dispersal + treat*dispersal + treat*sex,random = id, data = datasocial, Hess=TRUE)
  return(test)
}

Теперь я хочу провести анализ мощности, поэтому у меня есть еще один фрагмент кода.

tidy_output_clmm = function(fit){
  results = as.data.frame(coefficients(summary(fit)))
  colnames(results) = c("estimate","std.error","statistic","p.value")
  results %>% tibble::rownames_to_column("term")
}

sim_experiment_power <- function(rep) {
  s <- sim_experiment(n_sample = 1000,
                      prop_disp = 0.35,
                      prop_fem = 0.50,
                      disp_probability = 0.70,
                      nondisp_probability = 0.90,
                      fem_probability = 0.70,
                      mal_probability = 0.90)
  tidy_output_clmm(s) %>% mutate(rep=rep)
}

my_power <- map_df(1:10, sim_experiment_power)

ggplot(my_power, aes(estimate, color = term)) +
  geom_density() +
  facet_wrap(~term, scales = "free")

my_power %>% group_by(term) %>% summarise(power = mean(p.value < 0.05))

Как вы можете видеть, с этим кодом я в основном делаю 10 итерации одной и той же функции. Проблема в том, что иногда я получаю предупреждение «3: In sqrt (diag (v c)): произведено NaNs». И каждый раз я получаю предупреждение "clmm2 может не сойтись: оптимизатор 'ucminf' завершен с max |dient |: XXX". Мне вообще не понятно, почему это происходит, потому что, как я уже сказал, я немного новичок в R. У кого-нибудь есть представление о том, что происходит?