{

#Suppose your data is in X and Y
#X is the vector/matrix of indep. variables
#Y is the vector of dep. variables.

# 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

# 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..)
hist(x11,probability=T)          
library(sm)

sm.density(x11,add=T)                

#plot the residuals against the X value
plot(X[,i],x11,xlab="X",ylab="residuals") 

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

#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,1:2])
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)
}