标签: weibull

问题

我正在尝试在JAGS中建立一个生存模型,允许时变协变量.我希望它是一个参数模型 - 例如,假设生存遵循威布尔分布(但我想让危险变化,所以指数太简单了).因此,这基本上是可以在flexsurv包中完成的贝叶斯版本,它允许参数模型中的时变协变量.

因此,我希望能够以"计数过程"形式输入数据,其中每个主题有多行,每行对应于其协变量保持不变的时间间隔(如本pdf或此处所述).这是包或包允许的(start, stop]配方.survivalflexurv

不幸的是,关于如何在JAGS中进行生存分析的每一个解释似乎都假设每个主题一行.

我试图采用这种更简单的方法并将其扩展到计数过程格式,但模型没有正确估计分布.

失败的尝试:

这是一个例子.首先,我们生成一些数据:

library('dplyr')
library('survival')

## Make the Data: -----
set.seed(3)
n_sub <- 1000
current_date <- 365*2

true_shape <- 2
true_scale <- 365

dat <- data_frame(person = 1:n_sub,
                  true_duration = rweibull(n = n_sub, shape = true_shape, scale = true_scale),
                  person_start_time = runif(n_sub, min= 0, max= true_scale*2),
                  person_censored = (person_start_time + true_duration) > current_date,
                  person_duration = ifelse(person_censored, current_date - person_start_time, true_duration)
)

  person person_start_time …

我试图重新创建最大似然分布拟合,我已经可以在Matlab和R中做到这一点,但现在我想使用scipy.特别是,我想估计我的数据集的Weibull分布参数.

我试过这个:

import scipy.stats as s
import numpy as np
import matplotlib.pyplot as plt

def weib(x,n,a):
    return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a)

data = np.loadtxt("stack_data.csv")

(loc, scale) = s.exponweib.fit_loc_scale(data, 1, 1)
print loc, scale

x = np.linspace(data.min(), data.max(), 1000)
plt.plot(x, weib(x, loc, scale))
plt.hist(data, data.max(), normed=True)
plt.show()

得到这个:

(2.5827280639441961, 3.4955032285727947)

并且看起来像这样的分布:

我一直exponweib在阅读http://www.johndcook.com/distributions_scipy.html.我也尝试了scipy中的其他Weibull函数(以防万一!).

在Matlab(使用分布拟合工具 - 参见屏幕截图)和R(使用MASS库函数fitdistr和GAMLSS包)中,我得到(loc)和b(比例)参数更像1.58463497 5.93030013.我相信所有三种方法都使用最大似然法进行分布拟合.

如果你想去,我已经在这里发布了我的数据!为了完整起见,我使用的是Python 2.7.5,Scipy 0.12.0,R 2.15.2和Matlab 2012b.

为什么我会得到不同的结果!？

我正在尝试将Weibull模型拟合并绘制成生存数据.该数据只有一个协变量,同期,从2006年到2010年.所以,任何关于如何添加到两行代码的想法,以绘制2010年队列的生存曲线？

library(survival)
s <- Surv(subSetCdm$dur,subSetCdm$event)
sWei <- survreg(s ~ cohort,dist='weibull',data=subSetCdm)

使用Cox PH模型完成相同操作非常简单,具有以下几行.问题是survfit()不接受类型为幸存的对象.

sCox <- coxph(s ~ cohort,data=subSetCdm)
cohort <- factor(c(2010),levels=2006:2010)
sfCox <- survfit(sCox,newdata=data.frame(cohort))
plot(sfCox,col='green')

使用数据肺(来自生存包),这是我想要完成的.

#create a Surv object
s <- with(lung,Surv(time,status))

#plot kaplan-meier estimate, per sex
fKM <- survfit(s ~ sex,data=lung)
plot(fKM)

#plot Cox PH survival curves, per sex
sCox <- coxph(s ~ as.factor(sex),data=lung)
lines(survfit(sCox,newdata=data.frame(sex=1)),col='green')
lines(survfit(sCox,newdata=data.frame(sex=2)),col='green')

#plot weibull survival curves, per sex, DOES NOT RUN
sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)
lines(survfit(sWei,newdata=data.frame(sex=1)),col='red')
lines(survfit(sWei,newdata=data.frame(sex=2)),col='red')

