我想估计以下问题的力量。我有兴趣比较两个都遵循Weibull分布的两组。因此,组A具有两个参数(形状par = a1,比例par = b1),并且两个参数具有组B(a2,b2)。通过从感兴趣的分布中模拟随机变量(例如,假设使用不同的比例尺和形状参数,即a1 = 1.5 * a2和b1 = b2 * 0.5;或者组之间的差异只是形状或比例尺参数),应用对数似然比检验以检验a1 = a2和b1 = b2(或例如,当我们知道b1 = b2时a1 = a1),并估计检验的功效。
问题是完整模型的对数似然性是什么,以及当a)拥有准确数据,b)间隔检查的数据时,如何在R中对其进行编码?
也就是说,对于简化模型(当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
sim=1000
alpha=.05
a1=1
b1=20
a2=a1*1.5
b2=b1
#OR: a1=1, b1=20, a2=a1*1.5, b2=b1*0.5
# the main question is how to build this log-likelihood model, when a1!=a2, and b1=b2
# (or a1!=a2, and b1!=b2)
LL.full<-?????
LL.reduced <- function(par,data){sum(log(dweibull(data,shape=par[1],scale=par[2])))}
LR.test<-function(red,full,df) {
lrt<-(-2)*(red-full)
pvalue<-1-pchisq(lrt,df)
return(data.frame(lrt,pvalue))
}
rejections<-NULL
for (i in 1:sim) {
RV1<-rweibull (n, a1, b1)
RV2<-rweibull (n, a2, b2)
RV.Total<-c(RV1, RV2)
par.start<-c(1, 15)
mle.full<- ????????????
mle.reduced<-optim(par.start, LL, data=RV.Total, control=list(fnscale=-1))
LL.full<-?????
LL.reduced<-mle.reduced$value
LRT<-LR.test(LL.reduced, LL.full, 1)
rejections1<-ifelse(LRT$pvalue<alpha,1,0)
rejections<-c(rejections, rejections1)
}
table(rejections)
sum(table(rejections)[[2]])/sim # estimated power
是的,您可以将两组的对数似然相加(如果它们是分别计算的)。就像您将观察向量的对数似然相加一样,其中每个观察具有不同的生成参数。
我更愿意考虑一个大向量(即形状参数),该向量包含根据协变量的结构(即组成员)而变化的值。在线性模型上下文中,此向量可以等于线性预测变量(一旦通过链接函数进行了适当的转换):设计矩阵的点积与回归系数的向量。
这是一个(未功能化的)示例:
## setup true values
nobs = 50 ## number of observations
a1 = 1 ## shape for first group
b1 = 2 ## scale parameter for both groups
beta = c(a1, a1 * 1.5) ## vector of linear coefficients (group shapes)
## model matrix for full, null models
mm_full = cbind(grp1 = rep(c(1,0), each = nobs), grp2 = rep(c(0,1), each = nobs))
mm_null = cbind(grp1 = rep(1, nobs*2))
## shape parameter vector for the full, null models
shapes_full = mm_full %*% beta ## different shape parameters by group (full model)
shapes_null = mm_null %*% beta[1] ## same shape parameter for all obs
scales = rep(b1, length(shapes_full)) ## scale parameters the same for both groups
## simulate response from full model
response = rweibull(length(shapes_full), shapes_full, scales)
## the log likelihood for the full, null models:
LL_full = sum(dweibull(response, shapes_full, scales, log = T))
LL_null = sum(dweibull(response, shapes_null, scales, log = T))
## likelihood ratio test
LR_test = function(LL_null, LL_full, df) {
LR = -2 * (LL_null - LL_full) ## test statistic
pchisq(LR, df = df, ncp = 0, lower = F) ## probability of test statistic under central chi-sq distribution
}
LR_test(LL_null, LL_full, 1) ## 1 degrees freedom (1 parameter added)
要编写对数似然函数以找到Weibull模型的MLE,其中形状参数是协变量的某些线性函数,可以使用相同的方法:
## (negative) log-likelihood function
LL_weibull = function(par, data, mm, inv_link_fun = function(.) .){
P = ncol(mm) ## number of regression coefficients
N = nrow(mm) ## number of observations
shapes = inv_link_fun(mm %*% par[1:P]) ## shape vector (possibly transformed)
scales = rep(par[P+1], N) ## scale vector
-sum(dweibull(data, shape = shapes, scale = scales, log = T)) ## negative log likelihood
}
然后,您的功耗模拟可能如下所示:
## function to simulate data, perform LRT
weibull_sim = function(true_shapes, true_scales, mm_full, mm_null){
## simulate response
response = rweibull(length(true_shapes), true_shapes, true_scales)
## find MLE
mle_full = optim(par = rep(1, ncol(mm_full)+1), fn = LL_weibull, data = response, mm = mm_full)
mle_null = optim(par = rep(1, ncol(mm_null)+1), fn = LL_weibull, data = response, mm = mm_null)
## likelihood ratio test
df = ncol(mm_full) - ncol(mm_null)
return(LR_test(-mle_null$value, -mle_full$value, df))
}
## run simulations
nsim = 1000
pvals = sapply(1:nsim, function(.) weibull_sim(shapes_full, scales, mm_full, mm_null) )
## calculate power
alpha = 0.05
power = sum(pvals < alpha) / nsim
身份链接在上面的示例中工作正常,但是对于更复杂的模型,可能需要某种转换。
而且,您不必在对数似然函数中使用线性代数-显然,您可以按照您认为合适的任何方式构造形状的向量(只要您在向量中明确索引了适当的生成参数par
)。
间隔检查的数据
威布尔分布(在R中)的累积分布函数F(Tpweibull
)给出了在时间T之前发生故障的概率。因此,如果观察是在时间T [0]和T [1]之间进行间隔检查的,则对象在T [0]和T [1]之间失效的概率为F(T [1])-F(T [0 ]):对象在T [1]之前失败的概率减去它在T [0]之前失败的概率(PDF在T [0]和T [1]之间的积分)。由于Weibull CDF已在R中实现,因此修改上面的似然函数很简单:
LL_ic_weibull <- function(par, data, mm){
## 'data' has two columns, left and right times of censoring interval
P = ncol(mm) ## number of regression coefficients
shapes = mm %*% par[1:P]
scales = par[P+1]
-sum(log(pweibull(data[,2], shape = shapes, scale = scales) - pweibull(data[,1], shape = shapes, scale = scales)))
}
或者,如果您不想使用模型矩阵等,而只是限制自己只能按组索引形状参数矢量,则可以执行以下操作:
LL_ic_weibull2 <- function(par, data, nobs){
## 'data' has two columns, left and right times of censoring interval
## 'nobs' is a vector that contains the num. observations for each group (grp1, grp2, ...)
P = length(nobs) ## number of regression coefficients
shapes = rep(par[1:P], nobs)
scales = par[P+1]
-sum(log(pweibull(data[,2], shape = shapes, scale = scales) - pweibull(data[,1], shape = shapes, scale = scales)))
}
测试两个函数是否提供相同的解决方案:
## generate intervals from simulated response (above)
left = ifelse(response - 0.2 < 0, 0, response - 0.2)
right = response + 0.2
response_ic = cbind(left, right)
## find MLE w/ first LL function (model matrix)
mle_ic_full = optim(par = c(1,1,3), fn = LL_ic_weibull, data = response_ic, mm = mm_full)
mle_ic_null = optim(par = c(1,3), fn = LL_ic_weibull, data = response_ic, mm = mm_null)
## find MLE w/ second LL function (groups only)
nobs_per_group = apply(mm_full, 2, sum) ## just contains number of observations per group
nobs_one_group = nrow(mm_null) ## one group so only one value
mle_ic_full2 = optim(par = c(1,1,3), fn = LL_ic_weibull2, data = response_ic, nobs = nobs_per_group)
mle_ic_null2 = optim(par = c(1,3), fn = LL_ic_weibull2, data = response_ic, nobs = nobs_one_group)
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句