本文给出了一个简单线性回归模型t统计量的计算公式。
t=β1/SE(Beta1)
SE(beta1)=sqrt((RSS/var(x1))*(1/n-2))
如果我想用R来做一个简单的例子,我无法得到与R中的线性模型相同的结果。
x <- c(1,2,4,8,16)
y <- c(1,2,3,4,5)
mod <- lm(y~x)
summary(mod)
Call:
lm(formula = y ~ x)
Residuals:
1 2 3 4 5
-0.74194 0.01613 0.53226 0.56452 -0.37097
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.50000 0.44400 3.378 0.0431 *
x 0.24194 0.05376 4.500 0.0205 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6558 on 3 degrees of freedom
Multiple R-squared: 0.871, Adjusted R-squared: 0.828
F-statistic: 20.25 on 1 and 3 DF, p-value: 0.02049
如果我用手做这件事,我会得到另一种价值。
var(x)
37.2
sum(resid(mod)^2)
1.290323
beta1=0.24194
SE(Beta1)=sqrt(1.290323/37.2)*(1/3) SE(beta1)=0.1075269
所以t= 0.24194/0.1075269=2.250042
那么为什么我的计算精确到R值的一半呢?这和一到两个尾测试有关系吗?t(0.05/2)为3.18。
你好,简
发布于 2017-10-30 13:43:30
不同的结果是由se(beta)
公式中缺少的一个项引起的。它应该是:
se(beta) = sqrt((1 / (n - 2)) * rss / (var(x) * (n - 1)))
公式通常写成:
se(beta) = sqrt((1 / (n - 2)) * rss / sum((x - mean(x)) ^ 2))
而不是var(x)
。
为了完整起见,还需要进行计算检查:
reprex::reprex_info()
#> Created by the reprex package v0.1.1.9000 on 2017-10-30
x <- c(1, 2, 4, 8, 16)
y <- c(1, 2, 3, 4, 5)
n <- length(x)
mod <- lm(y ~ x)
summary(mod)
#>
#> Call:
#> lm(formula = y ~ x)
#>
#> Residuals:
#> 1 2 3 4 5
#> -0.74194 0.01613 0.53226 0.56452 -0.37097
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 1.50000 0.44400 3.378 0.0431 *
#> x 0.24194 0.05376 4.500 0.0205 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.6558 on 3 degrees of freedom
#> Multiple R-squared: 0.871, Adjusted R-squared: 0.828
#> F-statistic: 20.25 on 1 and 3 DF, p-value: 0.02049
mod_se_b <- summary(mod)$coefficients[2, 2]
rss <- sum(resid(mod) ^ 2)
se_b <- sqrt((1 / (n - 2)) * rss / (var(x) * (n - 1)))
all.equal(se_b, mod_se_b)
#> [1] TRUE
https://stackoverflow.com/questions/46991667
复制相似问题