Нэн и Инф в резюме модели суррег - PullRequest
0 голосов
/ 07 января 2019

Я использую анализ выживания, чтобы оценить относительное расстояние (а не время, как это обычно бывает в статистике выживания) до того, как произошло данное событие. Поскольку набор данных, с которым я работаю, довольно большой, вы можете скачать файл .rds моего набора данных здесь

При моделировании относительного расстояния с использованием survreg() я обнаружил значения NaN и Inf z и p (предположительно полученные из 0 значений Std Error) в сводке модели:

Call:
survreg(formula = Surv(RelDistance, Status) ~ Backshore + LowerBSize + 
    I(LowerBSize^2) + I(LowerBSize^3) + State, data = DataLong, 
    dist = "exponential")
                            Value Std. Error        z         p
(Intercept)                   2.65469   1.16e-01  22.9212 2.85e-116
BackshoreDune                -0.08647   9.21e-02  -0.9387  3.48e-01
BackshoreForest / Tree (>3m) -0.07017   0.00e+00     -Inf  0.00e+00
BackshoreGrass - pasture     -0.79275   1.63e-01  -4.8588  1.18e-06
BackshoreGrass - tussock     -0.14687   1.00e-01  -1.4651  1.43e-01
BackshoreMangrove             0.08207   0.00e+00      Inf  0.00e+00
BackshoreSeawall             -0.47019   1.43e-01  -3.2889  1.01e-03
BackshoreShrub (<3m)         -0.14004   9.45e-02  -1.4815  1.38e-01
BackshoreUrban / Building     0.00000   0.00e+00      NaN       NaN
LowerBSize                   -0.96034   1.96e-02 -49.0700  0.00e+00
I(LowerBSize^2)               0.06403   1.87e-03  34.2782 1.66e-257
I(LowerBSize^3)              -0.00111   3.84e-05 -28.8070 1.75e-182
StateNT                       0.14936   0.00e+00      Inf  0.00e+00
StateQLD                      0.09533   1.02e-01   0.9384  3.48e-01
StateSA                       0.01030   1.06e-01   0.0973  9.22e-01
StateTAS                      0.19096   1.26e-01   1.5171  1.29e-01
StateVIC                     -0.07467   1.26e-01  -0.5917  5.54e-01
StateWA                       0.24800   9.07e-02   2.7335  6.27e-03

Scale fixed at 1 

Exponential distribution
Loglik(model)= -1423.4   Loglik(intercept only)= -3282.8
    Chisq= 3718.86 on 17 degrees of freedom, p= 0 
Number of Newton-Raphson Iterations: 6 
n= 6350 

Я думал, что Inf и NaN могут быть вызваны разделением данных, и объединил несколько уровней Backshore:

