使用不同的分布拟合生存密度曲线
Fitting survival density curves using different distributions
我正在处理一些对数正态数据，自然我想证明对数正态分布的结果比其他可能的分布有更好的重叠。本质上，我想用我的数据复制下图：
其中拟合密度曲线并列在 log(time)。
链接图像所在的文本描述了拟合每个模型并获取以下参数的过程：
为此，我用上述分布拟合了四个朴素的生存模型：
survreg(Surv(time,event)~1,dist="family")
并提取形状参数(α)和系数(β)。
关于这个过程我有几个问题：
1) 这是正确的做法吗？我研究了几个 R 包，但找不到一个将密度曲线绘制为内置函数的包，所以我觉得我一定是忽略了一些明显的东西。
2) 对应对数正态分布（μ 和 σ$^2$）的值是否只是截距的均值和方差？
3) 如何在 R 中创建类似的 table？ （也许这更像是一个堆栈溢出问题）我知道我可以手动 cbind 它们，但我更感兴趣的是从拟合模型中调用它们。 survreg 对象存储系数估计值，但调用 survreg.obj$coefficients 会产生一个命名的数字向量（而不仅仅是一个数字）。
4) 最重要的是，如何绘制类似的图形？我认为如果我只提取参数并将它们绘制在直方图上会相当简单，但到目前为止运气不好。文本的作者说他根据参数估计了密度曲线，但我只是得到了一个点估计——我错过了什么？我应该在绘图之前根据分布手动计算密度曲线吗？
我不确定在这种情况下如何提供 mwe，但老实说，我只需要一个通用解决方案来为生存数据添加多条密度曲线。另一方面，如果您认为它会有所帮助，请随时推荐一个 mwe 解决方案，我会尝试制作一个。
感谢您的意见！
编辑：基于eclark的post，我取得了一些进展。我的参数是：
Dist = data.frame(
Exponential = rweibull(n = 10000, shape = 1, scale = 6.636684),
Weibull = rweibull(n = 10000, shape = 6.068786, scale = 2.002165),
Gamma = rgamma(n = 10000, shape = 768.1476, scale = 1433.986),
LogNormal = rlnorm(n = 10000, meanlog = 4.986, sdlog = .877)
)
然而，考虑到规模的巨大差异，这就是我得到的：
回到第3个问题，我应该这样获取参数吗？
目前我是这样做的（很抱歉弄得一团糟）：
summary(fit.exp)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "exponential")
Value Std. Error z p
(Intercept) 6.64 0.052 128 0
Scale fixed at 1
Exponential distribution
Loglik(model)= -2825.6 Loglik(intercept only)= -2825.6
Number of Newton-Raphson Iterations: 6
n= 397
summary(fit.wei)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "weibull")
Value Std. Error z p
(Intercept) 6.069 0.1075 56.5 0.00e+00
Log(scale) 0.694 0.0411 16.9 6.99e-64
Scale= 2
Weibull distribution
Loglik(model)= -2622.2 Loglik(intercept only)= -2622.2
Number of Newton-Raphson Iterations: 6
n= 397
summary(fit.gau)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "gaussian")
Value Std. Error z p
(Intercept) 768.15 72.6174 10.6 3.77e-26
Log(scale) 7.27 0.0372 195.4 0.00e+00
Scale= 1434
Gaussian distribution
Loglik(model)= -3243.7 Loglik(intercept only)= -3243.7
Number of Newton-Raphson Iterations: 4
n= 397
summary(fit.log)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "lognormal")
Value Std. Error z p
(Intercept) 4.986 0.1216 41.0 0.00e+00
Log(scale) 0.877 0.0373 23.5 1.71e-122
Scale= 2.4
Log Normal distribution
Loglik(model)= -2624 Loglik(intercept only)= -2624
Number of Newton-Raphson Iterations: 5
n= 397
我觉得我特别搞砸了对数正态分布，因为它不是标准的形状和系数串联，而是均值和方差。
试试这个；这个想法是使用随机分布函数生成随机变量，然后用输出数据绘制密度函数，这里有一个你需要的例子：
require(ggplot2)
require(dplyr)
require(tidyr)
SampleData <- data.frame(Duration=rlnorm(n = 184,meanlog = 2.859,sdlog = .246)) #Asume this is data we have sampled from a lognormal distribution
#Then we estimate the parameters for different types of distributions for that sample data and come up for this parameters
#We then generate a dataframe with those distributions and parameters
Dist = data.frame(
Weibull = rweibull(10000,shape = 1.995,scale = 22.386),
Gamma = rgamma(n = 10000,shape = 4.203,scale = 4.699),
LogNormal = rlnorm(n = 10000,meanlog = 2.859,sdlog = .246)
)
#We use gather to prepare the distribution data in a manner better suited for group plotting in ggplot2
Dist <- Dist %>% gather(Distribution,Duration)
#Create the plot that sample data as a histogram
G1 <- ggplot(SampleData,aes(x=Duration)) + geom_histogram(aes(,y=..density..),binwidth=5, colour="black", fill="white")
#Add the density distributions of the different distributions with the estimated parameters
G2 <- G1 + geom_density(aes(x=Duration,color=Distribution),data=Dist)
plot(G2)
我正在处理一些对数正态数据，自然我想证明对数正态分布的结果比其他可能的分布有更好的重叠。本质上，我想用我的数据复制下图：
其中拟合密度曲线并列在 log(time)。
链接图像所在的文本描述了拟合每个模型并获取以下参数的过程：
为此，我用上述分布拟合了四个朴素的生存模型：
survreg(Surv(time,event)~1,dist="family")
并提取形状参数(α)和系数(β)。
关于这个过程我有几个问题：
1) 这是正确的做法吗？我研究了几个 R 包，但找不到一个将密度曲线绘制为内置函数的包，所以我觉得我一定是忽略了一些明显的东西。
2) 对应对数正态分布（μ 和 σ$^2$）的值是否只是截距的均值和方差？
3) 如何在 R 中创建类似的 table？ （也许这更像是一个堆栈溢出问题）我知道我可以手动 cbind 它们，但我更感兴趣的是从拟合模型中调用它们。 survreg 对象存储系数估计值，但调用 survreg.obj$coefficients 会产生一个命名的数字向量（而不仅仅是一个数字）。
4) 最重要的是，如何绘制类似的图形？我认为如果我只提取参数并将它们绘制在直方图上会相当简单，但到目前为止运气不好。文本的作者说他根据参数估计了密度曲线，但我只是得到了一个点估计——我错过了什么？我应该在绘图之前根据分布手动计算密度曲线吗？
我不确定在这种情况下如何提供 mwe，但老实说，我只需要一个通用解决方案来为生存数据添加多条密度曲线。另一方面，如果您认为它会有所帮助，请随时推荐一个 mwe 解决方案，我会尝试制作一个。
感谢您的意见！
编辑：基于eclark的post，我取得了一些进展。我的参数是：
Dist = data.frame(
Exponential = rweibull(n = 10000, shape = 1, scale = 6.636684),
Weibull = rweibull(n = 10000, shape = 6.068786, scale = 2.002165),
Gamma = rgamma(n = 10000, shape = 768.1476, scale = 1433.986),
LogNormal = rlnorm(n = 10000, meanlog = 4.986, sdlog = .877)
)
然而，考虑到规模的巨大差异，这就是我得到的：
回到第3个问题，我应该这样获取参数吗？ 目前我是这样做的（很抱歉弄得一团糟）：
summary(fit.exp)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "exponential")
Value Std. Error z p
(Intercept) 6.64 0.052 128 0
Scale fixed at 1
Exponential distribution
Loglik(model)= -2825.6 Loglik(intercept only)= -2825.6
Number of Newton-Raphson Iterations: 6
n= 397
summary(fit.wei)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "weibull")
Value Std. Error z p
(Intercept) 6.069 0.1075 56.5 0.00e+00
Log(scale) 0.694 0.0411 16.9 6.99e-64
Scale= 2
Weibull distribution
Loglik(model)= -2622.2 Loglik(intercept only)= -2622.2
Number of Newton-Raphson Iterations: 6
n= 397
summary(fit.gau)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "gaussian")
Value Std. Error z p
(Intercept) 768.15 72.6174 10.6 3.77e-26
Log(scale) 7.27 0.0372 195.4 0.00e+00
Scale= 1434
Gaussian distribution
Loglik(model)= -3243.7 Loglik(intercept only)= -3243.7
Number of Newton-Raphson Iterations: 4
n= 397
summary(fit.log)
Call:
survreg(formula = Surv(duration, confterm) ~ 1, data = data.na,
dist = "lognormal")
Value Std. Error z p
(Intercept) 4.986 0.1216 41.0 0.00e+00
Log(scale) 0.877 0.0373 23.5 1.71e-122
Scale= 2.4
Log Normal distribution
Loglik(model)= -2624 Loglik(intercept only)= -2624
Number of Newton-Raphson Iterations: 5
n= 397
我觉得我特别搞砸了对数正态分布，因为它不是标准的形状和系数串联，而是均值和方差。
试试这个；这个想法是使用随机分布函数生成随机变量，然后用输出数据绘制密度函数，这里有一个你需要的例子：
require(ggplot2)
require(dplyr)
require(tidyr)
SampleData <- data.frame(Duration=rlnorm(n = 184,meanlog = 2.859,sdlog = .246)) #Asume this is data we have sampled from a lognormal distribution
#Then we estimate the parameters for different types of distributions for that sample data and come up for this parameters
#We then generate a dataframe with those distributions and parameters
Dist = data.frame(
Weibull = rweibull(10000,shape = 1.995,scale = 22.386),
Gamma = rgamma(n = 10000,shape = 4.203,scale = 4.699),
LogNormal = rlnorm(n = 10000,meanlog = 2.859,sdlog = .246)
)
#We use gather to prepare the distribution data in a manner better suited for group plotting in ggplot2
Dist <- Dist %>% gather(Distribution,Duration)
#Create the plot that sample data as a histogram
G1 <- ggplot(SampleData,aes(x=Duration)) + geom_histogram(aes(,y=..density..),binwidth=5, colour="black", fill="white")
#Add the density distributions of the different distributions with the estimated parameters
G2 <- G1 + geom_density(aes(x=Duration,color=Distribution),data=Dist)
plot(G2)