2

我正在尝试使用该modelsummary包制作统计模型汇总表。

我进行了 6 次回归。我对其中 2 个进行了联合假设检验，其中包括某些变量（特别是educ和exper）作为​​解释变量。我想将这些测试的 F 统计信息包含在表中，但无法从测试中提取信息并将其导入表中。

这里有更多细节。

(1) 我做了什么，第 1 部分：包、数据和分析

library(tidyverse)
library(wooldridge)   # Data
library(estimatr)     # Estimation
library(car)          # F test
library(modelsummary) # Table
library(kableExtra)   # Table

wage2 <- wage2 %>% 
  mutate(lwage = log(wage))

model1 <- lm_robust(lwage ~ educ, 
                    data = wage2, se_type = "HC1")
model2 <- lm_robust(lwage ~ exper, 
                    data = wage2, se_type = "HC1")
model3 <- lm_robust(lwage ~ educ + exper, 
                    data = wage2, se_type = "HC1")
F3 <- linearHypothesis(model3, test = "F", 
                       c("educ = 0", "exper = 0"))
model4 <- lm_robust(lwage ~ educ + tenure, 
                    data = wage2, se_type = "HC1")
model5 <- lm_robust(lwage ~ exper + tenure, 
                    data = wage2, se_type = "HC1")
model6 <- lm_robust(lwage ~ educ + exper + tenure, 
                    data = wage2, se_type = "HC1")
F6 <- linearHypothesis(model6, test = "F", 
                       c("educ = 0", "exper = 0"))

联合假设检验的结果如下。

F3

## Linear hypothesis test
## 
## Hypothesis:
## educ = 0
## exper = 0
## 
## Model 1: restricted model
## Model 2: lwage ~ educ + exper
## 
##   Res.Df Df      F    Pr(>F)    
## 1    934                        
## 2    932  2 67.715 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F6

## Linear hypothesis test
## 
## Hypothesis:
## educ = 0
## exper = 0
## 
## Model 1: restricted model
## Model 2: lwage ~ educ + exper + tenure
## 
##   Res.Df Df      F    Pr(>F)    
## 1    933                        
## 2    931  2 63.592 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(2) 我做了什么，第 2 部分：表格

regs <- list()
regs[['Model 1']] <- model1
regs[['Model 2']] <- model2
regs[['Model 3']] <- model3
regs[['Model 4']] <- model4
regs[['Model 5']] <- model5
regs[['Model 6']] <- model6

gm <- tribble(
  ~raw, ~clean, ~fmt,
  "adj.r.squared", "$\\bar{R}^2$", 3,
  "nobs", "Sample size", 0
  )

rows <- tribble(
  ~term, ~'Model 1', ~'Model 2', ~'Model 3', 
  ~'Model 4', ~'Model 5', ~'Model 6', 
  '$F$ statstics', '', '', '', '', '', '',
  '$H_0: \\beta_{\\rm{educ}} = 0, \\beta_{\\rm{exper}} = 0$', 
  '', '', '67.715***', '', '', '63.592***',
  '', '', '', '(0.000)', '', '', '(0.000)'
)
attr(rows, 'position') <- c(9:11)

table <- msummary(regs, 
         estimate = "{estimate}{stars}",
         gof_map = gm, 
         notes = list("*** p < 0.01, ** p < 0.05, * p < 0.1. 
                      Heteroskedasticity-consistent 
                      standard errors in brackets. 
                      P values in brackets for F statistics."),
         add_rows = rows) %>% 
  row_spec(c(1, 3, 5, 7, 9, 10), 
           extra_css = "border-bottom: hidden") 
table

上面的代码给了我一张漂亮的桌子。（不幸的是，我似乎无法发布它，因为我没有获得足够的声誉。）

(3) 我想知道的

在上面，我手动输入了 F 统计量、星数和 p 值rows来获得我想要的表，但我想知道是否有办法从表中提取信息F3并将F6其导入到表中。

4

2 回答 2

2

如果我理解正确，您想要：

  1. 在表格底部添加两行：F 统计量和关联的 p 值
  2. 这些行应仅添加到您手动选择的某些模型中

实现这一点的最简单方法是使用add_rows您在示例中使用的参数。你的方法是正确的！

如果你想要一个更通用和自动化的策略，你可以定义一个glance_custom.lm_robust方法。

你的例子很复杂，所以我冒昧地简化了它，把它变成了一个最小的可重现的例子。您会注意到我创建并检查了一个新属性以确定哪些模型将在表中进行线性假设检验。

（注意：回答你的问题让我发现了一个小错误，它混合了表格底部的行顺序。这个错误现在已经修复。下面粘贴的表格是使用modelsummary0.9.2.9000 的开发版本生成的，你可以通过调用安装remotes::install_github("vincentarelbundock/modelsummary") ：）

library(modelsummary)
library(estimatr)
library(tibble)
library(car)

mod1 <- lm_robust(mpg ~ am + vs, mtcars)
mod2 <- lm_robust(mpg ~ am + vs + hp, mtcars)
mod3 <- lm_robust(mpg ~ am + vs + wt, mtcars)
mod4 <- lm_robust(mpg ~ am + vs + qsec, mtcars)

# models 2 and 4 will have a linear hypothesis test
attr(mod2, "FTEST") <- TRUE
attr(mod4, "FTEST") <- TRUE

glance_custom.lm_robust <- function(x) {
  # check if we should include the linear hypothesis test in this specific model
  if (!isTRUE(attr(x, "FTEST"))) return(NULL)

  # conduct linear hypothesis test
  ftest <- linearHypothesis(x, test = "F", c("am = 0", "vs = 0"))

  # return a 1-row tibble with each statistic in separate columns
  out <- tibble(
    "$H_0: \\beta_{\\rm{vs}} = 0, \\beta_{\\rm{am}} = 0$" = ftest[["F"]][2],
    "     " = sprintf("(%.3f)", ftest[["Pr(>F)"]][2]))
  return(out)
}

modelsummary(list(mod1, mod2, mod3, mod4), escape = FALSE)

在此处输入图像描述

于 2021-10-15T11:40:07.287 回答
0

如果我理解正确，您想从模型中导出数据。直接喊出来吧：

F3$F

和

F3$`Pr(>F)`
于 2021-10-15T09:07:15.953 回答