r - 对来自不同线性模型的系数进行假设检验

Question

想象一下我有两个单独的lm对象

data(mtcars)

lm1 <- lm(mpg ~ wt, data = mtcars)
lm2 <- lm(mpg ~ wt + disp, data = mtcars)

在这种情况下，我想比较两个wt系数，并对两个模型中的系数相等的空值进行假设检验（出于技术原因，我实际上需要有两个模型，而不仅仅是包含交互）

Answer 1

由于您想对估计值进行假设检验，我建议使用完全贝叶斯模型，该模型将为您提供每个变量的完整后验分布。

rstanarm基于Stan, 并提供模仿通常lm,glm语法的便捷功能；如果您想了解更多关于Stan/的信息RStan，请参见此处。

基于每个变量的后验分布，我们可以执行例如t检验和 Kolmogorov-Smirnov 检验来比较每个变量的完整后验密度。

这是您可以执行的操作：

1. 执行模型拟合。

   library(rstanarm);
   lm1 <- stan_lm(mpg ~ wt, data = mtcars, prior = NULL);
   lm2 <- stan_lm(mpg ~ wt + disp, data = mtcars, prior = NULL);
   

   请注意使用 运行完全贝叶斯线性模型是多么容易rstanarm。

2. 提取所有共享系数的后验密度（在本例中为(Intercept)和wt）。

   library(tidyverse);
   shared.coef <- intersect(names(coef(lm1)), names(coef(lm2)));
   shared.coef;
   #[1] "(Intercept)" "wt"
   df1 <- lm1 %>%
       as.data.frame() %>%
       select(one_of(shared.coef)) %>%
       mutate(model = "lm1");
   df2 <- lm2 %>%
       as.data.frame() %>%
       select(one_of(shared.coef)) %>%
       mutate(model = "lm2");
   

   4000 次 MCMC 绘制的后验密度存储在两个data.frames 中。

3. 我们绘制后密度。

   # Plot posterior densities for all common parameters
   bind_rows(df1, df2) %>%
       gather(var, value, 1:length(shared.coef)) %>%
       ggplot(aes(value, colour = model)) +
           geom_density() +
           facet_wrap(~ var, scale = "free");
   

   

4. 我们比较了t检验和 KS 检验中每个共享参数的后验密度分布。在这里，我使用库broom来整理输出。

   # Perform t test and KS test
   library(broom);
   res <- lapply(1:length(shared.coef), function(i)
       list(t.test(df1[, i], df2[, i]), ks.test(df1[, i], df2[, i])));
   names(res) <- shared.coef;
   lapply(res, function(x) bind_rows(sapply(x, tidy)));
   #$`(Intercept)`
   #   estimate estimate1 estimate2 statistic p.value parameter  conf.low conf.high
   #1 -4.497093  30.07725  34.57434 -104.8882       0  7155.965 -4.581141 -4.413045
   #2        NA        NA        NA    0.7725       0        NA        NA        NA
   #                              method alternative
   #1            Welch Two Sample t-test   two.sided
   #2 Two-sample Kolmogorov-Smirnov test   two-sided
   #
   #$wt
   #   estimate estimate1 estimate2 statistic      p.value parameter  conf.low
   #1 0.1825202 -3.097777 -3.280297  9.120137 1.074479e-19  4876.248 0.1432859
   #2        NA        NA        NA  0.290750 0.000000e+00        NA        NA
   #  conf.high                             method alternative
   #1 0.2217544            Welch Two Sample t-test   two.sided
   #2        NA Two-sample Kolmogorov-Smirnov test   two-sided
   #
   #There were 12 warnings (use warnings() to see them)
   

   （警告来自绑定行时不相等的因子水平，可以忽略。）