我试图在R中重现最初在Stata中获得的考克斯生存模型结果。下面是Stata代码:
stset t, id(leadid) failure(c_coup)
stcox legislature lgdp_1 growth_1 exportersoffuelsmainlyoil_EL2008 ethfrac_FIXED communist mil cw age
下面是我在R
中编写的代码:
# Load survival package
library(survival)
# Set the survival object
surv_obj <- Surv(data$t, data$c_coup)
# Run model
m1 <- coxph(surv_obj ~ legislature + lgdp_1 + growth_1 + exportersoffuelsmainlyoil_EL2008 + ethfrac_FIXED + communist + mil + cw + age, data = data, method = "breslow")
# Examine hazard ratios
exp(coef(m1))
但无论如何,我无法得到同样的结果。例如,Stata结果(我想产生的原始结果)中对立法机构的估计是0.298。在R中,它是0.1688371。任何建议将不胜感激。
以下是用于复制的数据集预览:
structure(list(t = structure(c(1, 2, 3, 4, 5, 6), label = "Current time in office", format.stata = "%9.0g"),
c_coup = structure(c(0, 0, 0, 0, 0, 0), format.stata = "%9.0g"),
leadid = structure(c("A2.2-208", "A2.2-208", "A2.2-208",
"A2.2-208", "A2.2-208", "A2.2-208"), label = "Leader ID", format.stata = "%13s"),
legislature = structure(c(1, 1, 1, 1, 1, 1), format.stata = "%9.0g"),
lgdp_1 = structure(c(7.68524360656738, 7.69938945770264,
7.54960918426514, 7.57916784286499, 7.6033992767334, 7.67089462280273
), format.stata = "%9.0g"), growth_1 = structure(c(6.35386085510254,
1.42463231086731, -13.910285949707, 3, 2.45273375511169,
6.98254346847534), label = "annual growth, t-1, Maddison", format.stata = "%9.0g"),
exportersoffuelsmainlyoil_EL2008 = structure(c(0, 0, 0, 0,
0, 0), format.stata = "%8.0g"), ethfrac_FIXED = structure(c(NA_real_,
NA_real_, NA_real_, NA_real_, NA_real_, NA_real_), label = "eth. frac", format.stata = "%8.0g"),
communist = structure(c(0, 0, 0, 0, 0, 0), label = "Communist Leader", format.stata = "%8.0g"),
mil = structure(c(1, 1, 1, 1, 1, 1), format.stata = "%9.0g"),
cw = structure(c(1, 1, 1, 1, 1, 1), format.stata = "%9.0g"),
age = structure(c(52, 53, 54, 55, 56, 57), label = "Current age", format.stata = "%9.0g")), row.names = c(NA,
-6L), class = c("tbl_df", "tbl", "data.frame"))
1条答案
按热度按时间j0pj023g1#
您可以尝试为每一行设置一个
t0
:输出: