从 clogit 模型中检索摘要时出错

Question

我正在通过以下方式在多个数据集（具有数千个事件）上使用 clogit 训练条件逻辑回归模型：

library(survival) 
library(mgcv)

# load dataset 
df <- read.csv('1.csv')

model <- clogit(case ~ 
           var1 +
           # pspline(var2, df = 3) +
           strata(var3),
           data = df)

print(model)
summary(model)

列类型为：case: int, var1: factor, var2: int, var3: int。

如果我保留此行注释：pspline(var2, df = 3) +，当数据集在每个层中都有足够的案例时，摘要打印工作正常，否则我会收到以下警告和非常大的标准错误：

Warning message in fitter(X, Y, strats, offset, init, control, weights = weights, :
“Loglik converged before variable  1,2,3 ; beta may be infinite. ”

但是，如果我使用行 pspline(var2, df = 3) +，那么即使数据集在每个层中没有足够的案例，我也不会收到此类警告。 print(model) 行有效，但是当我尝试访问模型的 summary 时出现以下错误：

Error in pchisq(cox$score, df, lower.tail = FALSE): Non-numeric argument to mathematical function
Traceback:

1. print(summary(model))
2. summary(model)
3. summary.coxph(model)
4. pchisq(cox$score, df, lower.tail = FALSE)

我需要访问摘要，因为我正在将系数打印到 csv 文件以供以后处理：summary(model)$coefficients 因为我正在多个文件上训练模型。

我找不到此行为的原因，如有任何帮助，我们将不胜感激。

编辑：06.26 最小可重现示例：

num_cases = 100
var3 = rep((1:num_cases), each=3)
case = rep(c(0, 1, 1), num_cases)
var1 = factor(sample(c("Low", "Medium", "High"), num_cases, replace=TRUE, prob = c(0.5,0.35,0.25)))
var2 = runif(num_cases * 3, 10, 35)

generated_data <- data.frame(var3, case, var1, var2)

model <- clogit(case ~ 
           var1 +
           pspline(var2, df = 3) +
           strata(var3),
           data = generated_data)

print(model)
summary(model)$coefficients

结果：

在 case ~ var1 后添加逗号不会产生错误。代码现在打印系数，但它 returns 的系数与我删除逗号并使用 print(model).
时返回的系数不同
上述代码在num_cases = 200时无法收敛。
在 case ~ var1 之后添加逗号会产生另一个警告：

Warning message in clogit(case ~ var1, +pspline(var2, df = 3) + strata(var3), data = generated_data): “weights ignored: not possible for the exact method”

Answer 1

发生的事情是 class 'clogit' 的对象通过继承传递给 summary.cph，这显然不是为它设计的。您可以从 print.clogit 函数中获取系数，这是在您请求模型结果时隐式调用的函数：

 model
Call:
clogit(case ~ var1 + pspline(var2, df = 3) + strata(var3), data = generated_data)

               coef exp(coef) se(coef)      z     p
var1Low    -0.02063   0.97958  0.23019 -0.090 0.929
var1Medium -0.02171   0.97852  0.23831 -0.091 0.927
ps(var2)3  -0.05886   0.94284  0.40900 -0.144 0.886
ps(var2)4  -0.10752   0.89806  0.60868 -0.177 0.860
ps(var2)5  -0.16100   0.85129  0.65381 -0.246 0.805
ps(var2)6  -0.23156   0.79330  0.63652 -0.364 0.716
ps(var2)7  -0.26708   0.76561  0.61080 -0.437 0.662
ps(var2)8  -0.23270   0.79239  0.59903 -0.388 0.698
ps(var2)9  -0.23075   0.79394  0.59781 -0.386 0.700
ps(var2)10 -0.27852   0.75690  0.60117 -0.463 0.643
ps(var2)11 -0.26878   0.76431  0.64828 -0.415 0.678
ps(var2)12 -0.24330   0.78404  0.84486 -0.288 0.773

Likelihood ratio test=0.92  on 5.04 df, p=0.9702
n= 300, number of events= 200

作为奖励，您将获得 LLR 测试值和 p 值。如果您只想要通常由汇总函数返回的排序矩阵，则对 print.coxph:

中的代码部分进行明显的修改

{ coef <- model$coefficients
se <- sqrt(diag(model$var))
if (is.null(coef) | is.null(se)) 
    stop("Input is not valid")
if (is.null(model$naive.var)) {
    tmp <- cbind(coef, exp(coef), se, coef/se, pchisq((coef/se)^2, 
                                                      1, lower.tail = FALSE))
    dimnames(tmp) <- list(names(coef), c("coef", "exp(coef)", 
                                         "se(coef)", "z", "p"))} }

然后 tmp 将是所需的矩阵：

 tmp
                  coef exp(coef)  se(coef)           z         p
var1Low    -0.02062820 0.9795831 0.2301878 -0.08961464 0.9285935
var1Medium -0.02171186 0.9785221 0.2383051 -0.09110952 0.9274056
ps(var2)3  -0.05886243 0.9428365 0.4089999 -0.14391797 0.8855652
ps(var2)4  -0.10752217 0.8980566 0.6086828 -0.17664728 0.8597855
ps(var2)5  -0.16099704 0.8512946 0.6538079 -0.24624519 0.8054924
ps(var2)6  -0.23155828 0.7932965 0.6365191 -0.36378845 0.7160160
ps(var2)7  -0.26708193 0.7656103 0.6108000 -0.43726573 0.6619186
ps(var2)8  -0.23269795 0.7923929 0.5990301 -0.38845785 0.6976772
ps(var2)9  -0.23074825 0.7939393 0.5978122 -0.38598783 0.6995057
ps(var2)10 -0.27852015 0.7569030 0.6011671 -0.46329903 0.6431500
ps(var2)11 -0.26877900 0.7643121 0.6482824 -0.41460174 0.6784335
ps(var2)12 -0.24329853 0.7840374 0.8448572 -0.28797593 0.7733652

我不确定错误报告是否有必要，但如果您不这么认为，那么 Thomas Lumley 就是 clogit 的作者。 clogit 的帮助页面没有描述摘要方法，并且 print.clogit 和分派到 print.coxph 方法似乎被用于通常分配给 summary 的目的。

此外，系数本身可以通过 model$coef 获得，但不会返回完整的系数矩阵和辅助统计估计，

从 clogit 模型中检索摘要时出错

Error retrieving the summary from clogit model

regression

r

spline

survival-analysis

logistic-regression