### how errors in X affect slope of Y on X       jh bios601 2011.09.21


library(animation)

x = 0:25 
n=length(x)
y= 32 + (9/5)*x
min.y=min(y)
max.y=max(y)

x.min=-15; x.max=80 # limits for plotting

x.0=-10 ; dx=0.30   # for plotting b's
y.0=68; dy=4
B.hat=c(); I=0;

for (sigma.error in seq(0,22,0.1) ) {
    I=I+1;
    plot(c(-12,35), c(min.y,max.y), col="white",
    ylab="Degrees F", xlab="Degrees C", xlim=c(x.min,x.max))
    
    x.plus.error = x + rnorm(n,0,sigma.error)
    
    fit=lm(y~x.plus.error)
    b.hat=fit$coefficients[2]
    B.hat=c(B.hat,b.hat)
    segments(x.0,y.0,x.0+400*dx,y.0,lwd=0.5)
    segments(x.0,y.0+1.8*dy,x.0+400*dx,y.0+1.8*dy,col="green",lwd=0.5)
    text(x.0,y.0+1.8*dy,"1.8",adj=c(1,0.5), col="green") 
    text(x.0,y.0+0*dy,"0.0",adj=c(1,0.5)) 
    points(x,y,type="l", col="green", lwd=3)
    points(x,y, col="green", cex=1, pch=19)
    points(x.0+(1:I)*dx,y.0+B.hat*dy,cex=0.2,col="red")
    text(x.0+I*dx,y.0+b.hat*dy,
     paste("sigma(e_c)=",toString(sigma.error),", b=",toString(round(b.hat,2))),
       cex=0.75,col="red",adj=c(-0.25,0.5))
    points(x.plus.error,y,col="red")
    points(x.plus.error,fit$fitted.values,type="l", col="red")
    for (i in 1:n) segments(x[i],y[i],x.plus.error[i],y[i],col="red",lwd=0.5,lty=3)
    Sys.sleep(0.25)
}