2015-09-26 30 views
0

我試圖寫我自己的R函數類似於向前逐步選擇step,但不是使用AIC作爲選擇標準,我有幾個標準需要評估每個添加預測變量的時間。該模型的構造原理解釋如下。模型應從與因變量具有最高相關性的預測變量開始。然後根據新模型是否符合以下標準,每次添加另一個預測變量。基於預定義標準的多元線性模型逐步選擇

  1. 調整後的r2值必須增加1%以上;
  2. 新添加的變量和現有變量的係數必須爲正;
  3. 增加的變量必須是顯着的,即p值< 0.05。

重複此過程直到沒有其他變量滿足所有三個條件。我需要的輸出可能只是最終模型中所有預測變量的名稱,相應係數和最終模型的r2值。

我的示例數據(y是因變量,X1 - 5233是預測)

data = structure(list(y = c(23.6, 19.9, 40.7, 40.7, 40.7, 40.7, 40.2, 
41.7, 41.7, 28.8), x1 = c(0.1, 0, 0.3, 0.3, 0.3, 
0.3, 0.3, 0.3, 0.3, 0.1), x2 = c(0, 0.1, 0, 0, 0, 
0, 0, 0.1, 0.1, 0), x3 = c(2277.6, 3038.1, 7797.9, 7797.9, 
7797.9, 7797.9, 8392.2, 10127.2, 10127.2, 1799), x4 = c(34228.7, 
49815, 76917.1, 76917.1, 76917.1, 76917.1, 75981.4, 74881.1, 
74881.1, 56798.2), x5 = c(108786.5, 150465.5, 230397.1, 230397.1, 
230397.1, 230397.1, 239300.9, 238493.8, 238493.8, 188799.5), 
x6 = c(362.2, 198.2, 656.6, 656.6, 656.6, 656.6, 681, 
655.3, 655.3, 222.3)), .Names = c("y", "x1", 
"x2", "x3", "x4", "x5", "x6"), row.names = c(NA, 
10L), class = "data.frame") 

在我的模型選擇功能

modSel = function(data, var){ 
cor.result = cor(data[,var], df["y"]) #calculate correlation coeff for each variable against y 
max.cor = rownames(cor.result)[which.max(cor.result)] #identify the variable with max cor 
start.model = lm(as.formula(paste("y", max.cor, sep = "~")), data) 
if #my criteria?? 
else #??] 

第一次嘗試沒有編程多的背景下,我真不對於如何在未知的時間重複評估我的標準沒有任何意見。我意識到要達到這個目標可能需要相當多的編碼,但對於初學者,我很感激關於整個框架應該是什麼樣子的任何指導。

乾杯

+0

在測試所有變量之前算法會如何結束?似乎所有的解釋變量都需要進行測試。 – Whitebeard

回答

0

希望這可以幫助你開始

功能從算法運行算法

modSel <- function(data) { 
    # initial 
    cor.result <- cor(data$y, data[, -which(colnames(data) == "y")]) 
    vars.model <- colnames(cor.result)[which.max(cor.result)] 
    vars.remaining <- colnames(data)[!colnames(data) %in% c("y", max.cor)] 
    start.model <- lm(as.formula(paste("y", vars.model, sep = "~")), data) 
    adj.rsq <- summary(start.model)$adj.r.squared 

    # algorithm 
    for (var in vars.remaining) { 
    # model 
    vars.test <- paste(vars.model, var, sep="+") 
    fit <- lm(as.formula(paste("y", vars.test, sep="~")), data) 
    new.rsq <- summary(fit)$adj.r.squared 

    # check adj rsq 
    cond1 <- new.rsq > adj.rsq + .01 

    # check coefficients 
    cond2 <- coefficients(fit)[var] > 0 

    # new var significant 
    cf <- summary(fit)$coefficients[, 4] 
    cond3 <- cf[var] < .05 

    if (cond1 & cond2 & cond3) { 
     vars.model <- vars.test 
     adj.rsq <- new.rsq 
    } 
    } 
    lm(as.formula(paste("y", vars.model, sep="~")), data) 
} 

電話modSel回報的最佳模式

bestfit <- modSel(data) 

彙總型號爲

summary(bestfit) 
Call: 
lm(formula = as.formula(paste("y", vars.model, sep = "~")), data = data) 

Residuals: 
    Min  1Q Median  3Q  Max 
-0.8731 -0.3838 -0.3838 0.6273 1.5640 

Coefficients: 
      Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.331e+01 1.941e+00 6.856 0.000241 *** 
x1   5.417e+01 5.834e+00 9.285 3.48e-05 *** 
x4   1.498e-04 4.460e-05 3.359 0.012099 * 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9345 on 7 degrees of freedom 
Multiple R-squared: 0.9904, Adjusted R-squared: 0.9876 
F-statistic: 360.5 on 2 and 7 DF, p-value: 8.721e-08 
+0

謝謝,這工作出色。唯一的小問題是,對於條件2,它只評估新添加的變量的符號,但我需要用於評估所有變量的符號。但我已經通過用'[c(strsplit(vars.test,「\\ +」)[[1]])]替換了'[var] –