measurement_error_and_power__GxE.R

simfun<-function(N,a1,beta.inter=.1) {
    ff<-function(i,N,a1,a2,beta.inter) {
        add.error<-function(t,alpha) {
            ve<-var(t)/alpha-var(t)
            e<-rnorm(N,mean=0,sd=sqrt(ve))
            t+e
        }
        g<-rnorm(N)
        e<-rnorm(N)
        y<-.2*g+.5*e+beta.inter*g*e+rnorm(N)
        ge<-add.error(g,alpha=a1)
        ee<-add.error(e,alpha=a2)
        lm(y~ge*ee)->mod
        summary(mod)$coef[4,4]<.05->test
        #pow[i]
        test
    }
    ##
    co<-r2<-list()
    for (a1 in a1) {
        for (a2 in seq(.1,1,by=.1)) {
            print(a2)
            library(parallel)
            nsim<-list()
            for (i in 1:10000) i->nsim[[i]]
            mclapply(nsim,ff,N=N,a1=a1,a2=a2,mc.cores=30,mc.set.seed=TRUE,beta.inter=beta.inter)->pow
            c(a1,a2,mean(unlist(pow)))->co[[paste(a1,a2)]]
        }
    }
    do.call("rbind",co)->mat
    mat
}

L<-list()
for (a1 in c(.25,.5)) {
    mat<-list()
    for (N in c(500,1000,5000,10000)) {
        print(N)
        simfun(N=N,a1=a1,beta.inter=.1)->mat[[as.character(N)]]
    }
    mat->L[[as.character(a1)]]
}

par(mgp=c(2,1,0),mar=c(3,3,1,1),oma=rep(.5,4),mfrow=c(1,2))
cols<-colorRampPalette(c("blue", "red"))( length(L[[1]]) )
for (i in 1:length(L)) {
    names(L)[i]->a1
    plot(NULL,xlim=c(0,1),ylim=c(0,1),xlab="Alpha of environmental variable",ylab=paste("Power (for alpha=",a1," of PGS)",sep=""))
    for (j in 1:length(L[[i]]))  {
        L[[i]][[j]]->mat
        lines(mat[,2],mat[,3],type="b",lwd=2,pch=19,col=cols[j])
        }
    if (i==1) {
        text(.2,.95,"True Model",cex=.75)
        text(.2,.9,"y=.2*pgs+.5*env+.1*pgs*env+e",cex=.75)
        text(.2,.85,"e~Normal[0,1]",cex=.75)
    }
    if (i==2)  legend("topleft",fill=rev(cols),rev(paste("N=",names(L[[1]]))),bty="n",cex=.75)
    ##
    abline(h=.8,col="gray")
    col2rgb("gray")->cc
    rgb(cc[1],cc[2],cc[3],max=255,alpha=35)->cc
    polygon(c(.95,.95,1.05,1.05),c(-1,2,2,-1),col=cc)
}