JMa*_*ino 14 r cox-regression survival-analysis
我在使用coxph()时遇到了一些麻烦.我有两个分类变量:"tecnologia"和"pais",我想评估"pais"对"tecnologia"的可能的交互作用."tecnologia"是一个具有2个级别的变量因子:gps和convencional.而"pais"为2级:PT和ES.我不知道为什么这个警告不断出现.这是代码和输出:
cox_AC<-coxph(Surv(dados_temp$dias_seg,dados_temp$status)~tecnologia*pais,data=dados_temp)
Warning message:
In coxph(Surv(dados_temp$dias_seg, dados_temp$status) ~ tecnologia * :
X matrix deemed to be singular; variable 3
> cox_AC
Call:
coxph(formula = Surv(dados_temp$dias_seg, dados_temp$status) ~
tecnologia * pais, data = dados_temp)
coef exp(coef) se(coef) z p
tecnologiagps -0.152 0.859 0.400 -0.38 7e-01
paisPT 1.469 4.345 0.406 3.62 3e-04
tecnologiagps:paisPT NA NA 0.000 NA NA
Likelihood ratio test=23.8 on 2 df, p=6.82e-06 n= 127, number of events= 64
Run Code Online (Sandbox Code Playgroud)
我正在打开关于这个主题的另一个问题,虽然几个月前我做了一个类似的问题,因为我再次面临同样的问题,还有其他数据.而这次我确定这不是一个与数据相关的问题.
有人能帮助我吗?谢谢
更新: 问题似乎不是一个完美的分类
> xtabs(~status+tecnologia,data=dados)
tecnologia
status conv doppler gps
0 39 6 24
1 30 3 34
> xtabs(~status+pais,data=dados)
pais
status ES PT
0 71 8
1 49 28
> xtabs(~tecnologia+pais,data=dados)
pais
tecnologia ES PT
conv 69 0
doppler 1 8
gps 30 28
Run Code Online (Sandbox Code Playgroud)
dar*_*sco 13
这是一个简单的例子,似乎可以重现你的问题:
> library(survival)
> (df1 <- data.frame(t1=seq(1:6),
s1=rep(c(0, 1), 3),
te1=c(rep(0, 3), rep(1, 3)),
pa1=c(0,0,1,0,0,0)
))
t1 s1 te1 pa1
1 1 0 0 0
2 2 1 0 0
3 3 0 0 1
4 4 1 1 0
5 5 0 1 0
6 6 1 1 0
> (coxph(Surv(t1, s1) ~ te1*pa1, data=df1))
Call:
coxph(formula = Surv(t1, s1) ~ te1 * pa1, data = df1)
coef exp(coef) se(coef) z p
te1 -23 9.84e-11 58208 -0.000396 1
pa1 -23 9.84e-11 100819 -0.000229 1
te1:pa1 NA NA 0 NA NA
Run Code Online (Sandbox Code Playgroud)
现在让我们像这样寻找'完美分类':
> (xtabs( ~ s1+te1, data=df1))
te1
s1 0 1
0 2 1
1 1 2
> (xtabs( ~ s1+pa1, data=df1))
pa1
s1 0 1
0 2 1
1 3 0
Run Code Online (Sandbox Code Playgroud)
请注意,1
用于pa1
精确预测状态s1
等于的值0
.也就是说,根据您的数据,如果您知道,pa1==1
那么您可以肯定s1==0
.因此,在此设置中拟合Cox模型是不合适的,并且会导致数值误差.这可以看出来
> coxph(Surv(t1, s1) ~ pa1, data=df1)
Run Code Online (Sandbox Code Playgroud)
给
Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights, :
Loglik converged before variable 1 ; beta may be infinite.
Run Code Online (Sandbox Code Playgroud)
在拟合模型之前查看这些交叉表非常重要.在考虑涉及交互的那些模型之前,还值得从更简单的模型开始.
如果我们df1
手动添加交互项如下:
> (df1 <- within(df1,
+ te1pa1 <- te1*pa1))
t1 s1 te1 pa1 te1pa1
1 1 0 0 0 0
2 2 1 0 0 0
3 3 0 0 1 0
4 4 1 1 0 0
5 5 0 1 0 0
6 6 1 1 0 0
Run Code Online (Sandbox Code Playgroud)
然后检查一下
> (xtabs( ~ s1+te1pa1, data=df1))
te1pa1
s1 0
0 3
1 3
Run Code Online (Sandbox Code Playgroud)
我们可以看到它是一个无用的分类器,即它无助于预测状态s1
.
当组合所有3个项时,钳工确实设法产生数值te1
,pe1
即使pe1
是如上所述的完美预测器.然而,看看系数的值及其误差表明它们是不可信的.
编辑 @JMarcelino:如果您查看coxph
示例中第一个模型的警告消息,您将看到警告消息:
2: In coxph(Surv(t1, s1) ~ te1 * pa1, data = df1) :
X matrix deemed to be singular; variable 3
Run Code Online (Sandbox Code Playgroud)
这可能是你得到的同样的错误,并且是由于这个分类问题.此外,你的第三个交叉表xtabs(~ tecnologia+pais, data=dados)
不作为的表作为重要status
的interaction term
.您可以先手动添加交互项,如上例所示,然后检查交叉表.或者你可以说:
> with(df1,
table(s1, pa1te1=pa1*te1))
pa1te1
s1 0
0 3
1 3
Run Code Online (Sandbox Code Playgroud)
也就是说,我注意到你的第三个表中的一个单元格有一个零(conv
,PT
)意味着你没有观察到这种预测变量的组合.这会在试图适应时引起问题.
一般而言,结果应该对所有预测因子水平都有一些值,预测因子不应将结果分类为全部或全部或50/50.
编辑2 @ user75782131是的,一般来说xtabs
或类似的交叉表应该在结果和预测因子是离散的模型中执行,即具有有限的编号.水平.如果存在"完美分类",则预测模型/回归可能不合适.例如,逻辑回归(结果是二元)以及Cox模型都是如此.