# Show and tell on consequences of ignoring measurement error
####################################################################
# Model is   Y  = beta0 + beta1 X1 + beta2 X2 + epsilon
#            W1 = X1 + e1
#            W2 = W2 + e2
# Fit naive model 
#            Y  = beta0 + beta1 W1 + beta2 W2 + epsilon
####################################################################


rm(list=ls())
options(scipen=999) # To supress scientific notation
# Defining rmvn = function(nn,mu,sigma)
source('https://www.utstat.toronto.edu/brunner/Rfunctions/rmvn.txt')

set.seed(9999); n = 200
# True parameter values
beta0 = 1; beta1=1; beta2=0; phi12 = 0.50
# And also Var(epsilon) = Var(e1) = Var(e2) = Var(X1) = Var(X2) = 1
# Phi is the correlation matrix of the latent x variables
Phi = rbind(c(1,phi12),
            c(phi12,1))

mereg = function(n)
{
# Simulate one data set
e1 = rnorm(n); e2 = rnorm(n); epsilon = rnorm(n)
X = rmvn(n,c(0,0),Phi)
X1 = X[,1]; X2 = X[,2]
Y = beta0 + beta1*X1 + beta2*X2 + epsilon
W1 = X1 + e1
W2 = X2 + e2
# Fit the naive model
fit = lm(Y~W1+W2)
mereg <- summary(fit)$coefficients
return(mereg[3,4]) # Just the p-value for H0: beta2=0
} # End of function mereg

mereg(200)
mereg(200)
mereg(200)
mereg(200)
mereg(20)
mereg(20)
mereg(20)