levels(DataLong$Backshore)[levels(DataLong$Backshore)%in%c("Grass - 
pasture", "Grass - tussock", "Shrub (<3m)")] <- "Grass - pasture & tussock 
/ Shrub(<3m)"
levels(DataLong$Backshore)[levels(DataLong$Backshore)%in%c("Seawall", 
"Urban / Building")] <- "Anthropogenic"
levels(DataLong$Backshore)[levels(DataLong$Backshore)%in%c("Forest / Tree 
(>3m)", "Mangrove")] <- "Tree(>3m) / Mangrove"

но проблема сохраняется при повторном запуске модели (т.е. Backshore Tree(>3m) / Mangrove).

Call:
survreg(formula = Surv(RelDistance, Status) ~ Backshore + LowerBSize + 
    I(LowerBSize^2) + I(LowerBSize^3) + State, data = DataLong, 
    dist = "exponential")
                                              Value Std. Error       z         p
(Intercept)                                      2.6684   1.18e-01  22.551 1.32e-112
BackshoreDune                                   -0.1323   9.43e-02  -1.402  1.61e-01
BackshoreTree(>3m) / Mangrove                   -0.0530   0.00e+00    -Inf  0.00e+00
BackshoreGrass - pasture & tussock / Shrub(<3m) -0.2273   8.95e-02  -2.540  1.11e-02
BackshoreAnthropogenic                          -0.5732   1.38e-01  -4.156  3.24e-05
LowerBSize                                      -0.9568   1.96e-02 -48.920  0.00e+00
I(LowerBSize^2)                                  0.0639   1.87e-03  34.167 7.53e-256
I(LowerBSize^3)                                 -0.0011   3.84e-05 -28.713 2.59e-181
StateNT                                          0.2892   0.00e+00     Inf  0.00e+00
StateQLD                                         0.0715   1.00e-01   0.713  4.76e-01
StateSA                                          0.0507   1.05e-01   0.482  6.30e-01
StateTAS                                         0.1990   1.26e-01   1.581  1.14e-01
StateVIC                                        -0.0604   1.26e-01  -0.479  6.32e-01
StateWA                                          0.2709   9.05e-02   2.994  2.76e-03

Scale fixed at 1 

Exponential distribution
Loglik(model)= -1428.4   Loglik(intercept only)= -3282.8
    Chisq= 3708.81 on 13 degrees of freedom, p= 0 
Number of Newton-Raphson Iterations: 6 
n= 6350 

Я искал объяснение этому поведению практически везде в документации к пакету survival и в Интернете, но не смог найти ничего, что связано с этим.

Кто-нибудь знает, что может быть причиной Inf и NaN с в этом случае?

Ответы [ 2 ]

0 голосов
/ 07 января 2019

Ковариата LowerBSize отлично предсказывает исход Status; Status==0 только если LowerBSize==0 и Status==1 только если LowerBSize>0.

table(DataLong$LowerBSize, DataLong$Status)

          0    1   
  0    4996    0
  1.2     0  271
  2.4     0  331
  4.9     0  256
  9.6     0  155
  19.2    0  148
  36.3    0  193

Удобный способ рассмотреть LowerBSize в модели - включить двоичную переменную LowerBSize>0:

survreg(formula = Surv(RelDistance, Status) ~ Backshore +  State + 
        I(LowerBSize>0), data = DataLong,  dist = "exponential")

Coefficients:
                 (Intercept)                BackshoreDune BackshoreForest / Tree (>3m) 
                 22.97248461                  -0.04798348                  -0.27440059 
    BackshoreGrass - pasture     BackshoreGrass - tussock            BackshoreMangrove 
                 -0.33624746                  -0.07545700                   0.12020217 
            BackshoreSeawall         BackshoreShrub (<3m)    BackshoreUrban / Building 
                 -0.01008893                  -0.05115076                   0.29125024 
                     StateNT                     StateQLD                      StateSA 
                  0.15385826                   0.11617931                   0.08405151 
                    StateTAS                     StateVIC                      StateWA 
                  0.14914393                   0.08803225                   0.06395311 
       I(LowerBSize > 0)TRUE 
                -23.75967069 

Scale fixed at 1 

Loglik(model)= -316.5   Loglik(intercept only)= -3282.8
        Chisq= 5932.66 on 15 degrees of freedom, p= <2e-16 
n= 6350
0 голосов
/ 07 января 2019

@ MarcoSandri прав, что цензура путается с LowerBSize, но я не уверен, что это все решение. Это могло бы объяснить, почему модель настолько нестабильна, но сама по себе не должна (AFAICT) сделать модель некорректной. Если я заменю LowerBSize+ I(LowerBSize^2) + I(LowerBSize^3) на ортогональный полином (poly(LowerBSize,3)), я получу более разумные ответы:

ss3 <- survreg(formula = Surv(RelDistance, Status) ~ Backshore +
                   poly(LowerBSize,3) + State, data = DataLong, 
               dist = "exponential")

Call:
survreg(formula = Surv(RelDistance, Status) ~ Backshore + poly(LowerBSize, 
    3) + State, data = DataLong, dist = "exponential")
                                 Value Std. Error      z       p
(Intercept)                   2.18e+00   1.34e-01  16.28 < 2e-16
BackshoreDune                -1.56e-01   1.06e-01  -1.47 0.14257
BackshoreForest / Tree (>3m) -2.24e-01   2.01e-01  -1.11 0.26549
BackshoreGrass - pasture     -8.63e-01   1.74e-01  -4.97 6.7e-07
BackshoreGrass - tussock     -2.14e-01   1.13e-01  -1.89 0.05829
BackshoreMangrove             3.66e-01   4.59e-01   0.80 0.42519
BackshoreSeawall             -5.37e-01   1.53e-01  -3.51 0.00045
BackshoreShrub (<3m)         -2.08e-01   1.08e-01  -1.92 0.05480
BackshoreUrban / Building    -1.17e+00   3.22e-01  -3.64 0.00028
poly(LowerBSize, 3)1         -6.58e+01   1.41e+00 -46.63 < 2e-16
poly(LowerBSize, 3)2          5.09e+01   1.19e+00  42.72 < 2e-16
poly(LowerBSize, 3)3         -4.05e+01   1.41e+00 -28.73 < 2e-16
StateNT                       2.61e-01   1.93e-01   1.35 0.17557
StateQLD                      9.72e-02   1.12e-01   0.87 0.38452
StateSA                      -4.11e-04   1.15e-01   0.00 0.99715
StateTAS                      1.91e-01   1.35e-01   1.42 0.15581
StateVIC                     -9.55e-02   1.35e-01  -0.71 0.47866
StateWA                       2.46e-01   1.01e-01   2.44 0.01463

Если я подхожу к той же модели, но с poly(LowerBSize,3,raw=TRUE) (вызывающий результат ss4, см. Ниже), я снова получу ваши патологии. Кроме того, модель с ортогональными полиномами на самом деле подходит лучше (имеет более высокое логарифмическое правдоподобие):

logLik(ss4)
## 'log Lik.' -1423.382 (df=18)
logLik(ss3)
## 'log Lik.' -1417.671 (df=18)

В идеальном математическом / вычислительном мире это не должно быть правдой - это еще один признак того, что что-то нестабильно в определении эффектов LowerBSize таким образом. Я немного удивлен, что это происходит - число уникальных значений LowerBSize мало, но не должно быть патологическим, а диапазон значений не огромен или далек от нуля ...


Я до сих пор не могу сказать, что на самом деле вызывает это, но проксимальная проблема, вероятно, заключается в сильной корреляции между линейными / квадратичными / кубическими терминами: для чего-то более простого, например, линейной регрессии, корреляция 0,993 (между квадратическими и кубическими терминами) не Это не вызывает серьезных проблем, но чем сложнее численная проблема (например, анализ выживаемости или линейная регрессия), тем больше корреляции может быть проблемой ...

X <- model.matrix( ~ Backshore + LowerBSize + 
                       I(LowerBSize^2) + I(LowerBSize^3) + State,
                  data=DataLong)
print(cor(X[,grep("LowerBSize",colnames(X))]),digits=3)
library(corrplot)
png("survcorr.png")
corrplot.mixed(cor(X[,-1]),lower="ellipse",upper="number",
             tl.cex=0.4)
dev.off()

enter image description here

...