首页 文章

当y是r中的指标矩阵时,如何进行多元线性回归?

提问于
浏览
0

这是我第一次发帖提问,希望看起来不会混淆 . 非常感谢你的时间 .

我正在开发一个zipcode数据集,可以在这里下载: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不是指标矩阵,所以我迷失了如何将它们解释回数字并将它们与实际测试数据进行比较以进一步计算预测误差 .

1 回答

  • 0

    我终于想出了如何用分类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
    

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

相关问题