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

[celestialgod (攸藍)]

# 我創的例子
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)
→ celestialgod:我想第二種方法應該是你比較喜歡的
推 diplazium:推c大的認真用心!