#these commands illustrate the use of R as a calculator

#the specific commands carry out simple linear regression, including
#calculation of least squares estimators, construction of the ANOVA
#table, and CI for the intercept, for the mean of Y, and the prediction interval


#enter x and y co-ordinates
x=c(.02,.07,.11,.15)  
y=c(242,237,231,201)
#calculate sums of squares
SXX=sum((x-mean(x))^2)
SXY=sum((x-mean(x))*(y-mean(y)))
SYY=sum((y-mean(y))^2)
#calculate LSE
b1=SXY/SXX
b0=mean(y)-b1*mean(x)
#get predicted values, the residuals
yp=b0+b1*x
resids=y-yp
SSE=sum(resids^2)  #residual SS
SST=SYY            #total SS
SSR=SST-SSE        #regression SS
SS=c(SSR,SSE,SST)  #print the SS
n=length(y)        #how many data points
df=c(1,n-2,n-1)    #regression, error and total degrees of freedom
MS=SS/df           #calculate the mean squares
cbind(SS,df,MS)    #joint the SS, df and MS columns side by side

cat("Fobs = ", MS[1]/MS[2], "with 1 numerator and ",n-2, " denominator degrees of freedom", sep=" ")


#calculate the 95% confidence interval for the slope

MSE=SSE/(n-2)      #MSE
#upper .025'th percentile of a t with n-2 df.
t=qt(.975,n-2)  
t

SEb1=sqrt(MSE/SXX) #standard error of beta_1
c(b1-t*SEb1,b1+t*SEb1) #report the conidence interval

#calculate a 95% CI for the mean of Y when x=x0

x0=.10             #which x0 
muhat=b0+b1*x0     #estimate mean of Y at x=x0
muhat
SEmu=sqrt(MSE)*sqrt(1/n+(x0-mean(x))^2/SXX) #S.E. of muhat
SEmu
c(muhat-t*SEmu, muhat+t*SEmu)  #report the CI

#calculate a 95% prediction interval for Y at x=x0

SEmu=sqrt(MSE)*sqrt(1+1/n+(x0-mean(x))^2/SXX) #standard deviation of the prediction
c(muhat-t*SEmu, muhat+t*SEmu)    #report the prediction interval