最近有个小伙伴问到了这样的一个问题，他想对所有自变量的组合进行回归，然后从中筛选系数估计值均显著的。
为此我们使用 auto 数据集进行演示：
假如我们想考虑这些变量：headroom trunk length displacement 对 mpg 的影响：
reg mpg headroom trunk length displacement
|
自变量是 headroom trunk length displacement，可能的组合有下面这些：
reg mpg displacement
reg mpg length
reg mpg trunk
reg mpg headroom
reg mpg length displacement
reg mpg trunk displacement
reg mpg trunk length
reg mpg headroom displacement
reg mpg headroom length
reg mpg headroom trunk
reg mpg trunk length displacement
reg mpg headroom length displacement
reg mpg headroom trunk displacement
reg mpg headroom trunk length
reg mpg headroom trunk length displacement
|
这段代码是使用下面的代码生成的：
sysuse auto, clear
tuples headroom trunk length displacement
forvalues i = 1/`ntuples' {
di "reg mpg `tuple`i''"
}
|
逐行运行上面的代码，然后观测哪些是自变量系数估计值都显著的。对于组合较少的情况可以这样做，不过我们能不能自动把这些模型的结果保存起来，然后直接进行筛选呢？
提取保存模型结果的方法很多，这里我打算使用 esttab 输出。首先循环回归保存回归结果：
local m = ""
forvalues i = 1/`ntuples' {
qui reg mpg `tuple`i''
estadd local cmdline = e(cmdline), replace
est store m`i'
local m = "`m' m`i'"
}
esttab `m' using temp.csv, replace stats(cmdline) p(%6.4f)
|
这样就把所有模型的结果保存到了 temp.csv 中，把这个文件读入 Stata 中：
import delimited using temp.csv, clear
|
然后去除所有的 =””：
foreach i of varlist _all {
replace `i' = ustrregexs(1) if ustrregexm(`i', `"="(.*)""')
}
|
保留 3～13 行，然后把 p 值所在行的 v1 替换成 pvalue + 自变量的名字：
keep in 3/13
replace v1 = "pvalue_" + v1[_n - 1] if v1 == ""
|
由于我们只关注显著性，所以保留 pvalue_* 行即可：
*- 只关注显著性
keep if index(v1, "pvalue_") | v1 == "cmdline"
|
然后我们把数据转置下：
*- 安装 sxpose：ssc install sxpose
sxpose, clear
|
把第一行作为变量名：
*- 安装 nrow：ssc install nrow
nrow 1
*- 把 cmdline 放到第一列
|
然后去掉 p 值的括号：
foreach i of varlist pvalue_* {
replace `i' = ustrregexs(1) if ustrregexm(`i', "\((.*)\)")
}
destring, replace
|
然后我们就可以筛选自变量系数均显著的模型了：
egen ind = rowmax(pvalue_*)
replace ind = (ind < 0.05)
drop pvalue__cons
egen ind2 = rowmax(pvalue_*)
replace ind2 = (ind < 0.05)
|
这样我们就解决了这个问题。
不过如果组合特别多的时候（超过 300 个）使用这个方法会报错，因为 Stata 最多只能在内存中存储 300 个模型。为此我们再来看一下组合特别多的时候该怎么办？
首先我们生产多一些变量：
sysuse auto, clear
*- 生成一些随机变量
foreach i of varlist _all {
cap gen `i'2 = `i' * runiform()
cap gen `i'3 = `i' * runiform()
}
tuples headroom trunk length displacement headroom2 trunk2 length2 displacement2 headroom3
*- 一共 511 个模型
di "`ntuples'"
*> 511
|
由于 Stata 没办法直接保存 511 个模型，所以我们还是循环把每个模型直接保存成 csv 文件：
sysuse auto, clear
foreach i of varlist _all {
cap gen `i'2 = `i' * runiform()
cap gen `i'3 = `i' * runiform()
}
tuples headroom trunk length displacement headroom2 trunk2 length2 displacement2 headroom3
di "`ntuples'"
cap mkdir "res"
forvalues i = 1/`ntuples' {
di "`i'"
qui {
tempfile myfile
qui regress mpg `tuple`i''
estadd local cmdline = e(cmdline), replace
est store m`i'
esttab m`i' using "`myfile'.csv", replace stats(cmdline) p(%6.4f)
est clear
preserve
import delimited using "`myfile'.csv", clear
gen model = "m`i'"
save "res/`i'", replace
restore
}
}
|
循环合并所有的数据：
use res/1, clear
forval i = 2/511 {
append using res/`i'
}
|
简单处理下：
replace v1 = ustrregexs(1) if ustrregexm(v1, `"="(.*)""')
replace v2 = ustrregexs(1) if ustrregexm(v2, `"="(.*)""')
replace v2 = ustrregexs(1) if ustrregexm(v2, "\((.*)\)")
drop if v2 == "1"
drop if v1 == "p-values in parentheses"
drop if v1 == "* p<0.05, ** p<0.01, *** p<0.001"
drop if v2 == "mpg"
replace v1 = "pvalue_" + v1[_n - 1] if v1 == ""
|
只关注显著性：
keep if index(v1, "pvalue_") | v1 == "cmdline"
|
转换成宽数据：
*- 安装 spread：ssc install spread
spread v1 v2
|
如果不关心常数项可以删除 pvalue__cons：
drop pvalue__cons
destring, replace
|
然后就可以筛选自变量系数都显著的模型了：
egen ind = rowmax(pvalue_*)
replace ind = (ind < 0.05)
keep if ind == 1
|
这样也可以解决这个问题。
点击这里跳转到 RStata 短书平台获取附件：如何使用 Stata 生成所有自变量的组合进行回归并筛选自变量系数估计值均显著的模型？
评论