当 y 是 r 中的指示矩阵时，如何执行多元线性回归？

Question

这是我第一次发布问题，希望它看起来不会令人困惑。非常感谢您的宝贵时间。

我正在研究一个邮政编码数据集，可以在这里下载：http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/zip.train.gz http://statweb.stanford.edu/~tibs/ElemStatLearn/datasets/zip.test.gz

总的来说，我的目标是将主成分回归模型与训练数据集上的前 3 台 PC 拟合，因为这些响应变量是 2、3、5 和 8 的手写数字，然后使用测试数据进行预测.我的主要问题是在 X 矩阵上执行 PCA 后，我不确定我是否正确地进行了回归部分。我已将响应变量转换为 2487*4 指标矩阵，并希望拟合多元线性回归模型。但是预测结果不是二项式指标，所以我很困惑我应该如何将预测解释回原始响应变量，即预测为 2、3、5 或 8。还是我完全做了回归部分错误的？这是我的代码如下：

首先，我构建了那些响应变量等于 2、3、5 和 8 的子集：

zip_train <- read.table(gzfile("zip.train.gz")) 
zip_test <- read.table(gzfile("zip.test.gz"))
train <- data.frame(zip_train)
train_sub <- train[which(train$V1 == 2 | train$V1 == 3 | train$V1 == 5 | train$V1 == 8),]
test <- data.frame(zip_test)
test_sub <- test[which(test$V1 == 2 | test$V1 == 3 | test$V1 == 5 | test$V1 == 8),]    
xtrain <- train_sub[,-1]
xtest <- test_sub[,-1]
ytrain <- train_sub$V1
ytest <- test_sub$V1

其次，我将X矩阵居中，使用svd计算前3个主成分：

cxtrain <- scale(xtrain)
svd.xtrain <- svd(cxtrain)
cxtest <- scale(xtest)
svd.xtest <- svd(cxtest)

utrain.r3 <- svd.xtrain$u[,c(1:3)] # this is the u_r
vtrain.r3 <- svd.xtrain$v[,c(1:3)] # this is the v_r
dtrain.r3 <- svd.xtrain$d[c(1:3)]
Dtrain.r3 <- diag(x=dtrain.r3,ncol=3,nrow=3) # creat the diagonal matrix D with r=3
ztrain.r3 <- cxtrain %*% vtrain.r3 # this is the scores, the new components

utest.r3 <- svd.xtest$u[,c(1:3)] 
vtest.r3 <- svd.xtest$v[,c(1:3)] 
dtest.r3 <- svd.xtest$d[c(1:3)]
Dtest.r3 <- diag(x=dtest.r3,ncol=3,nrow=3) 
ztest.r3 <- cxtest %*% vtest.r3 

第三，这是我不确定我是否以正确的方式做的部分，我将响应变量转换为指标矩阵，并执行如下的多元线性回归：

ytrain.ind <-cbind(I(ytrain==2)*1,I(ytrain==3)*1,I(ytrain==5)*1,I(ytrain==8)*1)
ytest.ind <- cbind(I(ytest==2)*1,I(ytest==3)*1,I(ytest==5)*1,I(ytest==8)*1)

mydata <- data.frame(cbind(ztrain.r3,ytrain.ind))
model_train <- lm(cbind(X4,X5,X6,X7)~X1+X2+X3,data=mydata)
new <- data.frame(ztest.r3)
pred <- predict(model_train,newdata=new)

但是，pred 不是指标矩阵，所以我迷失了如何将它们解释回数字并将它们与真实测试数据进行比较以进一步计算预测误差。

Answer 1

我终于弄清楚了如何使用分类 y 执行多元线性回归。首先，我们需要将 y 转换为指示矩阵，然后我们可以将该矩阵中的 0 和 1 解释为概率。然后在x上回归y建立一个线性模型，最后用这个线性模型用x的测试集进行预测。结果是一个与我们的指标矩阵具有相同维度的矩阵。所有条目也应该被解释为概率，尽管它们可能大于 1 或小于 0（这就是为什么它之前让我感到困惑）。所以我们需要找到每行的最大数量，看看哪个预测的 y 具有最高的概率，这个 y 将是我们的最终预测。这样，我们就可以将连续的数字转换回类别，然后制作一个表格与 y 的测试集进行比较。所以我更新了我以前的代码如下。

