# 2012.11.15 with help from OS # (classical) measurement error in X attenuates slopes # can we de-attenuate them? #uniform=TRUE uniform=FALSE if(uniform) CELSIUS = seq(-10,10,2) if(!uniform) CELSIUS = qnorm( (1:11)/12, 0, 6) F = 32+ 1.8 * CELSIUS n=length(CELSIUS) lm(F~CELSIUS) ErrorsInXModelOS=readLines(,14) ################ OS model { for (i in 1:n) { MU[i] <- INTERCEPT + SLOPE * CELSIUS[i] CELSIUS[i] ~ dnorm(mu.C, tau.C) celsius[i] ~ dnorm(CELSIUS[i], 0.25) F[i] ~ dnorm(MU[i], tau.F ) } INTERCEPT ~ dnorm( 0,0.001) SLOPE ~ dnorm( 0,0.001) mu.C ~ dnorm( 0,0.001) tau.C ~ dgamma( 0.01,0.01) tau.F ~ dgamma( 0.01,0.01) } writeLines(ErrorsInXModelOS,'ErrorsInXModelOS.txt') library('rjags') ####### sims INITS = list( SLOPE = 1, INTERCEPT = 35 ) SIMS=25 naive.intercept=rep(NA,SIMS) naive.slope = rep(NA,SIMS) os.intercept=rep(NA,SIMS) os.slope = rep(NA,SIMS) os.mu.C=rep(NA,SIMS) os.tau.C = rep(NA,SIMS) os.tau.F = rep(NA,SIMS) for (sim in 1:SIMS) { celsius = CELSIUS + round(rnorm(n,0,2) , 3) ; naive.fit = lm(F~celsius) naive.intercept[sim] = naive.fit$coefficients[1] naive.slope[sim] = naive.fit$coefficients[2] jags <- jags.model("ErrorsInXModelOS.txt", data = list(n=n,F=F,celsius=celsius), inits=INITS, n.chains = 1, n.adapt = 10000) update(jags, 10000) results1=jags.samples(jags, c('SLOPE','INTERCEPT','mu.C','tau.C','tau.F','CELSIUS'), 1000) os.intercept[sim] = apply(results1$INTERCEPT,1,mean) os.slope[sim] = apply(results1$SLOPE,1,mean) os.mu.C[sim] = apply(results1$mu.C, 1,mean) os.tau.C[sim] = apply(results1$tau.C,1,mean) os.tau.F[sim] = apply(results1$tau.F,1,mean) } par(mfrow = c(2, 2)) plot(naive.intercept,naive.slope,xlim=c(29,35),ylim=c(1.2,2.2),col="white") abline(h=1.8,v=32,lwd=1,col="blue") points(naive.intercept,naive.slope,pch=19,cex=0.25) plot(os.intercept,os.slope,xlim=c(29,35),ylim=c(1.2,2.2),col="white") abline(h=1.8,v=32,lwd=1,col="blue") points(os.intercept,os.slope,pch=19,cex=0.25,col="blue") plot(os.mu.C,1/sqrt(os.tau.C),cex=0.25) summary(1/sqrt(os.tau.F)) par(mfrow = c(1,1)) XLIM=c(min(min(celsius),min(CELSIUS)),max(max(celsius),max(CELSIUS))) plot(celsius,F,xlab="Temperature on Celsius scale", ylab="Fahrenheit", ylim=c(0.7*min(F),max(F)), pch=19, cex=1,xlim=XLIM) lines(celsius,naive.fit$fitted.values) s = sample(1:1000,100) ; ns=length(s) for (i in 1:n) { imputed = mean(results1$CELSIUS[i,,1]) points( results1$CELSIUS[i,s,1], rep(F[i],ns), col="red",pch=19, cex=0.10) points( c(imputed), c(F[i]), col="red",pch=19, cex=1) points(c(CELSIUS[i]), c(0.75*min(F)),col="red",pch=19, cex=2) segments(CELSIUS[i], 0.75*min(F), celsius[i], 0.85*min(F),pch=19,cex=0.5) points(c(celsius[i]), c(0.85*min(F)),pch=19, cex=1) segments(imputed, 0.70*min(F), celsius[i], 0.85*min(F), col="blue",pch=19,cex=0.75, lty="dotted") points(imputed, 0.70*min(F),pch=19, cex=1, col="red") segments(imputed, 0.70*min(F), CELSIUS[i], 0.75*min(F),col="red") } lines(CELSIUS,os.intercept[sim] + os.slope[sim]*CELSIUS,col="red") text(XLIM[1],max(F), paste("Large Red Circles = TRUE CELSIUS Values (uniform=",toString(uniform),")" ), adj=c(0,1),col="red",cex=0.8) text(XLIM[1],0.95*max(F),"Larger Black Circles = recorded C values", adj=c(0,1),cex=0.8) text(XLIM[1],0.9*max(F),"Smaller red circles = ave. posterior C values", adj=c(0,1),cex=0.8,col="red") text(XLIM[1],0.85*max(F),"tiny red dots = posterior C values", adj=c(0,1),cex=0.65,col="red") ##########