pitakakariki / simr

Power Analysis of Generalised Linear Mixed Models by Simulation
70 stars 20 forks source link

Error: invalid subscript type 'language' estimating power for glmerMod #137

Open JCruk opened 6 years ago

JCruk commented 6 years ago

I am trying to do an observed power calculation for a glmer binomial model; but, the results are 0% power with the same error on every simulation:

stage index message 1 Fitting 1 invalid subscript type 'language'

Example:

Data <- structure(list(Subject = structure(c(4L, 4L, 4L, 4L, 4L, 4L, 
7L, 7L, 7L, 7L, 7L, 7L, 2L, 2L, 2L, 2L, 2L, 2L, 6L, 6L, 6L, 6L, 
6L, 6L, 8L, 8L, 8L, 8L, 8L, 8L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("BEr1390", 
"CFe1144", "CWa2543", "JAn443", "JGu81", "JOb28", "LAu556", "SSc3625"
), class = "factor"), Condition = structure(c(3L, 3L, 2L, 2L, 
1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 
2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L, 3L, 3L, 2L, 2L, 1L, 1L
), .Label = c("WBOCL", "OROCL - TNR", "OROCL + TNR"), class = "factor"), 
    Prop = c(0.3, 0.4, 0.6, 0.4, 0.3, 0.5, 0.7, 0.5, 0.7, 0.6, 
    0.5, 0.7, 0.7, 0.6, 0.5, 0.6, 0.5, 0.7, 0.3, 0.3, 0.4, 0.6, 
    0.5, 0.3, 0.8, 0.7, 0.4, 0.8, 0.9, 0.9, 0.7, 0.7, 0.6, 0.7, 
    0.5, 0.6), n = c(10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
    10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
    10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
    10L, 10L, 10L, 10L)), row.names = c(1L, 2L, 3L, 4L, 5L, 6L, 
7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 25L, 
26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 
39L, 40L, 41L, 42L), class = "data.frame", na.action = structure(c(19L, 
20L, 21L, 22L, 23L, 24L, 43L, 44L, 45L, 46L, 47L, 48L), .Names = c("19", 
"20", "21", "22", "23", "24", "43", "44", "45", "46", "47", "48"
), class = "omit"), .Names = c("Subject", "Condition", "Prop", 
"n"))

Model:

Mod <- glmer(cbind(Prop*n,n - Prop*n) ~ Condition + (1|Subject), family = binomial, 
                          data = Data)

Model Summary:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(Prop * n, n - Prop * n) ~ Condition + (1 | Subject)
   Data: Data

     AIC      BIC   logLik deviance df.resid 
   137.7    144.0    -64.8    129.7       32 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.14108 -0.63739  0.02805  0.61375  1.33598 

Random effects:
 Groups  Name        Variance Std.Dev.
 Subject (Intercept) 0.2012   0.4486  
Number of obs: 36, groups:  Subject, 6

Fixed effects:
                       Estimate Std. Error z value Pr(>|z|)
(Intercept)           3.166e-01  2.633e-01   1.202    0.229
ConditionOROCL - TNR -3.948e-06  2.673e-01   0.000    1.000
ConditionOROCL + TNR -7.118e-02  2.667e-01  -0.267    0.790

Correlation of Fixed Effects:
            (Intr) CORO-T
CnOROCL-TNR -0.508       
CnOROCL+TNR -0.509  0.501

Power estimation: powerSim(Mod, test = fixed("Condition"), nsim = 100)

Result:

Power [user defined], (95% confidence interval):================================================================================================================================|
       0.00% ( 0.00,  3.62)

Test: [user defined function]

Based on 100 simulations, (0 warnings, 100 errors)
alpha = 0.05, nrow = 72

Time elapsed: 0 h 0 m 5 s

nb: result might be an observed power calculation
Warning message:
In observedPowerWarning(sim) :
  This appears to be an "observed power" calculation

List of errors:

lastResult()$errors

      stage index                           message
1   Fitting     1 invalid subscript type 'language'
2   Fitting     2 invalid subscript type 'language'
3   Fitting     3 invalid subscript type 'language'
4   Fitting     4 invalid subscript type 'language'
...
100 Fitting   100 invalid subscript type 'language'
pitakakariki commented 6 years ago

It's probably having trouble with the calculations within cbind. Try adding a variable for Prop*n to your data frame rather than calculating it in the model formula.

JCruk commented 6 years ago

That was the problem - thanks!!