2014-01-07 52 views
13

我在使用coxph()時遇到了一些麻煩。我有兩個分類變量:「tecnologia」和「pais」,我想評估「pais」對「tecnologia」的可能互動效果。「tecnologia」是一個變量因子,具有2個級別:gps和convencional。而「pais」分爲2個等級:PT和ES。我不知道爲什麼這個警告不斷出現。 下面的代碼和輸出:coxph()X矩陣被認爲是單數;

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 

我打開另一個問題,關於這個問題,儘管我在幾個月前做了一個類似的一個,因爲我再次面臨同樣的問題,與其他數據。這一次,我確定這不是一個數據相關的問題。

有人可以幫我嗎? 謝謝

UPDATE: 似乎該問題不會是一個完美的分類

> 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 
+3

這看起來像'完美的分類'(即當查看交互時,至少有一個因素組合*所有*觀察結果都有一定的狀態)。你是否通過變量及其交互來查看'status'的交叉表? – dardisco

+0

你是什麼意思?我不明白我在尋找什麼.. – JMarcelino

+1

1.結果是:'xtabs(〜tecnologica + pais,data = dados)'? 2.爲什麼不「輸入」你的數據,讓人們檢查出來而不是推測? –

回答

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 

現在讓我們來看看對於'完美分類'就像這樣:

> (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 
那的 1pa1 正是值預測爲 s1等於 0狀態

注意。也就是說,根據你的數據,如果你知道pa1==1那麼你可以肯定的是s1==0。因此,在這種情況下,擬合Cox模型並不合適,並且會導致數字錯誤。 這可以用

> coxph(Surv(t1, s1) ~ pa1, data=df1) 

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

看擬合模型前,這些跨表是非常重要可見一斑。在考慮那些涉及交互的情況之前,也應該從簡單的模型開始。

如果我們的交互項添加到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 

然後用

> (xtabs(~ s1+te1pa1, data=df1)) 
    te1pa1 
s1 0 
    0 3 
    1 3 

檢查,我們可以看到,這是一個沒用分類,即它不能幫助預測狀態s1

當組合所有3個術語時,儘管pe1是一個完美的預測器,但是鉗工確實能夠產生te1pe1的數值。然而,看一下係數的值和他們的錯誤表明他們是不合理的。

編輯 @JMarcelino:如果你從本例中的第一coxph模型的警告信息,你會看到警告信息:

2: In coxph(Surv(t1, s1) ~ te1 * pa1, data = df1) : 
    X matrix deemed to be singular; variable 3 

這可能是你同樣的錯誤得到和正由於這個分類問題。此外,您的第三個交叉表xtabs(~ tecnologia+pais, data=dados)不像statusinteraction term那麼重要。您可以首先手動添加交互項,如上例所示,然後檢查交叉表。或者你可以說:

> with(df1, 
     table(s1, pa1te1=pa1*te1)) 
    pa1te1 
s1 0 
    0 3 
    1 3 

這麼說,我發現自己的第三個表格單元格中的一個具有零(convPT),這意味着你有一個預測的這個組合沒有意見。這會在嘗試適應時導致問題。

在一般情況下,結果應該是有的預測各級一些值和預測不應結果歸類爲準確全有或全無50/50

編輯2 @ user75782131是,一般來說xtabs或類似的橫表應在模型中,其中結果和預測是離散的,即具有有限的沒有被執行。的水平。如果存在「完美分類」,那麼預測模型/迴歸可能就不合適。例如,邏輯迴歸(結果是二進制)以及Cox模型就是如此。

+3

非常感謝你的解釋,現在我知道什麼是「完美分類」。這個問題似乎被人詬病,現在我只得到警告:「X矩陣被認爲是單數」,這就是爲什麼我改變了問題的標題。這可能是由於高相關性嗎? – JMarcelino

+1

@dardisco,如果我理解正確,'xtabs'可以用來確定可以包含在公式'coxph()'中的哪些協變量? – 2014-01-20 07:58:14

相關問題