setwd("~/AndrewFiles/research/philosophy")

# First just the original two plots, high power N = 3000, low power N = 50, true slope = .15

sims<-array(0,c(1000,4))
for (i in 1:1000) {
x <- rnorm(50,0,1)
y <- .15*x + rnorm(50,0,1)
sims[i,1]<-cor(x,y)
x<-x + rnorm(50,0,.5)
y<-y + rnorm(50,0,.5)
sims[i,2]<-cor(x,y)

x <- rnorm(3000,0,1)
y <- .15*x + rnorm(3000,0,1)
sims[i,3]<-cor(x,y)
x<-x + rnorm(3000,0,.5)
y<-y + rnorm(3000,0,.5)
sims[i,4]<-cor(x,y)

}

pdf("science_noise_fig1a.pdf", height=5, width=5)
par(mar=c(3,3,2,1), pty="s", mgp=c(1.7,.5,0), tck=-.01)
range <- range(sims[,3:4])
plot(range, range, ylab="Estimate with high measurement error",xlab="Estimate in ideal study", bty="l", type="n")
points(sims[,3], sims[,4], pty=19, cex=.5, col="red")
mtext("Large sample size", side=3, line=0, cex=1.2)
abline(0, 1)
dev.off()


pdf("science_noise_fig1b.pdf", height=5, width=5)
par(mar=c(3,3,2,1), pty="s", mgp=c(1.7,.5,0), tck=-.01)
range <- range(sims[,1:2])
plot(range, range, ylab="Estimate with high measurement error",xlab="Estimate in ideal study", bty="l", type="n")
points(sims[,1], sims[,2], pty=19, cex=.5, col="red")
mtext("Small sample size", side=3, line=0, cex=1.2)
abline(0, 1)
dev.off()

# But maybe preferable to show proportion of observations where the with error cor is greater
# than the without error cor as a function of sample size.  Perhaps something like:


propor <-numeric(120)
powers<-seq(25,3000,25)

for (j in 1:120)  {

sims<-array(0,c(1000,4))
for (i in 1:1000) {
x <- rnorm(powers[j],0,1)
y <- .15*x + rnorm(powers[j],0,1)
sims[i,1]<-cor(x,y)
x<-x + rnorm(powers[j],0,.5)
y<-y + rnorm(powers[j],0,.5)
sims[i,2]<-cor(x,y)

}

propor[j] <- table((abs(sims[,2])> abs(sims[,1])))[2]/1000

# Or could take the 100 largest observations and see proportion among those
# xx<-order(abs(sims[,2]),decreasing = T)
#  propor[j] <- table((abs(sims[xx[1:100],2])> abs(sims[xx[1:100],1])))[2]/100

propor[j] <- table((abs(sims[,2])> abs(sims[,1])))[2]/1000

}

plot(powers,propor,type="l",xlab="Sample Size",ylab="Prop where error cor greater",col="blue")