## 6.6节例题


################## 例6.6.1 ##################
library(car)
collinear = read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.6.1.txt", header=TRUE)
apply(collinear,2,mean) #各列均值
apply(collinear,2,sd)  #各列标准差
RX= cor(collinear[2:7]); kappa(RX, exact=TRUE)
lm.fit = lm(Y~., data=collinear)
round(vif(lm.fit), 3)





################## 例6.6.2 ##################
library(MASS)
library(car)
dataF <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.6.2.txt", header = T)
RX<-round(cor(dataF[,-4]),3) ##自变量的相关系数矩阵, 保留三位小数
mu<-apply(dataF,2,mean)  #各列均值
s<-apply(dataF,2,sd)    #各列标准差
lm.fit = lm(y~., data=dataF)
round(vif(lm.fit), 3)  ##方差扩大因子
eigen(RX)  #特征根和特征向量 

##条件数
kappa(RX, exact=TRUE)  # 677.6465=eigen(RX)$values[1]/eigen(RX)$values[3]

round(s,3)   #标准差
round(mu,3) #均值

## 岭估计
data1=scale(dataF)##    X和y标准化
lm.ridge(y~0+x1+x2+x3, data=data.frame(data1), lambda=c(seq(0,0.01,0.001),
                                                        seq(0.02,0.1,0.01), seq(0.2, 1,0.1))) ### 岭估计表6.6.3
beta<-lm.ridge(y~0+x1+x2+x3, data=data.frame(data1), 
               lambda=c(seq(0,0.2,0.01),seq(0.3,1,0.1)  )) ### 岭估计表6.6.3
beta
M<-as.matrix(data1)
M1<-M[1:11,]
mode(data1)    ####查看类型得到是 list 型
beta_x<-t(beta[[1]]);beta_x ##提取第一个列表并转置得到 beta_x
dim(beta_x)   #查看矩阵维数
beta_x[,1:3]  #提取这3列

RSS<-numeric(dim(beta_x)[1])
A<-M1
B<-beta_x
for (i in 1:dim(B)[1]){
  RSS[i]<-t(A[,4] - A[,-c(4)]%*%B[i,])%*%(A[,4] - A[,-c(4)]%*%B[i,])
}
RSS
round(cbind(B,RSS),3)   # 表6.6.3



par(mfrow=c(1,1), ann=TRUE, mar = c(5, 5, 1, 1))
####用matplot画岭迹点图6.6.1
ridgesol<-lm.ridge(y~0+x1+x2+x3, data=data.frame(data1),
                   lambda=seq(0,0.8,0.001))
matplot(x=ridgesol$lambda, y=t(ridgesol$coef),lwd=3, type="l",lty=c(1,2,3),
        xlab = expression(italic("k")),
        ylab=expression(italic(paste(widehat(beta)[i](k)))))
text(.75,0.15,expression(italic(paste(widehat(beta)[2](k)))))
text(.75,0.35,expression(italic(paste(widehat(beta)[1](k)))))
text(.75,0.53,expression(italic(paste(widehat(beta)[3](k)))))

## 取k=0.4 
lm.ridge(y~0+x1+x2+x3, data=data.frame(data1), lambda=0.4)
beta<-c(0.416,0.213, 0.531)
s<-round(s,3)
ave<-round(mu,3)
sr<-c(s[4]/s[1], s[4]/s[2], s[4]/s[3])
betaorigin<-round(sr*beta,4)        #回归系数 0.0630 0.5869 0.1169 
ybar<-ave[4]
xbar<-ave[1:3]
beta0<-ybar-t(sr)%*% (xbar*beta) # [1,] -8.647336








################## 例6.6.3 ##################
# 主成分回归完整流程（含标准化）

# 1. 数据读取与预处理
dataF <- read.table("/Users/macbook/Desktop/线性模型例题/Chapter6/data/例6.6.2.txt", header = T)
data1 <- scale(dataF)  # 数据标准化
RX <- round(cor(dataF[,-4]),3)  # 自变量相关系数矩阵
mu <- apply(dataF,2,mean)  # 各列均值
s <- apply(dataF,2,sd)    # 各列标准差

# 2. 主成分分析
X <- as.matrix(data1[,-4])  # 提取标准化自变量
pca_result <- prcomp(X, scale. = FALSE)  # 主成分分析
pca_scores <- pca_result$x  # 主成分得分

# 方差贡献分析
pca_sd <- pca_result$sdev  # 主成分标准差
pca_var <- pca_sd^2  # 特征值
pca_prop <- pca_var / sum(pca_var)  # 方差贡献率
pca_cumprop <- cumsum(pca_prop)  # 累计方差贡献率（前2个达99.91%）

# 主成分载荷
pca_loadings <- pca_result$rotation
round(pca_loadings, 4)  # 载荷矩阵

# 3. 主成分回归建模
y <- dataF[, 4]  # 因变量
pc1 <- pca_scores[, 1]  # 第一主成分得分
pc2 <- pca_scores[, 2]  # 第二主成分得分
X_pc <- data.frame(y, PC1 = pc1, PC2 = pc2)
pcr <- lm(y ~ PC1 + PC2, data = X_pc)  # 主成分回归模型
summary(pcr)  # 模型结果

# 4. 还原为原始变量方程
b0_pcr <- 21.8909  # 截距
b1_pcr <- 3.1350   # PC1系数
b2_pcr <- 0.8692   # PC2系数

# 提取载荷
loadings_pc1 <- pca_loadings[, "PC1"]
loadings_pc2 <- pca_loadings[, "PC2"]

# 计算标准化系数
beta_x1 <- b1_pcr * loadings_pc1["x1"] + b2_pcr * loadings_pc2["x1"]
beta_x2 <- b1_pcr * loadings_pc1["x2"] + b2_pcr * loadings_pc2["x2"]
beta_x3 <- b1_pcr * loadings_pc1["x3"] + b2_pcr * loadings_pc2["x3"]

# 计算原始变量方程参数
b0_original <- b0_pcr - (beta_x1 * mu["x1"]/s["x1"] + 
                           beta_x2 * mu["x2"]/s["x2"] + 
                           beta_x3 * mu["x3"]/s["x3"])

# 输出最终回归方程
cat("\n原始变量的回归方程：\n")
cat(sprintf("y = %.4f + %.4f*x1 + %.4f*x2 + %.4f*x3\n",
            b0_original,
            beta_x1 / s["x1"],
            beta_x2 / s["x2"],
            beta_x3 / s["x3"]))
# y = -9.1296 + 0.0728*x1 + 0.6092*x2 + 0.1063*x3


