Effect of measurement error mk2

effect of measurement error3.R

N=10000 #iterations
obs=20 #observations per iteration
beta=exp(rnorm(N,-2,1))

Bs=c() #bias in estimated Beta (standardised) for no error condition
Bsm=c() #bias in estimated Beta (standardised) for with error condition
betas2=c() #unstandardised true beta for no error iterations that came out significant
betasm2=c() #unstandardised true beta for with error iterations that came out significant

for(i in 1:N){
  x=rnorm(obs,0,1) #iv
b=rnorm(obs,0,1) #random variation in dv not due to iv
s=b+beta[i]*x #score of dv without measurement error
sm=b+rnorm(obs,0,1)+beta[i]*x #score of dv with measurement error
summs=summary(lm(scale(s)~scale(x))) #lm summary of no error regression
summsm=summary(lm(scale(sm)~scale(x)))#lm summary of with error regression
Betas=beta[i]/sd(s) #standardised true Beta for no error regression
Betasm=beta[i]/sd(sm) #standardised true Beta for with error regression

#select significant...
if(summs$coefficients[2,4]<.05) {
  Bs=c(Bs,summs$coefficients[2,1]-Betas) #bias of estimated Beta in no error regression
  betas2=c(betas2,beta[i]) #save current beta for later true t distribution
}
if(summsm$coefficients[2,4]<.05) {
  Bsm=c(Bsm,summsm$coefficients[2,1]-Betasm) #bias of estimated Beta in with error regression
  betasm2=c(betasm2,beta[i]) #save current beta for later true t distribution
}
}

#results
mean(Bs) #average bias of sig results without error
mean(Bsm)#average bias of sig results with error
par(mfrow=c(1,2))
plot(density(beta))
plot(density(Bs),main='bias in Beta')
points(density(Bsm),type='l',col='red')
legend(x='topleft',legend=c('no error', 'error'),text.col=c(1,2))


plot(betas2,Bs,pch=16,col=rgb(1,0,0,.2),xlab='True Beta',ylab='Bias in Beta')
points(betasm2,Bsm,pch=16,col=rgb(0,1,0,.2))
legend(x='topleft',legend=c('no error', 'error'),text.col=c(2,3))