Hello all,
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 and the glmer 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
