4

我正在使用生存分析来评估给定事件发生之前的相对距离（而不是时间，因为它通常是生存统计的情况）。由于我正在使用的数据集非常大，您可以在此处下载我的数据集的 .rds 文件

在使用 对相对距离进行建模时，我在模型摘要中survreg()遇到了zNaN和Infp 值（可能源自 的 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

我认为InfandNaN可能是由数据分离引起的，并将某些级别合并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和NaNs 的原因是什么？

4

2 回答 2

3

协变量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
于 2019-01-07T01:16:58.940 回答
2

@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()

在此处输入图像描述

于 2019-01-07T01:07:03.330 回答