# This code reproduces the figures in the Blogpost:
# http://jmbh.github.io/Deconstructing-ME/

# !!! DEFINE WHERE THE FIGURES SHOULD BE SAVED !!!
figDir <- '...'

# --------- 1) Simulation by the Authors ---------

set.seed(1)
nIter <- 1000
r <- .15
sims<-array(0,c(nIter,4))
xerror <- 0.5
yerror<-0.5

for (i in 1:nIter) {
  
  ## N = 50
  # No ME
  x <- rnorm(50,0,1)
  y <- r*x + rnorm(50,0,1)
  xx<-lm(y~x)
  sims[i,1]<-summary(xx)$coefficients[2,1]
  # ME
  x<-x + rnorm(50,0,xerror)
  y<-y + rnorm(50,0,yerror)
  xx<-lm(y~x)
  sims[i,2]<-summary(xx)$coefficients[2,1]
  
  ## N = 3000
  # No ME
  x <- rnorm(3000,0,1)
  y <- r*x + rnorm(3000,0,1)
  xx<-lm(y~x)
  sims[i,3]<-summary(xx)$coefficients[2,1]
  # ME
  x<-x + rnorm(3000,0,xerror)
  y<-y + rnorm(3000,0,yerror)
  xx<-lm(y~x)
  sims[i,4]<-summary(xx)$coefficients[2,1]
  
}


# --------- 2) Figure A: Densities of Sampling Distributions ---------

png(paste0(figDir, 'SamplingDistri.png'), width = 600, height = 400)

plot.new()
par(mar=c(4,2,1,1))
plot.window(xlim= c(-.3, .45), ylim=c(0,24))
box()
D_high <- density(sims[,2])
D_low <- density(sims[,4])
lines(D_high$x, D_high$y, lty = 2, lwd = 2)
lines(D_low$x, D_low$y, lty = 1, lwd = 2)
title(xlab = 'Coefficient Estimate', line = 2.5)
title(ylab = 'Density', line = 1)
axis(1, seq(-.40, .60, length = 11), cex.axis = 1)

# Mean of sampling distribution
abline(v = mean(sims[,2]), lty=2, col='red', lwd =2)
abline(v = mean(sims[,4]), lty=1, col='red', lwd = 2)
# True coefficient (without ME)
abline(v = .15, col='blue', lty=2, lwd = 2)
abline(v = mean(c(sims[,2], sims[,4])), col='green', lty=2, lwd = 2)

# Legend
legend(-.3, 22, c('Samp Dist: Low noise / high N', 'Samp Dist: High noise / small N'), lty = 1:2, lwd = c(2,2))
legend(-.3, 17, c('Mean Samp Dist: low noise', 'Mean Samp Dist: high noise'), lty = 1:2, col=c('red', 'red'), lwd=c(2,2))
legend(-.3, 12, c('True coefficient: no ME', 'True coefficient: ME'), lty = c(2,2), col=c('blue', 'green'), lwd=c(2,2))

dev.off()


# --------- 3) Figure B: Reproducing Figure in paper & rescaling ---------


png(paste0(figDir, 'ScalingIssue.png'), width = 700, height = 700)

par(mar=c(4,4,3,1), mfrow=c(2,2))

plot(sims[,4] ~ sims[,3], 
     xlab = 'No measurement error', 
     ylab = 'Measurement error', 
     main = 'Low noise / large N', 
     pch = 20, cex = .2, col = 'red',
     xlim=c(0.05, .2), 
     ylim=c(0.05, .2))
abline(0,1,col="black")
plot(sims[,2] ~ sims[,1], 
     xlab = 'No measurement error', 
     ylab = 'Measurement error', 
     main = 'High noise / small N', 
     pch = 20, cex = .2, col = 'red')
abline(0,1,col="black")

plot(sims[,4] ~ sims[,3], 
     xlab = 'No measurement error', 
     ylab = 'Measurement error', 
     main = 'Low noise / large N',
     xlim = c(-.4, .7), 
     ylim = c(-.4, .7), 
     pch = 20, cex = .2, col = 'red')
abline(0,1,col="black")
plot(sims[,2] ~ sims[,1], 
     xlab = 'No measurement error', 
     ylab = 'Measurement error', 
     main = 'High noise / small N',
     xlim = c(-.4, .7), 
     ylim = c(-.4, .7), 
     pch = 20, cex = .2, col = 'red')
abline(0,1,col="black")

dev.off()