我正在构建一个R包来显示R中的Weibull图(使用graphics::plot).该图具有对数变换的x轴和Weibull变换的y轴(缺少更好的描述).因此,双参数威布尔分布可以在该图上表示为直线.

x轴的对数变换就像将log="x"参数添加到plot()或中一样简单curve().如何以优雅的方式提供y轴变换,以便所有与图形相关的绘图都可以在我的轴变换图上使用？要演示我需要的内容,请运行以下示例代码:

## initialisation ##
beta     <- 2;eta <- 1000
ticks    <- c(seq(0.01,0.09,0.01),(1:9)/10,seq(0.91,0.99,0.01))
F0inv    <- function (p) log(qweibull(p, 1, 1))
    # this is the transformation function
F0       <- function (q) exp(-exp(q))
    # this is the inverse of the transformation function
weibull  <- function(x)pweibull(x,beta,eta)
    # the curve of this function represents the weibull distribution 
    # as a straight line on weibull paper
weibull2 <- function(x)F0inv(weibull(x))

首先用威布尔分布的例子 …

我试图模拟一些遵循S形曲线关系的数据.在我的工作领域(心理物理学)中,Weibull函数通常用于模拟这种关系,而不是概率.

我正在尝试使用R创建一个模型,并且正在努力学习语法.我知道我需要使用包中的vglm()功能VGAM,但我无法得到一个合理的模型.这是我的数据:

# Data frame example data
dframe1 <- structure(list(independent_variable = c(0.3, 0.24, 0.23, 0.16, 
0.14, 0.05, 0.01, -0.1, -0.2), dependent_variable = c(1, 1, 
1, 0.95, 0.93, 0.65, 0.55, 0.5, 0.5)), .Names = c("independent_variable", 
"dependent_variable"), class = "data.frame", row.names = c(NA, 
-9L))

这是dframe1中的数据图:

library(ggplot2)

# Plot my original data
ggplot(dframe1, aes(independent_variable, dependent_variable)) + geom_point()

这应该能够通过Weibull函数建模,因为数据符合S形曲线关系.以下是我对数据建模并生成代表性图的尝试:

library(VGAM)

# Generate model
my_model <- vglm(formula = dependent_variable ~ independent_variable, family = weibull, data = dframe1)

# Create a new …

我有以下幸存模型:

Call:
survreg(formula = Surv(time = (ev.time), event = ev) ~ age, 
    data = my.data, dist = "weib")
             Value Std. Error    z        p
(Intercept) 4.0961     0.5566 7.36 1.86e-13
age         0.0388     0.0133 2.91 3.60e-03
Log(scale)  0.1421     0.1208 1.18 2.39e-01
Scale= 1.15 

Weibull distribution