首先，我构建了那些响应变量等于 2、3、5 和 8 的子集（代码与我在问题中发布的代码保持不变）：

zip_train <- read.table(gzfile("zip.train.gz")) 
zip_test <- read.table(gzfile("zip.test.gz"))
train <- data.frame(zip_train)
train_sub <- train[which(train$V1 == 2 | train$V1 == 3 | train$V1 == 5 | train$V1 == 8),]
test <- data.frame(zip_test)
test_sub <- test[which(test$V1 == 2 | test$V1 == 3 | test$V1 == 5 | test$V1 == 8),]    
xtrain <- train_sub[,-1]
xtest <- test_sub[,-1]
ytrain <- train_sub$V1
ytest <- test_sub$V1

其次，我将 X 矩阵居中，并使用 eigen() 计算了前 3 个主成分。我更新了这部分代码，因为我在之前的代码中标准化了 x 而不是将其居中，导致 x 的协方差矩阵和 cov(x) 的特征向量的计算错误。

cxtrain <- scale(xtrain, center = TRUE, scale = FALSE) 
eigenxtrain <- eigen(t(cxtrain) %*% cxtrain / (nrow(cxtrain) -1)) # same as get eigen(cov(xtrain)), because I have already centered x before
cxtest <- scale(xtest, center = TRUE, scale = FALSE)
eigenxtest <- eigen(t(cxtest) %*% cxtest/ (nrow(cxtest) -1))
r=3 # set r=3 to get top 3 principles
vtrain <- eigenxtrain$vectors[,c(1:r)] 
ztrain <- scale(xtrain) %*% vtrain # this is the scores, the new componenets
vtest <- eigenxtrain$vectors[,c(1:r)] 
ztest <- scale(xtest) %*% vtest

第三，我将响应变量转化为指标矩阵，对训练集进行多元线性回归。然后用这个线性模型进行预测。

ytrain.ind <- cbind(I(ytrain==2)*1,I(ytrain==3)*1,I(ytrain==5)*1,I(ytrain==8)*1)
ytest.ind <- cbind(I(ytest==2)*1,I(ytest==3)*1,I(ytest==5)*1,I(ytest==8)*1)

mydata <- data.frame(cbind(ztrain,ytrain.ind))
model_train <- lm(cbind(X4,X5,X6,X7)~X1+X2+X3,data=mydata)
new <- data.frame(ztest)
pred<- predict(model_train,newdata=new)

pred 是一个包含所有概率条目的矩阵，因此我们需要将其转换回分类 y 列表。

pred.ind <- matrix(rep(0,690*4),nrow=690,ncol=4) # build a matrix with the same dimensions as pred, and all the entries are 0.
for (i in 1:690){
  j=which.max(pred[i,]) # j is the column number of the highest probability per row
  pred.ind[i,j]=1 # we set 1 to the columns with highest probability per row, in this way, we could turn our pred matrix back into an indicator matrix
}

pred.col1=as.matrix(pred.ind[,1]*2) # first column are those predicted as digit 2
pred.col2=as.matrix(pred.ind[,2]*3)
pred.col3=as.matrix(pred.ind[,3]*5)
pred.col4=as.matrix(pred.ind[,4]*8)
pred.col5 <- cbind(pred.col1,pred.col2,pred.col3,pred.col4) 

pred.list <- NULL
for (i in 1:690){
  pred.list[i]=max(pred.col5[i,])
} # In this way, we could finally get a list with categorical y

tt=table(pred.list,ytest)
err=(sum(tt)-sum(diag(tt)))/sum(tt) # error rate was 0.3289855

对于第三部分，我们还可以执行多项逻辑回归。但这样一来，我们就不需要将 y 转换为指示矩阵，我们只需将其分解即可。所以代码如下：

library(nnet)
trainmodel <- data.frame(cbind(ztrain, ytrain))
mul <- multinom(factor(ytrain) ~., data=trainmodel) 
new <- as.matrix(ztest)
colnames(new) <- colnames(trainmodel)[1:r]
predict<- predict(mul,new)
tt=table(predict,ytest)
err=(sum(tt)-sum(diag(tt)))/sum(tt) # error rate was 0.2627907

因此表明逻辑模型确实比线性模型表现更好。