作者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