我想根据上述估计绘制危险函数和生存函数.
我不想使用predict()或pweibull()(如此处所示参数生存或此处的问题.

我想用这个curve()功能.我有什么想法可以做到这一点？似乎幸存者的Weibull函数使用了比平常更多的尺度和形状定义(并且与例如rweibull不同).

更新:我想我真的需要它来表达危害/存活的估计的函数Intercept,age (+ other potential covariates),Scale不使用任何现成的*weilbull功能.

我试图使用从R中幸存的估计参数生成反Weibull分布.我的意思是我想,对于给定的概率(在MS Excel中实现的小型仿真模型中将是随机数),返回使用我的参数预计失败的时间.我理解逆Weibull分布的一般形式是:

X=b[-ln(1-rand())]^(1/a)

其中a和b分别是形状和比例参数,X是我想要的失败时间.我的问题在于解释来自幸存的截距和协变量参数.我有这些参数,时间单位是天:

             Value   Std. Error    z    p
(Intercept)     7.79    0.2288  34.051  0.000
Group 2        -0.139   0.2335  -0.596  0.551
Log(scale)     0.415    0.0279  14.88   0.000
Scale= 1.51 

Weibull distribution
Loglik(model)= -8356.7   Loglik(intercept only)= -8356.9 
Chisq = 0.37 on 1 degrees of freedom, p= 0.55 
Number of Newton-Raphson Iterations: 4 
n=1682 (3 observations deleted due to missing values)

我在帮助文件中读到R中的系数来自"极值分布",但我不确定这是什么意思以及我如何"回到"直接在公式中使用的标准比例参数.使用b = 7.79和a = 1.51给出了无意义的答案.我真的希望能够为基组和'Group 2'生成时间.我还应该注意,我自己没有进行分析,也无法进一步查询数据.

我的频率值随时间(x轴单位)而变化,如下图所示.在一些归一化之后,这些值可以被视为某些分布的密度函数的数据点.

问:假设这些频率点来自威布尔分布T,我如何将最佳威布尔密度函数拟合到这些点,以便从中推断出分布T参数？

sample <- c(7787,3056,2359,1759,1819,1189,1077,1080,985,622,648,518,
            611,1037,727,489,432,371,1125,69,595,624)

plot(1:length(sample), sample, type = "l")
points(1:length(sample), sample)

更新.为了防止被误解,我想补充一点解释.通过说我的频率值随时间变化(x轴单位)我的意思是我有数据说我有:

• 7787价值实现1
• 3056价值实现2
• 2359实现价值3 ......等

某种方式实现我的目标(我认为不正确)将创建一组这些实现:

# Loop to simulate values 
set.values <- c()
for(i in 1:length(sample)){
  set.values <<- c(set.values, rep(i, times = sample[i]))
}

hist(set.values)
lines(1:length(sample), sample)
points(1:length(sample), sample)

并使用fitdistr在set.values:

f2 <- fitdistr(set.values, 'weibull')
f2

为什么我认为这是不正确的方式以及为什么我在寻找更好的解决方案R？

• 在上面给出的分布拟合方法中,假设它set.values是从分布中完整的一组实现T

• 在我原来的问题中,我知道密度曲线第一部分的点 - …

我一直在R做一些数据分析,我试图找出如何使我的数据适合3参数Weibull分布.我找到了如何使用2参数Weibull来完成它,但是在找到如何使用3参数进行操作方面做得不够.

以下是我使用包中的fitdistr函数来拟合数据的方法MASS:

y <- fitdistr(x[[6]], 'weibull')

x[[6]] 是我的数据的子集,y是我存储拟合结果的地方.

我想估计以下问题的力量.我有兴趣比较两组都遵循Weibull分布.因此,A组有两个参数(形状par = a1,scale par = b1),两个参数有B组(a2,b2).通过模拟感兴趣分布的随机变量(例如假设不同的尺度和形状参数,即a1 = 1.5*a2,b1 = b2*0.5;或者组之间的差异只是形状或尺度参数),应用log-似然比测试以测试a1 = a2和b1 = b2(或者例如a1 = a1,当我们知道b1 = b2时),并估计测试的功效.

问题是什么是完整模型的对数似然,以及如何在a)具有精确数据时对R进行编码,以及b)对于区间删失数据？

也就是说,对于简化模型(当a1 = a2,b1 = b2时),精确和区间删失数据的对数似然是:

LL.reduced.exact <- function(par,data){sum(log(dweibull(data,shape=par[1],scale=par[2])))};
LL.reduced.interval.censored<-function(par, data.lower, data.upper) {sum(log((1-pweibull(data.lower, par[1], par[2])) – (1-pweibull(data.upper, par[1],par[2]))))}

什么是完整模型,当a1!= a2,b1!= b2时,考虑到两种不同的观测方案,即当必须估计4个参数时(或者,如果有兴趣查看形状参数的差异,必须估计3个参数)？

是否有可能估计它为不同的组建立两个对数似然并将它们加在一起(即LL.full <-LL.group1 + LL.group2)？

关于区间删失数据的对数似然,删失是非信息性的,所有观察都是间隔删失的.任何更好的想法如何执行此任务将不胜感激.

请在下面找到R代码以获取确切数据以说明问题.非常感谢你提前.

R Code:    
# n (sample size) = 500
# sim (number of simulations) = 1000
# alpha  = .05
# Parameters of Weibull distributions: 
   #group 1: a1=1, b1=20
   #group 2: a2=1*1.5 b2=b1

n=500 …