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項指標ではないので、予測を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 

第三に、私は、正しい方法で行った場合、私はわからなかった部分であり、Iは指標マトリックスに応答変数になって、多変量を行いますこのような線形回帰：

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）の固有ベクトルの間違った計算につながったので、このコード部分を更新しました。第3に、応答変数を指標行列に変換し、トレーニングセットに対して多変量線形回帰を実行しました。そして、この線形モデルを使って予測する。

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 

第3の部分については、代わりに多項ロジスティック回帰を実行することもできます。しかし、このようにして、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 

したがって、ロジスティックモデルは線形モデルよりも優れたパフォーマンスを示します。