pitakakariki / simr

Power Analysis of Generalised Linear Mixed Models by Simulation
69 stars 19 forks source link

Observed Power is 0 #136

Open djjohnson89 opened 5 years ago

djjohnson89 commented 5 years ago

I'm attempting to conduct an observed power calculation for the following model:

sM1 <- lmer(scale(strong) ~ scale(tHeight) + scale(tGrip) + scale(tBicep) + tRace +
  sex + race + scale(perStrong) + scale(height) +
  (1|subject) + (1|target),  
  data = subset(df, tSex == "male"),
  REML = F)

Which returns the following result:

summary(sM1, correlation = F)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: scale(strong) ~ scale(tHeight) + scale(tGrip) + scale(tBicep) +  
   tRace + sex + race + scale(perStrong) + scale(height) + (1 |subject) + (1 | target)
   Data: subset(df, tSex == "male")

     AIC      BIC   logLik deviance df.resid 
 25205.1  25309.5 -12588.6  25177.1    12728 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.8324 -0.6097  0.0473  0.6548  6.1893 

Random effects:
 Groups   Name        Variance Std.Dev.
 subject  (Intercept) 0.1321   0.3635  
 target   (Intercept) 0.2577   0.5077  
 Residual             0.3774   0.6143  
Number of obs: 12742, groups:  subject, 287; target, 157

Fixed effects:
                   Estimate Std. Error t value
(Intercept)       0.0323998  0.0767219   0.422
scale(tHeight)   -0.0788836  0.0441752  -1.786
scale(tGrip)      0.1643316  0.0440668   3.729
scale(tBicep)     0.2899332  0.0458503   6.323
tRaceasian       -0.4149583  0.1159914  -3.577
tRaceblack        0.3488566  0.1008313   3.460
sexmale           0.0418346  0.0572277   0.731
raceasian        -0.0414813  0.0542878  -0.764
raceblack        -0.1430989  0.0550271  -2.601
scale(perStrong) -0.0008812  0.0236644  -0.037
scale(height)    -0.0472586  0.0294134  -1.607

However, when I go to calculate the observed power for my data, I keep getting the following:

powerSim(sM1, test = fixed(xname = "tRace", method = c("z")), nsim = 3)

Power for predictor 'tRace', (95% confidence interval):
       0.00% ( 0.00, 70.76)

Test: z-test

Based on 3 simulations, (0 warnings, 3 errors)
alpha = 0.05, nrow = NA

Time elapsed: 0 h 0 m 27 s

nb: result might be an observed power calculation

This happens even when I greatly simplify the model:

sM2 <- lmer(strong ~ tRace +
  (1|target),
  data = subset(df, tSex == "male"),
  REML = F)

print(summary(sM2), correlation = F)

Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: strong ~ tRace + (1 | target)
   Data: subset(df, tSex == "male")

     AIC      BIC   logLik deviance df.resid 
 40408.0  40445.6 -20199.0  40398.0    13624 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.1330 -0.6100  0.0607  0.6932  3.4607 

Random effects:
 Groups   Name        Variance Std.Dev.
 target   (Intercept) 0.793    0.8905  
 Residual             1.081    1.0396  
Number of obs: 13629, groups:  target, 159

Fixed effects:
            Estimate Std. Error t value
(Intercept)   4.1262     0.1122   36.78
tRaceasian   -0.5666     0.1870   -3.03
tRaceblack    0.6995     0.1620    4.32

When I run powerSim() I get the same error for the simpler model:

powerSim(sM2, test = fixed(xname = "tRace", method = c("z")), nsim = 3)

Power for predictor 'tRace', (95% confidence interval):
       0.00% ( 0.00, 70.76)

Test: z-test

Based on 3 simulations, (0 warnings, 3 errors)
alpha = 0.05, nrow = 13629

Time elapsed: 0 h 0 m 59 s

nb: result might be an observed power calculation

In the full model, I get this error regardless of what predictor I am trying to measure observed power for. This clearly is not right. What am I missing here? Thanks.

djjohnson89 commented 5 years ago

The error messages I am getting are:

lastResult()$errors

       stage index                        message
1 Simulating     1 new levels detected in newdata
2 Simulating     2 new levels detected in newdata
3 Simulating     3 new levels detected in newdata

I believe this may refer to my predictors. tRace and race are contrast coded with Whites as the reference group and sex is effects coded (women = -.5, men = .5). How can I stop the simulations from creating new levels that don't exist?

After looking at some of the other problems logged in this forum, I'm not sure this is just a contrast issue. I am running rsim version 1.04. This problem occurs whether I have contrasts in my model (race and sex) or only a single continuous predictor.

pitakakariki commented 5 years ago

I usually see that error when there's missing data.

You could try e.g. mdata <- na.omit(subset(df, tSex == "male")) and see if that fixes things.

If that doesn't work, you might need to use droplevels instead of na.omit.

djjohnson89 commented 5 years ago

When running the very simple model:

sM1 <- lmer(strong ~ tRace +
  (1|target),
  data = na.omit(subset(df, tSex == "male")),
  REML = F)
powerSim(sM1, test = fixed(xname = "tRace", method = c("z")), nsim = 3)
lastResult()$errors

    stage index                 message
1 Testing     1 subscript out of bounds
2 Testing     2 subscript out of bounds
3 Testing     3 subscript out of bounds

I get the same result when I use drop levels instead of na.omit. This issue is happening even when I run a super simple model with only one continuous predictor (no factors):

sM2 <- lmer(strong ~ tBicep +
  (1|target),
  data = subset(df, tSex == "male"),
  REML = F)

When I run the powerSim() function I get the same errors. That is, I get a "new levels detected" when I use my data, and "subscript out of bounds" when I use na.omit().

pitakakariki commented 5 years ago

If you're using a z-test you'll need to specify the variable name from fixef(sM1), not just "tRace". It will lokk something like "tRaceabc".