Code modified from https://statmodeling.stat.columbia.edu/2017/02/11/measurement-error-replication-crisis/ https://statmodeling.stat.columbia.edu/wp-content/uploads/2017/02/graph-codes-to-share-for-science-paper-final.txt

```r

1. Load necessary library

library(svglite)

1. Set up SVG device

svg("three_plots_side_by_side.svg", width = 15, height = 5)

1. Set up 1x3 layout

par(mfrow = c(1, 3))

1. 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]

}

1. First plot

plot(sims[,2] ~ sims[,1], ylab="Observed with added error", xlab="Ideal Study") abline(0, 1, col="red")

1. Second plot

plot(sims[,4] ~ sims[,3], ylab="Observed with added error", xlab="Ideal Study") abline(0, 1, col="red")

1. third graph
2. 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)

}

1. Third plot

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

1. Close the SVG device

dev.off()