R 中 β1=β2 的 F 检验
F test for β1=β2 in R
如果我的模型看起来像这样,Y=β0+β1X1+β2X2+β3X3+β4X4,我想在 R 中对 β1=β2 执行 F 检验 (5%),我该怎么做?
我在网上能找到的唯一教程是处理 β1=β2=0 的,但这不是我在这里要找的。
这是 R 中的一个示例,用于测试 vs
的系数是否与 am
的系数相同:
data(mtcars)
mod <- lm(mpg ~ hp + disp + vs + am, data=mtcars)
library(car)
linearHypothesis(mod, "vs=am")
# Linear hypothesis test
#
# Hypothesis:
# vs - am = 0
#
# Model 1: restricted model
# Model 2: mpg ~ hp + disp + vs + am
#
# Res.Df RSS Df Sum of Sq F Pr(>F)
# 1 28 227.07
# 2 27 213.52 1 13.547 1.7131 0.2016
multcomp
包中的 glht
函数可以做到这一点(以及其他)。例如,如果您的模型是
mod1 <-lm( y ~ x1 + x2 + x3 + x4)
那么你可以使用:
summary(multcomp::glht(mod1, "x1-x2=0"))
运行 有约束和无约束的模型,然后使用 anova 来比较它们。没有使用包。
mod1 <- lm(mpg ~ cyl + disp + hp + drat, mtcars)
mod2 <- lm(mpg ~ I(cyl + disp) + hp + drat, mtcars) # constraint imposed
anova(mod2, mod1)
给予:
Analysis of Variance Table
Model 1: mpg ~ I(cyl + disp) + hp + drat
Model 2: mpg ~ cyl + disp + hp + drat
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 252.95
2 27 244.90 1 8.0513 0.8876 0.3545
基础计算如下。它给出了与上面相同的结果。
L <- matrix(c(0, 1, -1, 0, 0), 1) # hypothesis is L %*% beta == 0
q <- nrow(L) # 1
co <- coef(mod1)
resdf <- df.residual(mod1) # = nobs(mod1) - length(co) = 32 - 5 = 27
SSH <- t(L %*% co) %*% solve(L %*% vcov(mod1) %*% t(L)) %*% L %*% co
SSH/q # F value
## [,1]
## [1,] 0.8876363
pf(SSH/q, q, resdf, lower.tail = FALSE) # p value
## [,1]
## [1,] 0.3544728
如果我的模型看起来像这样,Y=β0+β1X1+β2X2+β3X3+β4X4,我想在 R 中对 β1=β2 执行 F 检验 (5%),我该怎么做? 我在网上能找到的唯一教程是处理 β1=β2=0 的,但这不是我在这里要找的。
这是 R 中的一个示例,用于测试 vs
的系数是否与 am
的系数相同:
data(mtcars)
mod <- lm(mpg ~ hp + disp + vs + am, data=mtcars)
library(car)
linearHypothesis(mod, "vs=am")
# Linear hypothesis test
#
# Hypothesis:
# vs - am = 0
#
# Model 1: restricted model
# Model 2: mpg ~ hp + disp + vs + am
#
# Res.Df RSS Df Sum of Sq F Pr(>F)
# 1 28 227.07
# 2 27 213.52 1 13.547 1.7131 0.2016
multcomp
包中的 glht
函数可以做到这一点(以及其他)。例如,如果您的模型是
mod1 <-lm( y ~ x1 + x2 + x3 + x4)
那么你可以使用:
summary(multcomp::glht(mod1, "x1-x2=0"))
运行 有约束和无约束的模型,然后使用 anova 来比较它们。没有使用包。
mod1 <- lm(mpg ~ cyl + disp + hp + drat, mtcars)
mod2 <- lm(mpg ~ I(cyl + disp) + hp + drat, mtcars) # constraint imposed
anova(mod2, mod1)
给予:
Analysis of Variance Table
Model 1: mpg ~ I(cyl + disp) + hp + drat
Model 2: mpg ~ cyl + disp + hp + drat
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 252.95
2 27 244.90 1 8.0513 0.8876 0.3545
基础计算如下。它给出了与上面相同的结果。
L <- matrix(c(0, 1, -1, 0, 0), 1) # hypothesis is L %*% beta == 0
q <- nrow(L) # 1
co <- coef(mod1)
resdf <- df.residual(mod1) # = nobs(mod1) - length(co) = 32 - 5 = 27
SSH <- t(L %*% co) %*% solve(L %*% vcov(mod1) %*% t(L)) %*% L %*% co
SSH/q # F value
## [,1]
## [1,] 0.8876363
pf(SSH/q, q, resdf, lower.tail = FALSE) # p value
## [,1]
## [1,] 0.3544728