## To conduct the F-test whether we can take the coefficient of the
## quadratic term in the cars model to be zero

> plot(cars,xlim=c(0,25))
> X <- cars[,1]
> Y <- cars[,2]
> n <- length(X)
> Design0 <- X
> Design <- cbind(X,X^2)

> Mod0 <- lm(Y~Design0-1)            # Fits Y_i = beta x_i + epsilon_i
> Mod <- lm(Y~Design-1)     # Fits Y_i = beta x_i + gamma x_i^2 + epsilon_i
> p <- 2
> p0 <- 1

> abline(Mod0)
> x <- seq(0,30,length=151)
> lines(x,Mod$coef[[1]]*x+Mod$coef[[2]]*x^2,lty=3)

> P0 <- Design0 %*% solve(t(Design0)%*%Design0) %*% t(Design0)
> P <- Design %*% solve(t(Design)%*%Design) %*% t(Design)

> RSS0 <- sum((Y - P0%*%Y)^2)  # Residual sum of squares 
> RSS0
> RSS <- sum((Y - P%*%Y)^2)
> RSS

> sum((P%*%Y - P0%*%Y)^2)
> RSS0 - RSS                  # Why are these two values the same?

> F <- ((RSS0 - RSS)/(p-p0))/(RSS/(n-p)) 
> F                           
 
> qf(0.95,p-p0,n-p)            # The upper 5% point of F_{p-p0,n-p}
> 1 - pf(F,p-p0,n-p)           # The p-value of the F-statistic

## Of course, it is possible to find out all this information directly in R:

> anova(Mod0,Mod)