作者celestialgod (攸藍)
看板R_Language
標題Re: [問題] 詢問有關迴歸分析 R的指令
時間Sat Aug 17 00:12:30 2013
# 我創的例子
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)
1F：→ celestialgod:我想第二種方法應該是你比較喜歡的 08/17 00:14
2F：推 diplazium:推c大的認真用心! 08/17 11:15