pitakakariki / simr

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

Error estimating power for detection of interaction effect #160

Open jonschaefer89 opened 5 years ago

jonschaefer89 commented 5 years ago

Hi there!

Like others in this forum, I'm also trying to estimate power for an interaction term. My data are drawn from a large, longitudinal twin study.

This is my model: model2<-lmer(PBF_E~polyvctzce18*zrpgsSCZ+(1|familyid), data=Erisk) summary(model2)

PBF_E represents a continuous measure of psychopathology standardized to mean=100, SD=15. polyvctzce18 represents a measure of victimization exposure with four levels (0,1,2,3+), but which I'm treating as continuous for the moment zrpgsSCZ represents a continuous measure of genetic risk for psychopathology familyid is the number assigned to each twin pair participating in the study. Each individual participant (twin) is nested within a single twin pair (or single family) with a single co-twin.

The results of this model are as follows:

> `Linear mixed model fit by REML ['lmerMod']
> Formula: PBF_E ~ polyvctzce18 * zrpgsSCZ + (1 | familyid)
>    Data: Erisk
> 
> REML criterion at convergence: 14705.9
> 
> Scaled residuals: 
>     Min      1Q  Median      3Q     Max 
> -2.1785 -0.6499 -0.0830  0.5524  3.3487 
> 
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  familyid (Intercept)  53.46    7.312  
>  Residual             114.27   10.690  
> Number of obs: 1860, groups:  familyid, 949
> 
> Fixed effects:
>                       Estimate Std. Error t value
> (Intercept)            95.5601     0.3972 240.591
> polyvctzce18            7.6976     0.3289  23.407
> zrpgsSCZ                1.0134     0.3802   2.666
> polyvctzce18:zrpgsSCZ  -0.3772     0.3200  -1.179
> 
> Correlation of Fixed Effects:
>             (Intr) plyv18 zrpSCZ
> polyvctzc18 -0.499              
> zrpgsSCZ     0.052 -0.022       
> plyvc18:SCZ -0.025  0.026 -0.494`

I'm trying to conduct a power analysis to see if my original model might be underpowered to detect an interaction between victimization exposure (polyvctzce18) and genetic risk (zrpgsSCZ). To see what sized interaction I am sufficiently well-powered to detect, I've tried setting the effect size of the interaction to various levels using code from this forum.

For example: fixef(model2)["polyvctzce18:zrpgsSCZ"]<-5 powerSim(model2,test=fixed("polyvctzce18:zrpgsSCZ","z"), nsim=100)

My problem is that it seems no matter what value I set the effect size to, the estimated power is always 0.00 (0.00, 3.62), probably because my 100 simulations generate 100 errors.

`Power for predictor 'polyvctzce18:zrpgsSCZ', (95% confidence interval):============================| 0.00% ( 0.00, 3.62)

Test: z-test Effect size for polyvctzce18:zrpgsSCZ is 5.0

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

Time elapsed: 0 h 0 m 12 s`

Can anyone tell what I'm doing to cause these errors? I'm not much of a statistician or programmer, so any help or insight is much appreciated!

pitakakariki commented 5 years ago

You can view the errors with lastResult()$errors.