讓我先說我確實認為這個問題是一個編碼問題，而不是一個統計問題。它幾乎肯定會在 Stats.SE 關閉。
R 中的leaps包有一個用于模型選擇的有用函式，稱為regsubsets該函式，對于任何給定大小的模型，它可以找到產生最小殘差平方和的變數。現在我正在閱讀Julian Faraway所著的Linear Models with R, 2nd Ed.一書。在第 154-5 頁，他有一個使用 AIC 進行模型選擇的示例。重現示例的完整代碼如下運行：
data(state)
statedata = data.frame(state.x77, row.names=state.abb)
require(leaps)
b = regsubsets(Life.Exp~.,data=statedata)
rs = summary(b)
rs$which
AIC = 50*log(rs$rss/50) (2:8)*2
plot(AIC ~ I(1:7), ylab="AIC", xlab="Number of Predictors")
該rs$which命令生成regsubsets函式的輸出，并允許您在繪制 AIC 并找到最小化 AIC 的引數數量后選擇模型。但問題是：雖然鍵入的示例運行良好，但當我嘗試使用此代碼并將其調整為其他資料時，我遇到了陣列中元素數量錯誤的問題。例如：
require(faraway)
data(odor, package='faraway')
b=regsubsets(odor~temp gas pack
I(temp^2) I(gas^2) I(pack^2)
I(temp*gas) I(temp*pack) I(gas*pack),data=odor)
rs=summary(b)
rs$which
AIC=50*log(rs$rss/50) (2:10)*2
產生警告資訊：
Warning message:
In 50 * log(rs$rss/50) (2:10) * 2 :
longer object length is not a multiple of shorter object length
果然length(rs$rss)=8，但是length(2:10)=9。現在我需要做的是模型選擇，這意味著我真的應該為每個模型大小設定一個 RSS 值。但是如果我b$rss在 AIC 公式中選擇，它不適用于原始示例！
所以這是我的問題：summary()函式的輸出在做什么regsubsets()？RSS值的數量不僅不一樣，而且值本身也不一樣。
uj5u.com熱心網友回復：
regsubsets有一個nvmax引數來控制“要檢查的子集的最大大小”。默認情況下為 8。如果將其增加到 9 或更高，則您的代碼可以正常作業。
但請注意，50您的 AIC 公式中的 是樣本量（即 50 個州statedata）。所以對于你的第二個例子，這應該是nrow(odor), 所以15。
uj5u.com熱心網友回復：
好的，所以你知道幫助頁面regsubsets說
regsubsets回傳一個"regsubsets"不包含用戶可維修部件的類物件。它被設計為由summary.regsubsets.
你即將找出原因。
中的代碼regsubsets呼叫 Alan Miller 的 Fortran 77 代碼進行子集選擇。也就是說，它不是我寫的，它在 Fortran 77 中。我理解演算法。在 1996 年寫作時leaps（以及在 2017 年進行重大修改時），我花了足夠的時間閱讀代碼以了解變數在做什么，但regsubsets主要遵循代碼隨附的 Fortran 驅動程式的結構。
該物件的rss欄位regsubsets具有該名稱，因為它存盤了一個RSS在 Fortran 代碼中呼叫的變數。該變數不是最佳模型的殘差平方和。 RSS在設定階段，在任何子集選擇完成之前，由 subrouteSSLEAPS計算，注釋為“從 AS75.1 的正交歸約中計算部分殘差平方和”。也就是說，RSS描述了設計矩陣中從左到右擬合的沒有選擇的模型的 RSS：只有最左邊的變數的模型，然后是最左邊的兩個變數，依此類推。如果他們不打算閱讀 Fortran，那么沒有理由任何人都需要知道這一點，因此沒有記錄。
中的代碼summary.regsubsets從物件的組件中提取輸出中的殘差平方和，該$ress組件來自RESSFortran 代碼中的變數。這是一個陣列，其元素是size 的第 -th 個最佳模型的[i,j]殘差平方和。ji
所有模型標準都是從$ress的同一個回圈中計算出來的summary.regsubsets，可以將其編輯為：
for (i in ll$first:min(ll$last, ll$nvmax)) {
for (j in 1:nshow) {
vr <- ll$ress[i, j]/ll$nullrss
rssvec <- c(rssvec, ll$ress[i, j])
rsqvec <- c(rsqvec, 1 - vr)
adjr2vec <- c(adjr2vec, 1 - vr * n1/(n1 ll$intercept -
i))
cpvec <- c(cpvec, ll$ress[i, j]/sigma2 - (n1 ll$intercept -
2 * i))
bicvec <- c(bicvec, (n1 ll$intercept) * log(vr)
i * log(n1 ll$intercept))
}
}
cpvec為您提供與 AIC 相同的資訊，但如果您想要 AIC，則執行相同的回圈并計算它會很簡單。
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/460021.html
上一篇：字串長度與字串長度串列的逐行比較