Re: [問題] 詢問有關迴歸分析 R的指令

看板R_Language作者celestialgod (攸藍)時間11年前 (2013/08/17 00:12)推噓1(1推 0噓 1→)

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

# 我創的例子 n = 50; p = 5 X = matrix(rnorm(n*p),n) data = data.frame(y = cbind(1,X) %*% 1:6 + rnorm(n), X) lm.fit = lm(y~., data) # 改自anova.lm這個function anova_f2 = function(lm.fit){ if(class(lm.fit)!="lm") stop("It just work for lm.") if(!is.null(lm.fit$weights)) stop("It does not work for weighted lm.") SSE = sum(resid(lm.fit)^2) SST = sum((lm.fit$model$y - mean(lm.fit$model$y))^2) ss = c(SST-SSE, SSE, SST) p = lm.fit$rank; n = length(resid(lm.fit)) df = c(p-1, n-p, n-1) ms = ss / df f = ms[1] / ms[2] p = pf(f, df[1], df[2], lower.tail = FALSE) table = data.frame(df, ss, ms, f, p) table[2:3, 4:5] = NA dimnames(table) <- list(c("model", "Residuals", "Total"), c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)")) table } anova_f2(lm.fit) # result Df Sum Sq Mean Sq F value Pr(>F) model 5 4548.74366 909.748731 870.5031 8.457429e-43 Residuals 44 45.98369 1.045084 NA NA Total 49 4594.72734 93.769946 NA NA # another way n = 50; p = 5 X = matrix(rnorm(n*p),n) y = cbind(1,X) %*% 1:6 + rnorm(n) anova(lm(y~X)) # covariate 用矩陣就可以顯示model的ANOVA table了 # result Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(>F) X 5 5840.0 1168.01 1268.5 < 2.2e-16 *** Residuals 44 40.5 0.92 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ※ 引述《sseeaann (屌哥)》之銘言： : 我知道ANOVA這個指令可以叫出變異數分析表 : 不過這指令跑出來的是一個MODEL的各變數的變異數分析表 : 而小弟我要的則是一個MODEL的SST,SSR,SSE的變異數分析表來分析整個MODEL : 知道指令的大大麻煩請解答一下吧 -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 218.164.72.75 ※ 編輯: celestialgod 來自: 218.164.72.75 (08/17 00:13)