# First just the original two plots, high power N = 3000, low power N = 50, true slope = .15 r <- .15 sims<-array(0,c(1000,4)) xerror <- 0.5 yerror<-0.5 for (i in 1:1000) { x <- rnorm(50,0,1) y <- r*x + rnorm(50,0,1) xx<-lm(y~x) sims[i,1]<-summary(xx)\$coefficients[2,1] x<-x + rnorm(50,0,xerror) y<-y + rnorm(50,0,yerror) xx<-lm(y~x) sims[i,2]<-summary(xx)\$coefficients[2,1] x <- rnorm(3000,0,1) y <- r*x + rnorm(3000,0,1) xx<-lm(y~x) sims[i,3]<-summary(xx)\$coefficients[2,1] x<-x + rnorm(3000,0,xerror) y<-y + rnorm(3000,0,yerror) xx<-lm(y~x) sims[i,4]<-summary(xx)\$coefficients[2,1] } plot(sims[,2] ~ sims[,1],ylab="Observed with added error",xlab="Ideal Study") abline(0,1,col="red") plot(sims[,4] ~ sims[,3],ylab="Observed with added error",xlab="Ideal Study") abline(0,1,col="red") # third graph # run 2000 regressions at points between N = 50 and N = 3050 r <- .15 propor <-numeric(31) powers<-seq(50,3050,100) xerror<-0.5 yerror<-0.5 for (j in 1:31) { sims<-array(0,c(1000,4)) for (i in 1:1000) { x <- rnorm(powers[j],0,1) y <- r*x + rnorm(powers[j],0,1) xx<-lm(y~x) sims[i,1:2]<-summary(xx)\$coefficients[2,1:2] x<-x + rnorm(powers[j],0,xerror) y<-y + rnorm(powers[j],0,yerror) xx<-lm(y~x) sims[i,3:4]<-summary(xx)\$coefficients[2,1:2] } # find significant observations (t test > 2) and then check proportion temp<-sims[abs(sims[,3]/sims[,4])> 2,] propor[j] <- table((abs(temp[,3]/temp[,4])> abs(temp[,1]/temp[,2])))[2]/length(temp[,1]) print(j) } plot(powers,propor,type="l",xlab="Sample Size",ylab="Prop where error slope greater",col="blue")