data=read.csv("http://chase.mathstat.dal.ca/~bsmith/stat3340/Data/cement.csv")
  #read data
 X=as.matrix(cbind(rep(1,13),data[,-c(1,2)])) #construct the X matrix
 y=matrix(data[,2],byrow=T,ncol=1)   #construct Y
  n=dim(X)[1]   #number of observations
  k=dim(X)[2]-1 #number of predictor variables in full model
  r=1           #number of parameters set to 0 under the null hypothesis

beta=solve(t(X)%*%X)%*%t(X)%*%y  #calculate least squares estimates
SSE=t(y)%*%y-t(beta)%*%t(X)%*%y  #get the error SS for full model
SSE
MSE=SSE/(n-k-1)                  #MSE under full model
MSE

X2=X[,c(1:3,5)]   #construct the X matrix for the reduced model
beta2=solve(t(X2)%*%X2)%*%t(X2)%*%y  #calculate least squares estimates
                                     #for the reduced model
SSE2=t(y)%*%y-t(beta2)%*%t(X2)%*%y   #get the error SS for the reduced model
SSE2

Fobs=((SSE2-SSE)/r)/MSE    #observed value of test statistic
Fobs
pv=1- pf(Fobs,r,n-k-1)     #p-value for the test
pv

lm.out=lm(y~x_1+x_2+x_3+x_4,data=data)
lm.out2=lm(y~x_1+x_2+x_4,   data=data)
anova(lm.out,lm.out2)

r2=2          #number parameters set to 0 under H_0
X3=X[,c(1:3)] #construct the X matrix for the reduced model
beta3=solve(t(X3)%*%X3)%*%t(X3)%*%y  #calculate least squares estimates
                                     #for the reduced model
SSE3=t(y)%*%y-t(beta3)%*%t(X3)%*%y   #get the SSE for the reduced model
Fobs2=((SSE3-SSE)/r2)/MSE    #observed value of test statistic
pv2=1- pf(Fobs2,r2,n-k-1)    #p-value for the test

c(SSE3, Fobs2, pv2)
lm.out3=lm(y~x_1+x_2,   data=data)
anova(lm.out,lm.out3)