Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者: 杨雨萱 (中山大学)
邮箱:yangyuxuanxuan1994@163.com
编者按:本文主要摘译自下文,特此致谢!
Source:Data Analysis Examples -Link-
目录
有些被解释变量只能取非负整数,例如专利个数、子女人数等。对于这一类计数数据的分析,经常会使用到泊松回归 (poisson regression)。接下来,我们以一个具体的例子来讲解。在该模型中,因变量 num_awards 为一所高中学生获得的奖励数量,自变量包括学生类型 prog (如职业、普通或学术) 以及他们的数学期末考试成绩 math。
. cnssc install lxhuse, lianxh replace
. lxhuse poisson_sim.dta, clear
. sum num_awards math
Variable | Obs Mean Std. dev. Min Max
-------------+--------------------------------------
num_awards | 200 .63 1.052921 0 6
math | 200 52.645 9.368448 33 75
. tabstat num_awards, by(prog) stats(mean sd n)
Summary for variables: num_awards
Group variable: prog (type of program)
prog | Mean SD N
---------+---------------------------
general | .2 .4045199 45
academic | 1 1.278521 105
vocation | .24 .5174506 50
---------+---------------------------
Total | .63 1.052921 200
-------------------------------------
. histogram num_awards, discrete freq scheme(s1mono)
. poisson num_awards i.prog math, vce(robust)
Iteration 0: log pseudolikelihood = -182.75759
Iteration 1: log pseudolikelihood = -182.75225
Iteration 2: log pseudolikelihood = -182.75225
Poisson regression Number of obs = 200
Wald chi2(3) = 80.15
Prob > chi2 = 0.0000
Log pseudolikelihood = -182.75225 Pseudo R2 = 0.2118
------------------------------------------------------------------------
| Robust
num_awards | Coefficient std. err. z P>|z| [95% conf. interval]
-----------+------------------------------------------------------------
prog |
academic | 1.084 0.322 3.37 0.001 0.453 1.715
vocation | 0.370 0.401 0.92 0.357 -0.417 1.157
|
math | 0.070 0.010 6.71 0.000 0.050 0.091
_cons | -5.247 0.648 -8.10 0.000 -6.516 -3.978
------------------------------------------------------------------------
math 的系数是 0.070,表明数学分数每增加一单位,得奖的平均次数是原来的
为了帮助评估模型拟合的效果,可以使用 estat gof
命令来获得拟合优度卡方检验。这不是对模型系数的检验,而是对模型形式的测试,即泊松模型形式是否适合我们的数据?
. estat gof
Deviance goodness-of-fit = 189.4496
Prob > chi2(196) = 0.6182
Pearson goodness-of-fit = 212.1437
Prob > chi2(196) = 0.2040
可以看出,模型拟合得相当好,因为拟合优度卡方检验在统计上不显著。如果该检验在统计上具有显著性,则表明该数据与模型不很吻合。在这种情况下,我们可以尝试确定是否有遗漏的自变量、线性假设是否成立、以及否存在过度分散的问题。
有时,我们可能希望将回归结果表示为事件发生比率,我们可以使用 irr
选项。这些 IRR 值与上面输出的系数之间存在如下换算关系:
. poisson, irr
Poisson regression Number of obs = 200
Wald chi2(3) = 80.15
Prob > chi2 = 0.0000
Log pseudolikelihood = -182.75225 Pseudo R2 = 0.2118
-------------------------------------------------------------------
| Robust
num_awards | IRR std. err. z P>|z| [95% conf. interval]
-----------+-------------------------------------------------------
prog |
academic | 2.956 0.951 3.37 0.001 1.573 5.555
vocation | 1.447 0.581 0.92 0.357 0.659 3.179
|
math | 1.073 0.011 6.71 0.000 1.051 1.095
_cons | 0.005 0.003 -8.10 0.000 0.001 0.019
-------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.
上述结果表明,保持其他变量不变,academic 的平均得奖率是参考组 general 的 2.956 倍。同样,保持其他变量不变, vocation 的平均得奖率是参考组 general 的 1.447 倍。数学成绩每增加一个单位,得奖事件发生率的百分比就增加 7.3%。
回顾我们的模型:
上述模型可以该写为:
可以看出,系数在 log(y) 尺度上具有 加法效应,IRR 在 y 尺度上具有乘法效应。为了更好地理解这个模型,下面使用 margin
命令来计算 prog 的每一层级的预测计数,在模型中保留所有其他变量的平均值。
. margins prog, atmeans
Adjusted predictions Number of obs = 200
Model VCE: Robust
Expression: Predicted number of events, predict()
At: 1.prog = .225 (mean)
2.prog = .525 (mean)
3.prog = .25 (mean)
math = 52.645 (mean)
---------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-----------+---------------------------------------------------------
prog |
general | 0.211 0.063 3.37 0.001 0.088 0.334
academic | 0.625 0.089 7.05 0.000 0.451 0.799
vocation | 0.306 0.083 3.69 0.000 0.144 0.468
---------------------------------------------------------------------
可以看出,在保持数学成绩在均值水平和 prog 取值为 1 水平上时,预测得到的平均得奖数大约是 0.211。同理,在保持数学成绩在均值水平和 prog 取值为 2 水平上时,预测得到的平均得奖数大约是为 0.625。在保持数学成绩在均值水平时和 prog 取值为 3 水平上时,预测得到的平均得奖数大约是 0.306。并且,prog 取值为 2 水平上的预测计数是 prog 取值为 1 水平上的 (.6249446/.211411)= 2.96 倍,这与我们在 IRR 输出表中看到的结果相一致。
下面我们将获得数学成绩从 35 到 75 (以 10 为增量) 的预测平均得奖数。
. margins, at(math=(35(10)75)) vsquish
Predictive margins Number of obs = 200
Model VCE: Robust
Expression: Predicted number of events, predict()
1._at: math = 35
2._at: math = 45
3._at: math = 55
4._at: math = 65
5._at: math = 75
--------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
-----+--------------------------------------------------------
_at |
1 | 0.131 0.036 3.66 0.000 0.061 0.201
2 | 0.264 0.048 5.57 0.000 0.171 0.358
3 | 0.533 0.058 9.27 0.000 0.421 0.646
4 | 1.076 0.122 8.82 0.000 0.837 1.315
5 | 2.170 0.412 5.27 0.000 1.363 2.976
--------------------------------------------------------------
上表结果显示,当 prog 变量观测值固定时,当数学成绩是 35 时,平均得奖数约为 0.131;当数学成绩是 75 分时,平均得奖数约为 2.170;其他同理。如果我们比较 math = 35 和 math = 45的平均得奖数时,计算得到的比率是 (0.2644714/0.1311326) = 2.017,即成绩每增加 10 单位,后者对应的得奖数是前者的 2.017 倍,与 IRR 值计算的结果 1.0727^10 = 2.017 一致。
如果想多获得一些其他指标,可以使用用户编写的 fitstat
命令,如果想比较模型,这些信息可能会有所帮助。
. fitstat
Measures of Fit for poisson of num_awards
Log-Lik Intercept Only: -231.864 Log-Lik Full Model: -182.752
D(195): 365.505 LR(3): 98.223
Prob > LR: 0.000
McFadden's R2: 0.212 McFadden's Adj R2: 0.190
Maximum Likelihood R2: 0.388 Cragg & Uhler's R2: 0.430
AIC: 1.878 AIC*n: 375.505
BIC: -667.667 BIC': -82.328
可以使用下面的命令来绘制事件的预测数量。
. twoway scatter c1 c2 c3 math, connect(l l l) sort ///
> ytitle("Predicted Count") ylabel( ,nogrid) legend(rows(3)) ///
> legend(ring(0) position(10)) scheme(s1mono)
=
该图表明,大多数奖项预计被学术计划 (prog=2) 中的学生拿走,尤其是在学生有更高数学成绩时。预计获奖人数最少的是普通计划课程的学生 (prog=1)。
负二项回归用于计数变量建模,通常用于过分散的计数结果变量。泊松回归的局限是泊松分布的期望与方差相等,被称为 “均等分散”,但这个特征与实际数据不符。
为此,Stata 命令提供了一个 LR 检验,其原假设为 “不存在过度分散,应使用泊松回归”。在完成负二项回归估计后,Stata 将自动输出检验结果 “LR test of alpha=0”,其中过度分散参数 “alpha=0” 对应于 “泊松回归”;即如果原假设 “alpha=0” 成立,则可使用泊松回归。如果使用选择项 “r” (稳健标准误),则 Stata 将只输出 alpha 的 95% 置信区间。根据这个置信区间,同样可以检验原假设 “alpha=0”。
下面我们以一个具体的例子,来看一下 Stata 具体操作流程,但由于和泊松回归过程类似,下面仅给出 Stata 命令。根据该数据构建的模型为,因变量为一两所高中学生的出勤天数 daysabs,自变量包括学生类型 prog (如职业、普通或学术) 以及他们的数学期末考试成绩 math。
lxhuse nb_data.dta, clear
summarize daysabs math
histogram daysabs, discrete freq scheme(s1mono)
nbreg daysabs math i.prog
test 2.prog 3.prog
nbreg, irr
margins prog, atmeans
margins, at(math=(0(20)100)) vsquish
fitstat
predict c
sort math
twoway (line c math if prog==1) ///
(line c math if prog==2) ///
(line c math if prog==3), ///
ytitle("Predicted Mean Value of Response") ylabel(0(3)15 ,nogrid) ///
xtitle("Math Score") legend(order(1 "prog=1" 2 "prog=2" 3 "prog=3")) scheme(s1mono)
如果计数数据中含有大量的 0 值,则可考虑使用 “零膨胀泊松回归” 或 “零膨胀负二项回归”。此外,理论表明,多余的零可以看作是由一个单独的过程产生的计数值,多余的零可以独立建模。因此,零膨胀模型有两部分,泊松计数模型和用于预测多余零的 logit 模型。
究竟应该使用标准的泊松回归还是零膨胀泊松回归?Stata 提供了一个 “Vuong” 统计量,其渐近分布为标准正态。如果 “Vuong” 统计量很大 (为正数),则应该选择零膨胀泊松回归 (或零膨胀负二项回归);反之,如果 “Vuong” 统计量很小 (为负数),则应该选择标准的泊松回归 (或标准的负二项回归)。
需要注意的是,不建议将零膨胀模型应用于小样本。另外,零膨胀负二项回归的 Stata 操作流程与零膨胀泊松回归类似,读者可自行研究。
断尾回归主要用于,部分因变量的观测值由于某些原因未被纳入分析中。例如,
. lxhuse truncreg.dta, clear
. summarize achiv langscore
Variable | Obs Mean Std. dev. Min Max
-----------+-------------------------------------
achiv | 178 54.23596 8.96323 41 76
langscore | 178 54.01124 8.944896 31 67
. tabstat achiv, by(prog) stats(n mean sd)
Summary for variables: achiv
Group variable: prog (type of program)
prog | N Mean SD
---------+------------------------------
general | 40 51.575 7.97074
academic | 101 56.89109 9.018759
vocation | 37 49.86486 7.276912
---------+------------------------------
Total | 178 54.23596 8.96323
----------------------------------------
. histogram achiv, bin(6) start(40) freq normal
下面我们使用 truncreg
命令来估计一个截断的回归模型。truncreg
命令中的 ll()
选项表示发生左截断的值,ul()
选项用于指示右截断值。
. truncreg achiv langscore i.prog, ll(40)
(0 obs truncated)
Fitting full model:
Iteration 0: log likelihood = -598.11669
Iteration 1: log likelihood = -591.68358
Iteration 2: log likelihood = -591.31208
Iteration 3: log likelihood = -591.30981
Iteration 4: log likelihood = -591.30981
Truncated regression
Limit: Lower = 40 Number of obs = 178
Upper = +inf Wald chi2(3) = 54.76
Log likelihood = -591.30981 Prob > chi2 = 0.0000
-------------------------------------------------------------------------
achiv | Coefficient Std. err. z P>|z| [95% conf. interval]
-----------+-------------------------------------------------------------
langscore | 0.713 0.114 6.22 0.000 0.488 0.937
|
prog |
academic | 4.065 2.055 1.98 0.048 0.038 8.093
vocation | -1.136 2.670 -0.43 0.671 -6.369 4.097
|
_cons | 11.302 6.773 1.67 0.095 -1.973 24.576
-----------+-------------------------------------------------------------
/sigma | 8.755 0.667 13.13 0.000 7.448 10.062
-------------------------------------------------------------------------
可以看出,langscore 系数为 0.713,表明语言成绩每增加一个单位,预期成绩就会增加 0.713 分。academic 系数为 4.065,表明在保持其他条件一定时,与 general 项目的同学相比,academic 项目的同学的预测成绩高了约 4.07 分。为了确定 prog 变量本身是否具有统计意义,我们可以使用 test
命令来获得该变量的二自由度检验。
. test 2.prog 3.prog
(1) [eq1]2.prog = 0
(2) [eq1]3.prog = 0
chi2( 2) = 7.19
Prob > chi2 = 0.0274
卡方检验表明,拒绝原假设,表明 prog 是具有统计学意义的预测因子。我们还可以使用 margins
命令来获得期望的边际效应。
. margins prog
Predictive margins Number of obs = 178
Model VCE: OIM
Expression: Linear prediction, predict()
-------------------------------------------------------------------
| Delta-method
| Margin std. err. z P>|z| [95% conf. interval]
----------+--------------------------------------------------------
prog |
general | 49.789 1.897 26.24 0.000 46.070 53.507
academic | 53.854 1.150 46.83 0.000 51.600 56.108
vocation | 48.653 2.140 22.73 0.000 44.458 52.848
--------------------------------------------------------------------
. marginsplot
在上表中的结果,我们可以看到,在 prog 等于 general 时,avchiv 的期望均值为 49.789。同理,在 prog 等于 academic 时,avchiv 的期望均值为 53.854。在 prog 等于 vocation 时,avchiv 的期望均值为 48.65。
truncreg
命令输出既不包含
. predict p
. correlate p achiv
. display r(rho)^2
.30519203
计算结果为 0.305,这个值是对回归中
tobit 模型又称归并回归 (censored regression) 模型,用于估计因变量中存在左归并或右归并时变量之间的线性关系(也称为上下归并)。具体来说,当某个值大于或等于某一阈值时,就会出现上述归并,因此真实值可能等于某一阈值,但也可能更高。那些高于或低于某个阈值的值将被归并。例如,在 20 世纪 80 年代,联邦法律限制速度计读数不超过每小时 85 英里。所以,如果你想通过马力和引擎大小来预测一辆车的最高时速,你得到的读数不会高于 85,即使这辆车的实际行驶速度更快。这是一个典型的数据右截 (右截) 的例子。我们唯一能确定的是这些车辆的时速至少为 85 英里。
我们以某个具体例子来说明,其中因变量为学业能力分数 apt,自变量为阅读成绩 read、数学成绩 math、以及学生类型 prog。
. lxhuse tobitreg.dta, clear
. summarize apt read math
Variable | Obs Mean Std. dev. Min Max
----------+-------------------------------------------
apt | 200 640.035 99.21903 352 800
read | 200 52.23 10.25294 28 76
math | 200 52.645 9.368448 33 75
. tabulate prog
type of |
program | Freq. Percent Cum.
------------+-----------------------------------
academic | 45 22.50 22.50
general | 105 52.50 75.00
vocational | 50 25.00 100.00
------------+-----------------------------------
Total | 200 100.00
. histogram apt, normal bin(10) xline(800)
从上面的直方图中,我们可以看到数据中的归并 (censored)。需要注意的是,重点考察 y 在哪里被截断的,要有明确的来源和根据。接下来我们使用 tobit
命令进行回归分析,在 tobit
命令中的 ul()
选项表示右归并开始的值 (即上限),ll()
选项用于指示左归并值 (下限),在本例中不需要这个值。具体模型解释和 Truncated Regression 类似,下面仅提供代码,不做详细解释。
lxhuse tobitreg.dta, clear
tobit apt read math i.prog, ul(800)
test 2.prog 3.prog
test 2.prog = 3.prog
predict yhat
correlate apt yhat
fitstat
Note:产生如下推文列表的 Stata 命令为:
lianxh logit probit
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh