Stack Exchange
Generalised Linear Model - Response error is constant; boundary (singular) fit; leading leading minor of order 4 is not positive definite
I am trying to run the code below in order to simulate P-values using a generalised linear model. (I am running a logistic regression for power analysis): 
#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
}

Everytime I try to run the code...
I am getting the error: Error: response is constant.
Additionally, the code seems to take a very long time to run before it then prints the error...
Does anybody know why?
The debugger gives me this script:
Lapack routine dgesv: system is exactly singular: U[6,6] = 0

Which from my understanding could either indicate collinearity or another problem such as when a random effect variance is estimated to be very near zero and thus very loosely to the data that it is not sufficiently informative to drag the estimate away from the zero starting value.
Additionally, I get the error: the leading minor of order 4 is not positive definite
And a message along the lines of:
not positive definite or contains NA values: falling back to var-cov estimated from RXModel is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

Here is a snippet of my data output: 
 df
         x1 y groupcategory id
1   0.60000 1             0  1
2   0.90000 1             0  1
3   1.35000 1             0  1
4   2.02500 1             0  1
5   3.03750 1             0  1
6   4.55625 1             0  1
7   0.60000 0             1  1
8   0.90000 0             1  1
9   1.35000 0             1  1
10  2.02500 0             1  1
11  3.03750 0             1  1
12  4.55625 0             1  1
13  0.60000 0             2  1
14  0.90000 0             2  1
15  1.35000 0             2  1
16  2.02500 0             2  1
17  3.03750 0             2  1
18  4.55625 0             2  1
19  0.60000 1             0  1
20  0.90000 0             0  1
21  1.35000 1             0  1
22  2.02500 1             0  1
23  3.03750 0             0  1
24  4.55625 1             0  1
25  0.60000 0             1  1
26  0.90000 0             1  1
27  1.35000 0             1  1
28  2.02500 0             1  1
29  3.03750 0             1  1
30  4.55625 0             1  1
31  0.60000 0             2  1
32  0.90000 0             2  1
33  1.35000 0             2  1
34  2.02500 0             2  1
35  3.03750 0             2  1
36  4.55625 0             2  1

mod_group_glmer <-  glmer(y ~ x1 + groupcategory + (1+x1|id), data = df, family = "binomial")

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: y ~ x1 + groupcategory + (1 + x1 | id)
   Data: df
      AIC       BIC    logLik  deviance  df.resid 
 3648.796  3693.445 -1818.398  3636.796     12594 
Random effects:
 Groups Name        Std.Dev. Corr 
 id     (Intercept) 4.0727        
        x1          0.3478   -0.95
Number of obs: 12600, groups:  id, 20
Fixed Effects:
  (Intercept)             x1  groupcategory  
       1.2671        -0.1408        -6.8005  


The error would imply that your the response variable in your model is all the same. You can't fit a model to tell the difference between outcomes when all your observed outcomes are the exact same.

@MrFlick the value of my response variable does differ which is why I'm confused. How should I proceed?

Probably by providing some data.

@Edward, I have edited the post above
