R stargazer 包输出：缺少 felm 回归的 F 统计（lfe 包）
R stargazer package output: Missing F statistic for felm regression (lfe package)
我正在尝试使用 stargazer 包来输出我的回归结果。我使用 lfe 包中的 felm 执行回归。 stargazer 输出表正确显示所有内容，但 F 统计值仍为空白。 lm 结果不会出现此问题。
原因是什么以及如何让 felm 回归的 F 统计值出现在 stargazer 输出中？
我知道我可以手动添加一行来显示 F 值，但如果可能的话我更喜欢更自动的方法。
下面是使用提供的数据的示例代码here
library(foreign)
temp_dat <- read.dta("http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.dta")
temp_lm <- lm(y ~ x, temp_dat)
temp_felm <- felm(y ~ x, temp_dat)
library(stargazer)
stargazer(temp_lm, temp_felm, type = "text")
输出：
====================================================================
Dependent variable:
------------------------------------
y
OLS felm
(1) (2)
--------------------------------------------------------------------
x 1.035*** 1.035***
(0.029) (0.029)
Constant 0.030 0.030
(0.028) (0.028)
--------------------------------------------------------------------
Observations 5,000 5,000
R2 0.208 0.208
Adjusted R2 0.208 0.208
Residual Std. Error (df = 4998) 2.005 2.005
F Statistic 1,310.740*** (df = 1; 4998)
====================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
似乎没有在 stargazer 中实现自动化的方法，这是一个很好的包，但不可扩展。选项 keep.stat = "f" 不会为 felm 对象生成 f-stat。但是，texreg 有一个选项，其中包括 felm 个对象的 f-stat。
library(foreign);library(texreg);library(lfe)
temp_dat <- read.dta("http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.dta")
temp_lm <- lm(y ~ x, temp_dat)
temp_felm <- felm(y ~ x, temp_dat)
screenreg(list(temp_lm, temp_felm), include.fstatistic = T)
产生：
==================================================
Model 1 Model 2
--------------------------------------------------
(Intercept) 0.03 0.03
(0.03) (0.03)
x 1.03 *** 1.03 ***
(0.03) (0.03)
--------------------------------------------------
R^2 0.21
Adj. R^2 0.21
Num. obs. 5000 5000
F statistic 1310.74
RMSE 2.01
R^2 (full model) 0.21
R^2 (proj model) 0.21
Adj. R^2 (full model) 0.21
Adj. R^2 (proj model) 0.21
F statistic (full model) 1310.74
F (full model): p-value 0.00
F statistic (proj model) 1310.74
F (proj model): p-value 0.00
==================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
createTexreg 功能允许您选择要提取和显示的特定统计数据。您首先需要编写一个小函数来从 summary.felm 对象中提取对象，然后将其转换为 texreg 对象。
extract.felm <- function(model, include.f.full = TRUE,
include.f.proj = TRUE,
include.rsquared = TRUE,
include.adjrs = TRUE,
include.nobs = TRUE, ...) {
s <- summary(model, ...)
names <- rownames(s$coefficients)
co <- s$coefficients[, 1]
se <- s$coefficients[, 2]
pval <- s$coefficients[, 4]
gof <- numeric()
gof.names <- character()
gof.decimal <- logical()
if (include.rsquared == TRUE) {
rs <- s$r.squared
gof <- c(gof, rs)
gof.names <- c(gof.names, "R$^2$")
gof.decimal <- c(gof.decimal, TRUE)
}
if (include.adjrs == TRUE) {
adj <- s$adj.r.squared
gof <- c(gof, adj)
gof.names <- c(gof.names, "Adj.\ R$^2$")
gof.decimal <- c(gof.decimal, TRUE)
}
if (include.nobs == TRUE) {
n <- s$N
gof <- c(gof, n)
gof.names <- c(gof.names, "Num.\ obs.")
gof.decimal <- c(gof.decimal, FALSE)
}
if (include.f.full == TRUE) {
ffs <- s$fstat
ffpval <- round(s$F.fstat[4],4)
gof <- c(gof, ffs, ffpval)
gof.names <- c(gof.names, "F statistic (Full model)", "F (full model): p-value")
gof.decimal <- c(gof.decimal, TRUE, TRUE)
}
if (include.f.proj == TRUE) {
fps <- s$P.fstat[5]
fppval <- s$P.fstat[4]
gof <- c(gof, fps, fppval)
gof.names <- c(gof.names, "F statistic (proj. model)", "F (proj. model): p-value") #Modify the names as you see fit
gof.decimal <- c(gof.decimal, TRUE, TRUE)
}
tr <- createTexreg(
coef.names = names,
coef = co,
se = se,
pvalues = pval,
gof.names = gof.names,
gof = gof,
gof.decimal = gof.decimal
)
return(tr)
}
setMethod("extract", signature = className("felm", "stats"),
definition = extract.felm)
现在，运行 函数，设置参数 include.f.prof = F 并将其发送到 screenreg:
> m <- extract.felm(temp_felm, include.f.proj = F)
> screenreg(list(temp_lm, m))
==================================================
Model 1 Model 2
--------------------------------------------------
(Intercept) 0.03 0.03
(0.03) (0.03)
x 1.03 *** 1.03 ***
(0.03) (0.03)
--------------------------------------------------
R^2 0.21 0.21
Adj. R^2 0.21 0.21
Num. obs. 5000 5000
RMSE 2.01
F statistic (Full model) 1310.74
F (full model): p-value 0.00
==================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
对于您的问题，这是一个不完全令人满意的解决方案。
希望对你有所帮助
从lfe:::summary.felm获取F统计量：
(Fstat <- summary(temp_felm)$F.fstat)
###############
F df1 df2 p
1.310740e+03 1.000000e+00 4.998000e+03 4.252163e-255
然后在使用stargazer生成的table的注释中添加F统计：
stars <- ifelse(Fstat[4]<0.01,"***",ifelse(Fstat[4]<0.05,"**",ifelse(Fstat[4]<0.1,"*","")))
notes <- c(paste("F Statistic(full) \hfill ", round(Fstat[1],3),stars," (df = ",Fstat[2],"; ", Fstat[3],")",sep=""),
"Note:\hfill *p<0.1; **p<0.05; ***p<0.01")
stargazer(temp_felm, notes=notes, notes.append=F, notes.label="", notes.align="l")
最终的 LaTeX table 是：
我正在尝试使用 stargazer 包来输出我的回归结果。我使用 lfe 包中的 felm 执行回归。 stargazer 输出表正确显示所有内容，但 F 统计值仍为空白。 lm 结果不会出现此问题。
原因是什么以及如何让 felm 回归的 F 统计值出现在 stargazer 输出中？
我知道我可以手动添加一行来显示 F 值，但如果可能的话我更喜欢更自动的方法。
下面是使用提供的数据的示例代码here
library(foreign)
temp_dat <- read.dta("http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.dta")
temp_lm <- lm(y ~ x, temp_dat)
temp_felm <- felm(y ~ x, temp_dat)
library(stargazer)
stargazer(temp_lm, temp_felm, type = "text")
输出：
====================================================================
Dependent variable:
------------------------------------
y
OLS felm
(1) (2)
--------------------------------------------------------------------
x 1.035*** 1.035***
(0.029) (0.029)
Constant 0.030 0.030
(0.028) (0.028)
--------------------------------------------------------------------
Observations 5,000 5,000
R2 0.208 0.208
Adjusted R2 0.208 0.208
Residual Std. Error (df = 4998) 2.005 2.005
F Statistic 1,310.740*** (df = 1; 4998)
====================================================================
Note: *p<0.1; **p<0.05; ***p<0.01
似乎没有在 stargazer 中实现自动化的方法，这是一个很好的包，但不可扩展。选项 keep.stat = "f" 不会为 felm 对象生成 f-stat。但是，texreg 有一个选项，其中包括 felm 个对象的 f-stat。
library(foreign);library(texreg);library(lfe)
temp_dat <- read.dta("http://www.kellogg.northwestern.edu/faculty/petersen/htm/papers/se/test_data.dta")
temp_lm <- lm(y ~ x, temp_dat)
temp_felm <- felm(y ~ x, temp_dat)
screenreg(list(temp_lm, temp_felm), include.fstatistic = T)
产生：
==================================================
Model 1 Model 2
--------------------------------------------------
(Intercept) 0.03 0.03
(0.03) (0.03)
x 1.03 *** 1.03 ***
(0.03) (0.03)
--------------------------------------------------
R^2 0.21
Adj. R^2 0.21
Num. obs. 5000 5000
F statistic 1310.74
RMSE 2.01
R^2 (full model) 0.21
R^2 (proj model) 0.21
Adj. R^2 (full model) 0.21
Adj. R^2 (proj model) 0.21
F statistic (full model) 1310.74
F (full model): p-value 0.00
F statistic (proj model) 1310.74
F (proj model): p-value 0.00
==================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
createTexreg 功能允许您选择要提取和显示的特定统计数据。您首先需要编写一个小函数来从 summary.felm 对象中提取对象，然后将其转换为 texreg 对象。
extract.felm <- function(model, include.f.full = TRUE,
include.f.proj = TRUE,
include.rsquared = TRUE,
include.adjrs = TRUE,
include.nobs = TRUE, ...) {
s <- summary(model, ...)
names <- rownames(s$coefficients)
co <- s$coefficients[, 1]
se <- s$coefficients[, 2]
pval <- s$coefficients[, 4]
gof <- numeric()
gof.names <- character()
gof.decimal <- logical()
if (include.rsquared == TRUE) {
rs <- s$r.squared
gof <- c(gof, rs)
gof.names <- c(gof.names, "R$^2$")
gof.decimal <- c(gof.decimal, TRUE)
}
if (include.adjrs == TRUE) {
adj <- s$adj.r.squared
gof <- c(gof, adj)
gof.names <- c(gof.names, "Adj.\ R$^2$")
gof.decimal <- c(gof.decimal, TRUE)
}
if (include.nobs == TRUE) {
n <- s$N
gof <- c(gof, n)
gof.names <- c(gof.names, "Num.\ obs.")
gof.decimal <- c(gof.decimal, FALSE)
}
if (include.f.full == TRUE) {
ffs <- s$fstat
ffpval <- round(s$F.fstat[4],4)
gof <- c(gof, ffs, ffpval)
gof.names <- c(gof.names, "F statistic (Full model)", "F (full model): p-value")
gof.decimal <- c(gof.decimal, TRUE, TRUE)
}
if (include.f.proj == TRUE) {
fps <- s$P.fstat[5]
fppval <- s$P.fstat[4]
gof <- c(gof, fps, fppval)
gof.names <- c(gof.names, "F statistic (proj. model)", "F (proj. model): p-value") #Modify the names as you see fit
gof.decimal <- c(gof.decimal, TRUE, TRUE)
}
tr <- createTexreg(
coef.names = names,
coef = co,
se = se,
pvalues = pval,
gof.names = gof.names,
gof = gof,
gof.decimal = gof.decimal
)
return(tr)
}
setMethod("extract", signature = className("felm", "stats"),
definition = extract.felm)
现在，运行 函数，设置参数 include.f.prof = F 并将其发送到 screenreg:
> m <- extract.felm(temp_felm, include.f.proj = F)
> screenreg(list(temp_lm, m))
==================================================
Model 1 Model 2
--------------------------------------------------
(Intercept) 0.03 0.03
(0.03) (0.03)
x 1.03 *** 1.03 ***
(0.03) (0.03)
--------------------------------------------------
R^2 0.21 0.21
Adj. R^2 0.21 0.21
Num. obs. 5000 5000
RMSE 2.01
F statistic (Full model) 1310.74
F (full model): p-value 0.00
==================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
对于您的问题，这是一个不完全令人满意的解决方案。
希望对你有所帮助
从lfe:::summary.felm获取F统计量：
(Fstat <- summary(temp_felm)$F.fstat)
###############
F df1 df2 p
1.310740e+03 1.000000e+00 4.998000e+03 4.252163e-255
然后在使用stargazer生成的table的注释中添加F统计：
stars <- ifelse(Fstat[4]<0.01,"***",ifelse(Fstat[4]<0.05,"**",ifelse(Fstat[4]<0.1,"*","")))
notes <- c(paste("F Statistic(full) \hfill ", round(Fstat[1],3),stars," (df = ",Fstat[2],"; ", Fstat[3],")",sep=""),
"Note:\hfill *p<0.1; **p<0.05; ***p<0.01")
stargazer(temp_felm, notes=notes, notes.append=F, notes.label="", notes.align="l")
最终的 LaTeX table 是：