Подпрограмма Лапака dgesv: система точно единственная: U [6,6] = 0

Question

Я пытаюсь запустить приведенный ниже код, чтобы смоделировать набор P-значений с использованием обобщенной линейной модели

Однако я продолжаю получать ошибку : процедура Лапака dgesv: система точно единственное число: U [6,6] = 0

Вот код, который я пытаюсь запустить:

#which_p_value = "x1"
which_p_value = "groupcategory"
#which_p_value = "x1:groupcategory"

run_anova = FALSE 
simulate_mixed_effect = TRUE 
mixed_effect_sd = 3.094069
mixed_effect_sd_slope = 3.098661

library(tidyverse)
n_people <- c(2,5,10,15,20)
coef1 <- 1.61
coef2 <- -0.01
#coef3 <- 5
#coef4 <- 0

g1 = 0
g2 = 1
g3 = 2 
distances <- c(60,90,135,202.5,303.75,455.625)/100
n_trials <- 35
oneto1000 <- 25
n_track_lengths <- length(distances)
groupcategory = c(rep(g1, n_track_lengths), rep(g2, n_track_lengths),rep(g3,n_track_lengths))
z = c(n_people)
emptydataframeforpowerplots = NULL

coef3s <- c(-5, -4, -3, -2,-1, 0, 1, 2, 3, 4, 5)
coef4s <- c(-1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 1)

Datarray <- array(dim=c(length(coef3s), length(coef4s),length(n_people)))

coef3_counter =1
for (coef3 in coef3s) {
  coef4_counter =1
  for (coef4 in coef4s) {
    z1_g2 <- coef1 + coef2*distances + coef3*g2 + coef4*g2*distances
    z1_g3 <- coef1 + coef2*distances + coef3*g3 + coef4*g3*distances
    d = NULL
    pr1 = 1/(1+exp(-z1_g2))
    pr2 = 1/(1+exp(-z1_g3))

    counter=1
    for (i in n_people) {
      for (j in 1:oneto1000){
        df <- c()
        for (k in 1:i){

          # random effect from drawing a random intercept with sd = x

          if (simulate_mixed_effect){
            coef1_r = rnorm(1, mean=coef1, sd=mixed_effect_sd)
            coef2_r = rnorm(1, mean=coef1, sd=mixed_effect_sd_slope)
          } else {
            coef1_r = coef1
            coef2_r = coef2
          }

          z_g1 <- coef1_r + coef2*distances + coef3*g1 + coef4*g1*distances
          pr = 1/(1+exp(-z_g1))

          z1_g2 <- coef1_r + coef2*distances + coef3*g2 + coef4*g2*distances
          pr1 = 1/(1+exp(-z1_g2))

          if (run_anova) {
            df <- rbind(df, data.frame(x1 = c(rep(distances, 3)),
                                     y = c(rbinom(n_track_lengths,n_trials,pr), rbinom(n_track_lengths,n_trials,pr1),rbinom(n_track_lengths,n_trials,pr2)),
                                     groupcategory = groupcategory, id = c(rep(k,18))))

          } else { # this is for glmer data organisation
            for (m in 1:n_trials) {
            df <- rbind(df, data.frame(x1 = c(rep(distances, 3)),
 y = c(rbinom(n_track_lengths,1,pr),rbinom(n_track_lengths,1,pr1),rbinom(n_track_lengths,1,pr2)),groupcategory = groupcategory,id = c(rep(k,18))))

            }
          }
        }

        if (run_anova) {
          #df_aov <- aov(y~x1*groupcategory+Error(id/(x1*groupcategory)),data=df)
          #df_aov_sum <- summary(df_aov)
          #pvalue <- df_aov_sum[[5]][[1]][which_p_value,"Pr(>F)"]

          df_aov <- aov(y~x1*groupcategory+Error(id),data=df)
          df_aov_sum <- summary(df_aov)
          pvalue <- df_aov_sum[[2]][[1]][which_p_value, "Pr(>F)"]

        } else { # glmer
          mod_group_glmer <-  glmer(y ~ x1 + groupcategory + (1+x1|id), data = df, family = "binomial")
          sum <- summary(mod_group_glmer)
          pvalue <- sum$coefficients[which_p_value, "Pr(>|z|)"]
        }

        d = rbind(d,data.frame(pvalue))
      }
      count <- plyr::ldply(d,function(c) sum(c<=0.05))
      Datarray[coef3_counter,coef4_counter,counter] <- count$V1/oneto1000
      counter = counter +1
      d = NULL
    }
    coef4_counter = coef4_counter + 1
  }
  coef3_counter = coef3_counter + 1
}

Ниже приведен скрипт отладчика:

Lapack routine dgesv: system is exactly singular: U[6,6] = 0

8. stopifnot(length(value <- as.numeric(value)) == 1L)

7. nM$newf(fn(nM$xeval()))

6. (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf, n), control = list()) { n <- length(par) ...

5. do.call(optfun, arglist)

4. withCallingHandlers(do.call(optfun, arglist), warning = function(w) { curWarnings <<- append(curWarnings, list(w$message)) })

3. optwrap(optimizer, devfun, start, rho$lower, control = control, adj = adj, verbose = verbose, ...)

2. optimizeGlmer(devfun, optimizer = control$optimizer[[2]], restart_edge = control$restart_edge, boundary.tol = control$boundary.tol, control = control$optCtrl, start = start, nAGQ = nAGQ, verbose = verbose, stage = 2, calc.derivs = control$calc.derivs, use.last.params = control$use.last.params)

1. glmer(y ~ x1 + groupcategory + (1 + x1 id), data = df, family = "binomial")

Кто-нибудь сможет протянуть руку помощи относительно того, как я могу действовать отсюда?