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]) (TSS <- sum((Mean_Jan - mean(Mean_Jan))^2)) (n <- length(Mean_Jan)) one_n <- as.matrix(rep(1,n),ncol=1) X <- cbind(ELEV_C,LAT,LONG) (xbar <- t(X) %*% one_n/n) XBAR <- cbind(xbar[1,1]*one_n, xbar[2,1]*one_n, xbar[3,1]*one_n) (SSCP_XX <- t(X-XBAR) %*% (X-XBAR)) (S_XX <- SSCP_XX/(n-1)) (S_XX_m12 <- solve(diag(sqrt(diag(S_XX))))) (R_XX <- S_XX_m12 %*% S_XX %*% S_XX_m12) (Lambda_X <- diag(eigen(R_XX)$values)) (Q <- eigen(R_XX)$vectors) Q %*% Lambda_X %*% t(Q) Z <- (X - XBAR) %*% S_XX_m12/sqrt(n-1) PC1 <- Z %*% Q[,1] PC2 <- Z %*% Q[,2] PC3 <- Z %*% Q[,3] PC <- cbind(PC1,PC2,PC3) plot(PC1,PC2, xlim=c(-.10,.20), ylim=c(-.10,.20)) W <- Z %*% Q Wstar <- cbind(one_n, W) (gammahat <- solve(t(Wstar) %*% Wstar) %*% t(Wstar) %*% Mean_Jan) (s2.Ws <- (t(Mean_Jan-Wstar%*%gammahat) %*% (Mean_Jan-Wstar%*%gammahat))/ (n-4)) P_Ws <- Wstar %*% solve(t(Wstar) %*% Wstar) %*% t(Wstar) P_Ws_1 <- Wstar[,-2] %*% solve(t(Wstar[,-2]) %*% Wstar[,-2]) %*% t(Wstar[,-2]) P_Ws_2 <- Wstar[,-3] %*% solve(t(Wstar[,-3]) %*% Wstar[,-3]) %*% t(Wstar[,-3]) P_Ws_3 <- Wstar[,-4] %*% solve(t(Wstar[,-4]) %*% Wstar[,-4]) %*% t(Wstar[,-4]) (SSR_W1 <- t(Mean_Jan) %*% (P_Ws - P_Ws_1) %*% Mean_Jan) (SSR_W2 <- t(Mean_Jan) %*% (P_Ws - P_Ws_2) %*% Mean_Jan) (SSR_W3 <- t(Mean_Jan) %*% (P_Ws - P_Ws_3) %*% Mean_Jan) (F_W1 <- SSR_W1/s2.Ws); (P_W1 <- 1-pf(F_W1,1,n-4)) (F_W2 <- SSR_W2/s2.Ws); (P_W2 <- 1-pf(F_W2,1,n-4)) (F_W3 <- SSR_W3/s2.Ws); (P_W3 <- 1-pf(F_W3,1,n-4)) Q_g <- Q[,1:2] W_g <- Z %*% Q_g Wstar_g <- cbind(one_n, W_g) (gammahat_g <- solve(t(Wstar_g) %*% Wstar_g) %*% t(Wstar_g) %*% Mean_Jan) (betahat_g_PC <- Q_g %*% gammahat_g[2:3]) (V_bh_g_PC <- s2.Ws[1,1] * (Q_g %*% solve(Lambda_X[1:2,1:2]) %*% t(Q_g))) Yhat_g <- Wstar_g %*% gammahat_g Yhat_b <- mean(Mean_Jan) + Z %*% betahat_g_PC cbind(Yhat_g, Yhat_b) plot(Mean_Jan ~ Yhat_b) abline(0,1, col="red")