## Use simulation to illustrate the effects of covariate measurement
## error on the sampling variance and bias of beta_hat in a simple
## linear regression.

## Sample sizes.
N = seq(50, 500, by=50)

## The measurement error standard deviation.
Esd = c(0, 0.5, 1)

## Storage for bias and variance.
QB = array(0, c(length(N),length(Esd)))
QV = array(0, c(length(N),length(Esd)))

for (k in 1:length(N)) {
    for (j in 1:length(Esd)) {

        ## The sample size.
        n = N[k]

        ## Replications.
        B = NULL
        for (r in 1:100) {

          ## Simulate the exact covariate and the response.
          X = rnorm(n)
          Y = X + rnorm(n)

          ## Simulate the observed covariate.
          Xobs = X + Esd[j]*rnorm(n)

          ## Fit the model and store the slope estimate.
          m = lm(Y ~ Xobs)
          B[r] = m$coeff[2]
        }

    QB[k,j] = (mean(B) - 1)^2
    QV[k,j] = var(B)
  }
} 

## Plot the squared bias curves in blue.
ymax = 1.5*max(QB[,2:3])
par(mar=c(5,4,4,4))
plot(N, QB[,1], ylim=c(0, ymax), xlab='Sample size', ylab="", 
     col='blue', t='b', pch=1)
for (j in 2:length(Esd)) {
    points(N, QB[,j], t='b', col='blue', pch=j)
}

## Plot the variance curves in red.
ymax = 1.5*max(QV[,2:3])
par(new=TRUE)
plot(N, QV[,1], ylim=c(0, ymax), xlab='', ylab='', axes=F, col='red', 
     t='b', pch=1)
for (j in 2:length(Esd)) {
    points(N, QV[,j], t='b', col='red', pch=j)
}
axis(4)
par(new=FALSE)

mtext("Variance",side=4,line=2,col='red')
mtext("Squared bias",side=2,line=2,col='blue')

legend('topright', c('Low ME', 'Medium ME', 'High ME', "Squared bias", 
       "Variance"), col=c('blue', 'blue', 'blue', 'blue', 'red'), 
       lty=c(-1,-1,-1,1,1), pch=c(1,2,3,-1,-1), ncol=2)