Re: [問題] 平均VIF

看板Statistics作者celestialgod (攸藍)時間10年前 (2014/03/12 22:27)推噓3(3推 0噓 30→)

留言33則, 2人參與討論串2/2 (看更多)

※ 引述《anovachen (　)》之銘言： : 例如下列R程式碼模擬結果：我不太懂你模擬的意義在哪裏？我提供一個反例問你 library(car) n = 100; p = 5 X = matrix(NA, n, p) X[,1:(p-2)] = rnorm((p-2) * n) X[,p-1] = X[,1] + X[,2] + rnorm(n,0,0.5) X[,p] = X[,3] + X[,4] + rnorm(n,0,0.5) cor(X) # 一般而言，我們認為cor > 0.7，可能存在multicollinearity : [,1] [,2] [,3] [,4] [,5] : [1,] 1.00000000 -0.19604544 0.06225474 0.63444937 0.5073527 : [2,] -0.19604544 1.00000000 -0.06263743 0.54292224 0.3159969 : [3,] 0.06225474 -0.06263743 1.00000000 0.01226171 0.6073165 : [4,] 0.63444937 0.54292224 0.01226171 1.00000000 0.7277750 : [5,] 0.50735268 0.31599693 0.60731647 0.72777496 1.0000000 # 根據correlation matrix可知X4, X5還有 X1,X2,X4之間高相關 # 我們應該可以猜測存在multicollinearity Beta = rep(1, p) # 公平比較，因此五組變數的斜率設為一樣 Y = 3 + X %*% Beta + rnorm(n, 0, 0.1) dat = data.frame(Y, X) summary(lm(Y ~ ., dat))$coef : Coefficients: : Estimate Std. Error t value Pr(>|t|) : (Intercept) 3.01484 0.00872 345.75 <2e-16 *** : X1 0.97862 0.01928 50.76 <2e-16 *** : X2 0.99301 0.02074 47.89 <2e-16 *** : X3 0.96549 0.01675 57.64 <2e-16 *** : X4 0.97222 0.02425 40.09 <2e-16 *** : X5 1.03393 0.01442 71.71 <2e-16 *** # 由上表可知檢定結果都是顯著，而估計值都跟真實值靠近 (VIF = vif(lm(Y ~ ., dat))) : X1 X2 X3 X4 X5 : 5.526498 4.764776 4.258915 13.936465 9.189339 mean(VIF) : 7.535199 # 而VIF 跟平均VIF如上所示，確實有multicollinearity存在 # 但是係數的估計跟檢定都未受影響，請問multicollinearity到底有什麼影響？ # 自問自答把剩下的補完 # VIF_j可以解釋成其方根倍的se(beta_j)為你所看到的SE... # 上句話很繞口，重說一次，VIF_j的方根可以解釋為你實際得到的se與真實se的倍數 # 一般解決這種問題有數種方式，可以自行上網尋找 # 其中一種就是ridge regression，此處以ridge regression算出其se為何？ lambda = 10 ^ (seq(-3,1.3,length = 100)) # Given lambda MSE_CV_f = function(lamb, fold){ X_tr = cbind(1, X[index != fold,]) X_te = cbind(1, X[index == fold,]) Y_tr = Y[index != fold] Y_te = Y[index == fold] mean((Y_te - X_te %*% solve(t(X_tr) %*% X_tr + lamb * diag(ncol(X_tr)), t(X_tr) %*% Y_tr))^2) } # 用cross-validation去找最恰當的lambda fold = 10 index = sample(rep(1:fold, n/fold), n) MSE_CV = sapply(lambda, function(lamb){ mean(sapply(1:5, function(fold) MSE_CV_f(lamb, fold))) }) (lambda_hat = lambda[which.min(MSE_CV)]) : [1] 0.001 beta_ridge = solve(t(cbind(1, X)) %*% cbind(1, X) + lambda_hat * diag(ncol(cbind(1, X))), t(cbind(1, X)) %*% Y) sd_beta_ridge = diag(solve(t(cbind(1, X)) %*% cbind(1, X) + lambda_hat * diag(ncol(cbind(1, X)))) * sum((Y - cbind(1, X) %*% beta_ridge)^2/(n-2))) table_ridge = cbind(beta_ridge, sd_beta_ridge, beta_ridge / sd_beta_ridge) dimnames(table_ridge) = list(c("Intercept", paste0("X",1:5)),c("Coef.", "SE", "t")) table_ridge : Coef. SE t : Intercept 3.0148113 7.292708e-05 41340.077 : X1 0.9785612 3.565188e-04 2744.767 : X2 0.9929565 4.123900e-04 2407.809 : X3 0.9654641 2.690656e-04 3588.211 : X4 0.9722463 5.640260e-04 1723.762 : X5 1.0339477 1.994020e-04 5185.241 # 這樣看不明顯，我補上兩個SE的比值 summary(lm(Y ~ ., dat))$coef[,2] / sd_beta_ridge : (Intercept) X1 X2 X3 X4 X5 : 119.56581 54.08005 50.28349 62.24999 42.99675 72.31116 # 結論：multicollinearity影響se的大小，並非影響斜率的估計(其估計還是BLUE) # 這篇已經偏離原問題，此處僅是回應前篇模擬的結果與結論 # 順帶提到有關多重共線性的影響 -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 36.238.199.20

推