pitakakariki / simr

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

error "Couldn't find object's data." when using powerSim on lmer #274

Closed biancaserio closed 8 months ago

biancaserio commented 8 months ago

Hello,

I am having trouble using your powerSim() function applied to the following linear fixed effects model:


## for sex * total surface area interaction effect in lmer = Gradient_Eigenvalues ~ Sex + Age + tot_SA + Sex*tot_SA + random nested effect(family relatedness/twin status) 

family_id = HCP_demographics_cleaned_final$Family_ID
twin_status = HCP_demographics_cleaned_final$TwinStatus

gradient_val <- scale(HCP_hemi_array_aligned_fc_G1[[1]])
sex <- HCP_demographics_cleaned_final$Gender
age <- scale(HCP_demographics_cleaned_final$Age_in_Yrs)
totSA <- scale(HCP_demographics_cleaned_final$tot_SA) 

lmer_fit <- lmer(gradient_val ~ sex + age + totSA + sex*totSA +  (1 | family_id/twin_status), REML = FALSE)

I am specifically interested in finding out the power of the sex by tot_SA interaction, therefore;

fixef(lmer_fit)["sexM:totSA"] <- -0.01
power <- powerSim(lmer_fit, nsim = 200, test= fixed("sexM:totSA", method = "z"))
power

I get the following output, with 200 errors (in 200 simulations):


Power for predictor 'sexM:totSA', (95% confidence interval):
       0.00% ( 0.00,  1.83)

Test: z-test
      Effect size for sexM:totSA is -0.010

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

Time elapsed: 0 h 0 m 0 s

Furthermore, when checking lastResult()$errors, I get (for each simulation) "Couldn't find object's data.".

I would be extremely grateful if you could help me identify what I am doing wrong in specifying the model?

Please find below the output of the lmer (after manually changing the fixed effect size for "sexM:totSA"])


Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: gradient_val ~ sex + age + totSA + sex * totSA + (1 | family_id/twin_status)

     AIC      BIC   logLik deviance df.resid 
  2793.7   2832.9  -1388.8   2777.7      992 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1896 -0.5659 -0.1373  0.4019  5.1027 

Random effects:
 Groups                Name        Variance Std.Dev.
 twin_status:family_id (Intercept) 0.0726   0.2694  
 family_id             (Intercept) 0.2059   0.4538  
 Residual                          0.7083   0.8416  
Number of obs: 1000, groups:  twin_status:family_id, 647; family_id, 430

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)
(Intercept)   0.02204    0.05745 748.44186   0.384    0.701
sexM         -0.05878    0.08264 927.57203  -0.711    0.477
age           0.01184    0.03331 890.78772   0.356    0.722
totSA         0.08786    0.05793 913.12276   1.517    0.130
sexM:totSA   -0.01000    0.07969 964.13100  -0.125    0.900

Correlation of Fixed Effects:
           (Intr) sexM   age    totSA 
sexM       -0.639                     
age        -0.040  0.081              
totSA       0.563 -0.407  0.079       
sexM:totSA -0.378 -0.069 -0.035 -0.680

Thank you very much in advance.

Bianca

pitakakariki commented 8 months ago

Hi @biancaserio, just checking that you closed this because you solved the problem already?

biancaserio commented 8 months ago

Hi @pitakakariki, Yes - apologies for the confusion. I solved this problem by putting all the data in one dataframe (df) and then specifying data=df). Thanks for checking in! Have a nice day:)