如何使用stargazer报告coxph回归的exp(coefs)？

Question

假如我有一个my.model

My.model <- coxph(Surv(stop, event) ~ (rx + size + number) * strata(enum),
      cluster = id, bladder1)

我想创建一个模型报告表格，其中包含 exp(coefs) 而不是 coefs。

stargazer(my.model)

是否有类似于exponentiate = TRUE的参数，可以报告exp(coefs)而不是coefs？还是说我需要在传递给stargazer()之前转换模型结果？

- Macosso

---

我认为唯一的方法是先手动转换系数。 - Rsquare

Answer 1

为了获得指数化系数, 需要添加参数 apply.coef = exp, p.auto = FALSE, t.auto = FALSE。

My.model <- coxph(Surv(stop, event) ~ rx + size + number,
                  cluster = id, bladder)

原始模型未转换系数

stargazer(My.model, align=TRUE, 
          type="text",  digits = 3)
================================================
                         Dependent variable:    
                     ---------------------------
                                stop            
------------------------------------------------
rx                             -0.540*          
                               (0.200)          
                                                
size                           -0.055           
                               (0.070)          
                                                
number                        0.193***          
                               (0.046)          
                                                
------------------------------------------------
Observations                     340            
R2                              0.064           
Max. Possible R2                0.971           
Log Likelihood                -588.104          
Wald Test                12.510*** (df = 3)     
LR Test                  22.321*** (df = 3)     
Score (Logrank) Test     25.183*** (df = 3)     
================================================
Note: se in parenthesis *p<0.1; **p<0.05; ***p<0.01

使用参数apply.coef = exp进行指数运算。

stargazer(My.model, align=TRUE, apply.coef = exp,
          type="text",  digits = 3)

================================================
                         Dependent variable:    
                     ---------------------------
                                stop            
------------------------------------------------
rx                            0.583***          
                               (0.200)          
                                                
size                          0.947***          
                               (0.070)          
                                                
number                        1.213***          
                               (0.046)          
                                                
------------------------------------------------
Observations                     340            
R2                              0.064           
Max. Possible R2                0.971           
Log Likelihood                -588.104          
Wald Test                12.510*** (df = 3)     
LR Test                  22.321*** (df = 3)     
Score (Logrank) Test     25.183*** (df = 3)     
================================================
Note: se in parenthesis *p<0.1; **p<0.05; ***p<0.01

然而，正如您所看到的，星号提供了误导性推断，因为 t.stat = coef/se，但在这种情况下，指数化的系数被用作计算 t 统计量和 p 值的分子。

解决方案是添加参数 p.auto = FALSE 和 t.auto = FALSE，这将允许使用原始系数计算模型的 t 统计量和 p 值。

stargazer(My.model, align=TRUE, 
          type="text", apply.coef = exp, p.auto = FALSE, 
          t.auto = FALSE, digits = 3)


================================================
                         Dependent variable:    
                     ---------------------------
                                stop            
------------------------------------------------
rx                             0.583*           
                               (0.200)          
                                                
size                            0.947           
                               (0.070)          
                                                
number                        1.213***          
                               (0.046)          
                                                
------------------------------------------------
Observations                     340            
R2                              0.064           
Max. Possible R2                0.971           
Log Likelihood                -588.104          
Wald Test                12.510*** (df = 3)     
LR Test                  22.321*** (df = 3)     
Score (Logrank) Test     25.183*** (df = 3)     
================================================
Note: se in parenthesis *p<0.1; **p<0.05; ***p<0.01

此外，为了避免给读者带来困惑，您可以报告 t 统计量或 p 值，而不是标准误差。

stargazer(My.model, align=TRUE, 
          type="text", apply.coef = exp, p.auto = FALSE, 
          t.auto = FALSE, digits = 3, report=('vc*p'))

================================================
                         Dependent variable:    
                     ---------------------------
                                stop            
------------------------------------------------
rx                             0.583*           
                              p = 0.070         
                                                
size                            0.947           
                              p = 0.535         
                                                
number                        1.213***          
                              p = 0.005         
                                                
------------------------------------------------
Observations                     340            
R2                              0.064           
Max. Possible R2                0.971           
Log Likelihood                -588.104          
Wald Test                12.510*** (df = 3)     
LR Test                  22.321*** (df = 3)     
Score (Logrank) Test     25.183*** (df = 3)     
================================================
Note:                *p<0.1; **p<0.05; ***p<0.01

- Rsquare

Answer 2

您可以使用 `stargazer` 的 `apply.coef` 参数，像这样：

stargazer(model, apply.coef = exp)

根据以下帖子中的链接，您可能需要执行一些更复杂的操作来处理标准误差。

链接：在stargazer() LaTeX输出中使用比值而非logit

get.or.se <- function(model) {
  broom::tidy(model) %>%
    mutate(or = exp(estimate),
           var.diag = diag(vcov(model)),
           or.se = sqrt(or^2 * var.diag)) %>%
    select(or.se) %>% unlist %>% unname
}

代码来自链接，可能需要稍作调整以适应生存。

- dcsuka