{

#Suppose your data is in X and Y
#X is the vector/matrix of indep. variables
#Y is the vector of dep. variables.
test=matrix(scan("Avg-SpringFlows-April1SWE-49-04.txt", skip = 1), ncol=2, byrow=T)

Y=test[,1]		#Spring Flows
X=test[,2]		#April 1 SWE
# do a linear regression..
zz=lsfit(X,Y)


#get the residuals of the model..
x11=residuals(zz)

# residuals = Y - Yestimate ==> Yestimate = Y - residuals
yest=Y-x11

plot(Y, yest)
abline(a=0, b=1)

# model diagnostics..
# plot the histogram of the residuals and add the pdf. This is to
# see how good the residuals fit a normal distribution..
# you can do a KS test on this, or a qqplot..)
par(mfrow=c(2,3))
hist(x11,probability=T)
library(sm)

sm.density(x11,add=T)

qqnorm(x11)
qqline(x11) # Create a qq line
#plot the residuals against the X value
plot(X,x11,xlab="X",ylab="residuals")

# plot the Y estimates agains the residuals..
plot(yest,x11,xlab="estiimate of Y", ylab="residuals")

# plot the observed Y against the residuals
plot(Y, x11)

#plot the autocorrelation function - to see if the residuals are related
# to each other..
z1=acf(x11)


#----------------------------------------

#get the cross validated estimates..

n=length(Y)
yxval=1:n

index=1:n
for(i in 1:n){
index1=index[index != i]
Xval=X[index1]
Yval=Y[index1]

zz=lsfit(Xval,Yval)

#now estimate at the point that was dropped
xpred=c(1,X[i])
yest[i]=sum(zz$coef * xpred)

}
#now surface the x-validated estimates..
plot(Y,yest, xlim=range(Y,yest), ylim=range(Y,yest), xlab="True value", ylab="x-validated estimate")
lines(Y,Y)
}