#====================================================================== # Figure 12.4.R. # R script to calculate intercept b[1] and slope b[2] for the # least squares prediction line. Script produces a scatterplot # of the data along with overlayed prediction line. # # Data file should have two columns: response variable to be # predicted labeled "y" and predictor variable labeled "x". # Change the R working directory to the location of the data file. # Re-label axes in the plot( ) function if not using the Geyser data. #====================================================================== Geyser=read.table("old_faithful.txt",header=TRUE) # Change file name if # necessary. attach(Geyser) n=length(y) # Number of observations is n. X=matrix(1,n,2) # Form the X[,2]=x; # X matrix. b=solve((t(X)%*%X))%*%t(X)%*%y # Least squares estimates in b; # t() is transpose. plot(x,y,type="p",xlab="duration of eruption (minutes)", ylab="time until next eruption (minutes)") #Scatterplot of data. yhat1=b[1]+b[2]*min(x) # Calculate predicted y values at yhat2=b[1]+b[2]*max(x) # smallest and largest values of x. yhat=rbind(yhat1,yhat2) xvals=rbind(min(x),max(x)) points(xvals,yhat,type="l") # Connect two predicted values with line. "least squares intercept and slope:" # Print the intercept and slope b # to the console. detach(Geyser)