在给定事件发生之前,我正在使用生存分析来评估相对距离(而不是时间,通常是生存统计的情况)。由于我正在使用的数据集非常大,您可以下载我的数据集here的.rds文件
当使用survreg()
对相对距离建模时,我在模型摘要中遇到了NaN
和Inf
z和p值(可能来自Std Error
的0值):
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
s可能是什么原因?
@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()
协变量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