我正在使用R
:svolik
和中的两个数据集est
。就背景而言,我开发了一种新的概念衡量标准(立法权力共享),并用它来复制之前的一项研究:Svolik (2012)。练习的目的是看看使用我的测量结果是否不同。
这是svolik
数据:https://drive.google.com/file/d/1nCBhRXNcBrLEr6-R2pkyuQ9mCtJKkdmm/view ?usp=sharing
这是est
数据:https://drive.google.com/file/d/1D-UmHSi9LIEsmY5VBvU8nxu8u1gix7Ay/view ?usp=sharing
我从 Svolik 用于生成结果的数据集开始。我成功地重现了他的结果(图中的模型1、3、5)。然后,我将他的数据集与包含我的新度量的数据集合并,丢弃任何不完全匹配的观察结果:
# load original data (the data used to produce original results)
svolik <- read_dta("svolik.dta")
# load data containing my new measure
est <- read.csv("Merging with Svolik.csv")
# merge
final <- merge(svolik, est, by = c("ccode", "year"), all = FALSE)
Run Code Online (Sandbox Code Playgroud)
接下来,我再次运行他的模型,但将他的立法机构变量替换为我的立法权力共享变量(图中的模型 2、4 和 6)。请注意,尽管数据涵盖同一时间段,但原始模型和我自己的模型包含的观测值数量略有不同(2,903 个而不是 2,934 个)。
我一生都无法弄清楚为什么我会得到这些额外的观察结果。我的猜测是它与合并/重复或类似的事情有关。您觉得这可能是个问题吗?如果是这样,您知道有什么方法可以找出这些观察结果吗?解决方案可能很简单,我可能只是想太多了。任何意见,将不胜感激!请注意,我尝试使用不同的合并策略 --- left_join
in dplyr()
--- 但这不起作用。
请注意,我正在 Stata 中运行结果。以下是原始结果(即模型 1、3 和 5)的 Stata 代码:
* SURVIVAL ANALYSIS
use "leaders, institutions, covariates, updated tvc.dta"
* NATURAL DEATHS
gen c_natural=censoring
replace c_natural=0 if exit!="natural"
replace c_natural=. if exit==""
tab c_natural
stset t, id(leadid) failure(c_natural)
stcox legislature lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age
outreg2 using survival, replace ctitle(natural, leg) tex nonotes bdec(3) e(all) ef
* COUPS
gen c_coup= censoring
replace c_coup=0 if exit!="coup"
replace c_coup=. if exit==""
stset t, id(leadid) failure(c_coup)
* REMOVE SOM DUPLICATE OBSERVATIONS
* drop if (t[_n-1]==t & leadid[_n-1]== leadid)
stset t, id(leadid) failure(c_coup)
stcox legislature lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age
outreg2 using survival, ctitle(coups, leg) tex nonotes bdec(3) e(all) ef
* REVOLTS
gen c_revolt= censoring
replace c_revolt=0 if exit!="revolt"
replace c_revolt=. if exit==""
tab c_revolt
stset t, id(leadid) failure(c_revolt)
* COMMUNIST LEFT OUT BECAUSE IT IS A PERFECT PREDICTOR
stcox legislature lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED mil cw age
outreg2 using survival, ctitle(revolt, leg) tex nonotes bdec(3) e(all) ef
Run Code Online (Sandbox Code Playgroud)
以下是新结果的 Stata 代码(即模型 2、4 和 6):
* SURVIVAL ANALYSIS
use "merged_test.dta"
* NATURAL DEATHS
gen c_natural=censoring
replace c_natural=0 if exit!="natural"
replace c_natural=. if exit==""
tab c_natural
stset t, id(leadid) failure(c_natural)
stcox estimate lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age
outreg2 using survival, replace ctitle(natural, leg) tex nonotes bdec(3) e(all) ef
* COUPS
gen c_coup= censoring
replace c_coup=0 if exit!="coup"
replace c_coup=. if exit==""
stset t, id(leadid) failure(c_coup)
* REMOVE SOM DUPLICATE OBSERVATIONS
* drop if (t[_n-1]==t & leadid[_n-1]== leadid)
stset t, id(leadid) failure(c_coup)
stcox estimate lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age
outreg2 using survival, ctitle(coups, leg) tex nonotes bdec(3) e(all) ef
* REVOLTS
gen c_revolt= censoring
replace c_revolt=0 if exit!="revolt"
replace c_revolt=. if exit==""
tab c_revolt
stset t, id(leadid) failure(c_revolt)
* COMMUNIST LEFT OUT BECAUSE IT IS A PERFECT PREDICTOR
stcox estimate lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED mil cw age
outreg2 using survival, ctitle(revolt, leg) tex nonotes bdec(3) e(all) ef
Run Code Online (Sandbox Code Playgroud)
您的问题不是svolik
有 2903 个观察值,final
而是有 2934 个观察值,因此是由合并中的一些重复行引起final
的超集。svolik
您永远不会有 2903 行可以在两个回归中使用。这是因为并非 Svolik 回归中使用的所有完整观测值在数据集中都有相应的行est
。首先我们来了解一下 2903 个观测值是svolik
从哪里来的:
svolik_reg_cols <- c("legislative", "lgdp_1", "growth_1", "exportersoffuelsmainlyoil_EL2008", "ethfrac_FIXED", "communist", "mil", "cw", "age")\nsvolik_is_complete <- complete.cases(svolik[, svolik_reg_cols])\nsum(svolik_is_complete) # 2903\n
Run Code Online (Sandbox Code Playgroud)\n正如您所看到的,它是回归中所有列的完整案例数。现在让我们final
使用您的 join 方法对 执行相同的操作:
final <- merge(svolik, est, by = c("ccode", "year"), all = FALSE)\nfinal_reg_cols <- svolik_reg_cols\nfinal_reg_cols[final_reg_cols == "legislative"] <- "estimate"\nfinal_is_complete <- complete.cases(final[, final_reg_cols])\nsum(final_is_complete) # 2934\n
Run Code Online (Sandbox Code Playgroud)\n同样,2934 是任何协变量都没有缺失数据的观测值数量。
\n不过,让我们看看您要加入的数据集。中ccode
和中未出现的有 278 组。year
svolik
est
# How many ccode and year are in svolik but not est\ndplyr::anti_join(\n svolik,\n est,\n by = c("ccode", "year")\n) |>\n group_by(ccode, cabb, year) |>\n summarise(n = n()) |>\n arrange(desc(n)) |>\n print(n = 2)\n\n# # A tibble: 278 \xc3\x97 3\n# # Groups: ccode [39]\n# ccode year n\n# <dbl> <dbl> <int>\n# 1 990 1982 4\n# 2 947 2001 3\n# # \xe2\x80\xa6 with 276 more rows\n
Run Code Online (Sandbox Code Playgroud)\n这意味着根据您拥有的数据,不可能比较所有观察结果。
\n您有三个选择:
\n您将知道 1. 或 2. 是否可能。然而,由于您的分析目的似乎是将您的新指标与 Svolik 进行比较,因此 3. 似乎是一种合理的方法,特别是当您最终不会删除很多行时。
\n首先找到公共行(有2830行)并保存到dta
:
all_complete <- complete.cases(final[, c("estimate", svolik_reg_cols)])\nsum(all_complete) # 2830\nfinal_complete <- final[all_complete, ]\nwrite_dta(final_complete, "./tmp/svolik_est_merged.dta")\n
Run Code Online (Sandbox Code Playgroud)\n您现在可以在 Stata 中运行回归。首先像之前一样加载并准备数据:
\nuse svolik_est_merged.dta, clear\n\n* NATURAL DEATHS\ncap drop c_natural c_coup c_revolt _d _t _t0\ngen c_natural=censoring\nreplace c_natural=0 if exit!="natural"\nreplace c_natural=. if exit==""\ntab c_natural\n\nstset t, id(leadid) failure(c_natural)\n
Run Code Online (Sandbox Code Playgroud)\n现在运行 Svolik 回归。您可以看到有 2830 个观测值:
\nstcox legislative lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age\n
Run Code Online (Sandbox Code Playgroud)\nCox regression with Breslow method for ties\n\nNo. of subjects = 383 Number of obs = 2,830\nNo. of failures = 40\nTime at risk = 3,098\n LR chi2(9) = 28.46\nLog likelihood = -157.48569 Prob > chi2 = 0.0008\n\n-------------------------------------------------------------------------------------\n _t | Haz. ratio Std. err. z P>|z| [95% conf. interval]\n--------------------+----------------------------------------------------------------\n legislative | 1.006541 .0088251 0.74 0.457 .9893923 1.023988\n lgdp_1 | 1.437144 .3138694 1.66 0.097 .9366983 2.204962\n growth_1 | 1.010814 .0283629 0.38 0.701 .956725 1.067962\nexportersoffue~2008 | 2.487166 1.205382 1.88 0.060 .9620061 6.430308\n ethfrac_FIXED | 1.011694 .00645 1.82 0.068 .9991306 1.024415\n communist | 2.0526 1.610128 0.92 0.359 .4411573 9.550262\n mil | 1.06844 .3944057 0.18 0.858 .5182463 2.202744\n cw | 4.15784 2.325053 2.55 0.011 1.389562 12.44106\n age | 1.057077 .0172812 3.40 0.001 1.023744 1.091496\n-------------------------------------------------------------------------------------\n
Run Code Online (Sandbox Code Playgroud)\n然后运行回归:
\nstcox estimate lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age\n
Run Code Online (Sandbox Code Playgroud)\n输出:
\nCox regression with Breslow method for ties\n\nNo. of subjects = 383 Number of obs = 2,830\nNo. of failures = 40\nTime at risk = 3,098\n LR chi2(9) = 28.00\nLog likelihood = -157.71273 Prob > chi2 = 0.0010\n\n-------------------------------------------------------------------------------------\n _t | Haz. ratio Std. err. z P>|z| [95% conf. interval]\n--------------------+----------------------------------------------------------------\n estimate | .9742007 .1278445 -0.20 0.842 .7532603 1.259946\n lgdp_1 | 1.506868 .3265272 1.89 0.058 .9854309 2.304222\n growth_1 | 1.007996 .028074 0.29 0.775 .954447 1.06455\nexportersoffue~2008 | 2.147553 1.257702 1.31 0.192 .6814636 6.767761\n ethfrac_FIXED | 1.011719 .0070275 1.68 0.093 .9980384 1.025587\n communist | 2.064115 1.619767 0.92 0.356 .4433766 9.609369\n mil | 1.018648 .3747256 0.05 0.960 .4953321 2.094845\n cw | 3.961413 2.202203 2.48 0.013 1.332464 11.77727\n age | 1.054575 .0174756 3.21 0.001 1.020873 1.089389\n-------------------------------------------------------------------------------------\n
Run Code Online (Sandbox Code Playgroud)\n又是 2830 个观察值。结果似乎与我非常相似:相同的两个协变量 (cw
和age
) 具有较小的 p 值,并且所有系数都接近 Svolik。如果您正在尝试开发一个衡量标准,该指标会告诉您一些新信息,但可能不是您想听到的。然而,如果您试图通过与既定指标进行比较来确定您的指标是否稳健,也许这是更好的消息。
当然,需要注意的是,要考虑无法连接的行是否完全随机丢失。如果存在系统性问题导致与因变量相关的数据丢失(例如,专制国家更有可能阻止在某些年份发布数据),则推断此处的类似结果可以推断为丢失可能是不合理的数据。查看此问题的一种方法是将包含 2830 行的 Svolik 结果与包含 2903 行的结果进行比较,看看它们是否存在显着差异。
\n在评论中回答您的问题。您想查看c_natural
其中包含哪些失败,svolik
但没有final
:
use svolik_est_merged.dta, clear\n\n* NATURAL DEATHS\ncap drop c_natural c_coup c_revolt _d _t _t0\ngen c_natural=censoring\nreplace c_natural=0 if exit!="natural"\nreplace c_natural=. if exit==""\ntab c_natural\n\nstset t, id(leadid) failure(c_natural)\n
Run Code Online (Sandbox Code Playgroud)\n您可以对c_coup
或执行相同的操作c_revolt
。