#x = read.table("obp_al.txt")
x=read.table("obp_nl.txt")
names=x[,2]
OBP2002 = x[,4]
OBP2003 = x[,6]
cor(OBP2002,OBP2003)
#par(mfrow=c(2,1))
plot(OBP2002, OBP2003)
text(x=OBP2002, y = OBP2003, labels=names)

mean(OBP2002); mean(OBP2003)
sd(OBP2002); sd(OBP2003)
plot(OBP2002, OBP2003)
abline(0,1)
z = lm(OBP2003 ~ OBP2002)
z$coef
abline(z$coef)
summary(z)
sd(OBP2003)
plot(OBP2002, OBP2003 - OBP2002)
abline(h=0)
abline(v=.34)

plot(OBP2002, OBP2003)
abline(z$coef)
zp = predict(z,se.fit=T)
points(OBP2002, zp$fit + 2*zp$se.fit, pch="+")
points(OBP2002, zp$fit - 2*zp$se.fit, pch="+")