txjan <- read.csv("http://www.stat.ufl.edu/~winner/sta6167/TX_jantemp.csv")
attach(txjan); names(txjan)

plot(txjan[,2:5])
cor(txjan[,2:5])

### Matrix Form 
(n <- length(Mean_Jan))
Y <- Mean_Jan
X0 <- rep(1,n)
X1 <- ELEV_C
X2 <- LAT
X3 <- LONG

# Total (Uncorrected) Sum of Squares
I <- diag(n)
(TSSU <- t(Y) %*% I %*% Y)
# Correction for the Mean P0 = J/n  SS(Int)
P0 <- X0 %*% solve(t(X0) %*% X0) %*% t(X0)
(R0 <- t(Y) %*% P0 %*% Y)
# Total (Corrected) Sum of Squares
(TSS <- t(Y) %*% (I-P0) %*% Y)
# SS(ELEV | Int)
X01 <- cbind(X0,X1)
P01 <- X01 %*% solve(t(X01) %*% X01) %*% t(X01)
(R1.0 <- t(Y) %*% (P01 - P0) %*% Y)
# SS(LAT | Int, ELEV)
X012 <- cbind(X0,X1,X2)
P012 <- X012 %*% solve(t(X012) %*% X012) %*% t(X012)
(R2.01 <- t(Y) %*% (P012 - P01) %*% Y)
# SS(LONG | Int, ELEV, LAT)
X0123 <- cbind(X0,X1,X2,X3)
P0123 <- X0123 %*% solve(t(X0123) %*% X0123) %*% t(X0123)
(R3.012 <- t(Y) %*% (P0123 - P012) %*% Y)
# SS(LAT | Int)
X02 <- cbind(X0,X2)
P02 <- X02 %*% solve(t(X02) %*% X02) %*% t(X02)
(R2.0 <- t(Y) %*% (P02 - P0) %*% Y)
# SS(ELEV | Int, LAT)
(R1.02 <- t(Y) %*% (P012 - P02) %*% Y)
# SS(ELEV | Int, LAT, LONG)
X023 <- cbind(X02,X3)
P023 <- X023 %*% solve(t(X023) %*% X023) %*% t(X023)
(R1.023 <- t(Y) %*% (P0123 - P023) %*% Y)
# SS(LAT | Int, ELEV, LONG)
X013 <- cbind(X01,X3)
P013 <- X013 %*% solve(t(X013) %*% X013) %*% t(X013)
(R2.013 <- t(Y) %*% (P0123 - P013) %*% Y)
# SSE
(SSE <- t(Y) %*% (I - P0123) %*% Y)
(df_E <- n-ncol(X0123))
(MSE <- SSE / df_E)

## Sequential Sums of Squares, F-tests, and Coeff of Partial Determination
F1.0 <- R1.0/MSE; F2.01 <- R2.01/MSE; F3.012 <- R3.012/MSE
P1.0 <- 1-pf(F1.0,1,df_E)
P2.01 <- 1-pf(F2.01,1,df_E)
P3.012 <- 1-pf(F3.012,1,df_E)
Rsq1.0 <- R1.0/TSS
Rsq2.01 <- R2.01/(TSS-R1.0)
Rsq3.012 <- R3.012/(TSS-R1.0-R2.01)

SS_seq_123 <- rbind(R1.0,R2.01,R3.012)
F_seq_123  <- rbind(F1.0,F2.01,F3.012)
P_seq_123  <- rbind(P1.0,P2.01,P3.012)
R2_seq_123 <- rbind(Rsq1.0,Rsq2.01,Rsq3.012)

SS_seq_out <- cbind(SS_seq_123, F_seq_123, P_seq_123, R2_seq_123)
colnames(SS_seq_out) <- c("Sequential SS", "F", "Pr(>F)","Partial R2")
rownames(SS_seq_out) <- c("ELEV_C", "LAT", "LONG")
round(SS_seq_out,4)

## Partial Sums of Squares, F-tests, and Coeff of Partial Determination
F1.023 <- R1.023/MSE; F2.013 <- R2.013/MSE; F3.012 <- R3.012/MSE
P1.023 <- 1-pf(F1.023,1,df_E)
P2.013 <- 1-pf(F2.013,1,df_E)
P3.012 <- 1-pf(F3.012,1,df_E)
R12.0 <- t(Y) %*% (P012 - P0) %*% Y
R13.0 <- t(Y) %*% (P013 - P0) %*% Y 
R23.0 <- t(Y) %*% (P023 - P0) %*% Y
Rsq1.023 <- R1.023/(TSS-R23.0)
Rsq2.013 <- R2.013/(TSS-R13.0)
Rsq3.012 <- R3.012/(TSS-R12.0)

SS_par_123 <- rbind(R1.023,R2.013,R3.012)
F_par_123  <- rbind(F1.023,F2.013,F3.012)
P_par_123  <- rbind(P1.023,P2.013,P3.012)
R2_par_123 <- rbind(Rsq1.023,Rsq2.013,Rsq3.012)

SS_par_out <- cbind(SS_par_123, F_par_123, P_par_123, R2_par_123)
colnames(SS_par_out) <- c("Partial SS", "F", "Pr(>F)","Partial R2")
rownames(SS_par_out) <- c("ELEV_C", "LAT", "LONG")
round(SS_par_out,4)

## Sequential and Partial SS/F-tests using lm function

txjan.mod1 <- lm(Mean_Jan ~ ELEV_C + LAT + LONG)
summary(txjan.mod1)
anova(txjan.mod1)               ## Sequential SS / F-tests
drop1(txjan.mod1, test="F")     ## Partial SS / F-tests