R-код для проверки разности коэффициентов регрессоров от одной регрессии
Я хочу проверить, отличаются ли коэффициенты в одной линейной регрессии друг от друга, или же хотя бы один из них значительно отличается от одного определенного значения, скажем, 0, это кажется довольно интуитивно понятным в Stata. Например
webuse iris
reg iris seplen sepwid petlen
seplen==sepwid==petlen
seplen==sepwid==petlen==0
Интересно, как я могу это сделать, если я хочу проверить это в R?
2 ответа
car
Пакет имеет простую функцию, чтобы сделать это.
Во-первых, подойдет ваша модель:
model <- lm(Sepal.Length ~ Sepal.Width + Petal.Length, data = iris)
Чем вы можете проверить различные линейные гипотезы, используя linearHypothesis
функция, например:
library(car)
# tests if the coefficient of Sepal.Width = Petal.Length
linearHypothesis(model, "Sepal.Width = Petal.Length")
Linear hypothesis test
Hypothesis:
Sepal.Width - Petal.Length = 0
Model 1: restricted model
Model 2: Sepal.Length ~ Sepal.Width + Petal.Length
Res.Df RSS Df Sum of Sq F Pr(>F)
1 148 16.744
2 147 16.329 1 0.4157 3.7423 0.05497 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Other examples:
# tests if the coefficient of Sepal.Width = 0
linearHypothesis(model, "Sepal.Width = 0")
# tests if the coefficient of Sepal.Width = 0.5
linearHypothesis(model, "Sepal.Width = 0.5")
# tests if both are equal to zero
linearHypothesis(model, c("Sepal.Width = 0", "Petal.Length = 0"))
Вы можете сравнить список коэффициентов от каждой соответствующей модели (скажем, mod1 и mod2), как в:
diff=merge(mod1$coefficients, mod2$coefficients, by=0, all=TRUE)
diff[is.na(diff)]=0
diff$error=abs(diff$x-diff$y)
diff[order(diff$error, decreasing=TRUE),]
Это создает кадр данных, отсортированный по абсолютной величине разницы в коэффициентах, то есть:
Row.names x y error
1 (Intercept) -0.264189182 -0.060450853 2.037383e-01
6 id 0.003402056 0.000000000 3.402056e-03
3 b -0.001804978 -0.003357193 1.552215e-03
2 a -0.049900767 -0.049417150 4.836163e-04
4 c 0.013749907 0.013819799 6.989203e-05
5 d -0.004097366 -0.004110830 1.346320e-05
Если уклоны не те, что вам нужны, вы можете получить доступ к другим коэффициентам, используя функцию coef():
coef(summary(model))
Например, чтобы получить Pr(>|z|), используйте:
coef(summary(model))[,"Pr(>|z|)"]