# Scott Merrill
# Oct 24 distribution to line


rm = list(ls = (all=TRUE))

sample.size = 40
data.error = rnorm(mean=0,sd = 2.5*sample.size^.5, n = sample.size)

# hidden slope line for compiling below
									slope.real = .5
# hidden slope line for compiling below
x = seq(1:sample.size)

y = x*slope.real + data.error
 plot(x,y, xlim = c(0,sample.size), lwd = 2)
abline(0,slope.real, col = 3, lwd = 3)

abline(0,slope.real, col = "white", lwd = 3)

mod.linear = lm(y~x)
summary(mod.linear)
abline(mod.linear$coefficients, lwd = 3, col = 4)

SSE = sum(mod.linear$residuals^2)
SSE

se.mod.linear = (SSE/(sample.size-1))^.5
se.mod.linear

df = sample.size - 2

se.1 = SSE # this is error in the y term
se.2 = sum((x-mean(x))^2) # variation in the x term - equivalent of SSE in x
se.slope = (1/df*(se.1/se.2))^.5 # basically SSE/(SSE in x) divided by the df
se.slope
abline(mod.linear$coefficients[1],mod.linear$coefficients[2]+ se.slope*1.96, col = 4, lwd = 3, lty = 3)
abline(mod.linear$coefficients[1],mod.linear$coefficients[2]- se.slope*1.96, col = 4, lwd = 3, lty = 3)

# intercept uncertainty
manual.se.intercept =-4.4710 
abline(mod.linear$coefficients[1]-manual.se.intercept,mod.linear$coefficients[2]+ se.slope*1.96, col = 4, lwd = 3, lty = 7)
abline(mod.linear$coefficients[1]+manual.se.intercept,mod.linear$coefficients[2]- se.slope*1.96, col = 4, lwd = 3, lty = 7)

poly.x = c(0,sample.size,sample.size)
poly.y = c(mod.linear$coefficients[1],(mod.linear$coefficients[1] + sample.size*(mod.linear$coefficients[2]+se.slope*1.96)),(mod.linear$coefficients[1] + sample.size*(mod.linear$coefficients[2]-se.slope*1.96)))
polygon(poly.x,poly.y, xpd=xpd, col="grey90", lty=2, lwd=2, border="grey95")
abline(mod.linear$coefficients, lwd = 3, col = 4)
points(x,y,lwd=2)

summary(mod.linear)