bp_power.R

gistfile1.txt

out<-list()
for (s1 in seq(0,.2,by=.01)) {
    print(s1)
    test<-numeric()
    for (i in 1:250) {
        N=1000
        h=sqrt(.6)
        e=sqrt(1-h^2)
        m=.35
        s0=sqrt(1-m^2)
        Mod=runif(N,-1.75,1.75)
        PGS=rnorm(N,0,1)
        Env=rnorm(N,0,1)
        ystar=h*PGS+e*Env #Eq. 4b
        s=s0+s1*Mod #s part of Eq. 4a (this is how we will program it in Mplus, so I'm illustrating that here).
        y=m*Mod+s*ystar #Eq. 4a (slightly different notation from word doc: here the main effect is s0 and the moderating effect is s1)
        library(lmtest)
        test[i]<-bptest(y~Mod)$p.value
    }
    out[[as.character(s1)]]<-mean(ifelse(test<.05,1,0))
}

plot(names(out),unlist(out),xlab='s1',ylab='power',type='b',lwd=2)