pitakakariki / simr

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

Errors in all simulations and 'nrow=NA' #143

Closed Hayafh closed 5 years ago

Hayafh commented 5 years ago

Hi, I am trying to do a power analysis for the following model:
BPD_FINAL ~ SEX + AGE + SES + ethnicity

I am trying to predict offspring BPD (BPD_FINAL) with parent substance use (DAD_DRUGS_SX) and seeing if the adoptive status (idabd) of the offspring moderate the association (DAD_DRUGS_SX*idabd). The data has sibling pairs, nested within families (IDYrFamf). I am trying to use SIMR to figure out what effect size I’m powered to detect. Here's what I have tried, and what I get:

library(simr) library(lme4) model1<- glmer(BPD_FINAL ~ SEX + AGE + SES + ethnicity

    • idabd + AGE_M + AGE_D +
  • DAD_DRUGS_SX + DAD_DRUGS_SXidabd +(1|IDYrFamf), data=zinb, na.action = na.exclude) Warning message: In glmer(BPD_FINAL ~ SEX + AGE + SES + ethnicity + idabd + AGE_M + : calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is deprecated; please call lmer() directly summary(model1) Linear mixed model fit by REML ['lmerMod'] Formula: BPD_FINAL ~ SEX + AGE + SES + ethnicity + idabd + AGE_M + AGE_D +
    DAD_DRUGS_SX + DAD_DRUGS_SX
    idabd + (1 | IDYrFamf) Data: zinb

REML criterion at convergence: 7403.1

Scaled residuals: Min 1Q Median 3Q Max -2.44186 -0.65810 -0.02018 0.59392 3.00327

Random effects: Groups Name Variance Std.Dev. IDYrFamf (Intercept) 13.11 3.620
Residual 72.40 8.509
Number of obs: 1019, groups: IDYrFamf, 520

Fixed effects: Estimate Std. Error t value (Intercept) 44.20132 4.56210 9.689 SEX -2.26402 0.59287 -3.819 AGE 0.02689 0.15163 0.177 SES 0.01861 0.72896 0.026 ethnicityKorean -1.82528 0.96996 -1.882 ethnicityOther -2.93594 1.71184 -1.715 idabd -3.24498 0.99987 -3.245 AGE_M -0.05144 0.11631 -0.442 AGE_D 0.04405 0.10230 0.431 DAD_DRUGS_SX 0.08398 0.50120 0.168 idabd:DAD_DRUGS_SX 1.26421 0.57741 2.189

powerSim(model1) Power for predictor 'SEX', (95% confidence interval):================================================================| 0.00% ( 0.00, 0.37)

Test: Kenward Roger (package pbkrtest) Effect size for SEX is -2.3

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

Time elapsed: 0 h 1 m 12 s


I have tried na.omit, and specifying a certain variable for fixed effects (fixef(model1)["SEX"] fixef(model1)["SEX"]<- -2.26402) but I sitll get the same error, and power of 0. What am I missing? :)