如何在 R 中编写多参数对数似然函数

Posted 2023-02-16

【中文标题】如何在 R 中编写多参数对数似然函数

【英文标题】：How to code a multiparameter log-likelihood function in R

【发布时间】：2014-01-14 04:37:58

【问题描述】：

我想估计以下问题的功效。我有兴趣比较两个都遵循 Weibull 分布的组。所以，A组有两个参数（形状par = a1，比例par = b1），两个参数有B组（a2，b2）。通过模拟感兴趣分布的随机变量（例如假设不同的尺度和形状参数，即 a1=1.5*a2 和 b1=b2*0.5；或者组之间的差异只是在形状或尺度参数中），应用 log-似然比检验来检验 a1=a2 和 b1=b2（或者例如 a1=a1，当我们知道 b1=b2 时），并估计检验的功效。

问题是完整模型的对数似然是什么，以及如何在 R 中编写代码 a) 有准确的数据， 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）？

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

请在下面找到确切数据的 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

【问题讨论】：

这个问题似乎离题了，因为它是关于如何得出对数似然的，因此不在 Stack Overflow 的范围内。它应该迁移到 stats.stackexchange.com。

这个问题可以通过一个小的改写来成为主题，例如“”

【参考方案1】：

是的，您可以将两组的对数似然相加（如果它们是单独计算的）。就像您将观察向量的对数似然相加一样，其中每个观察具有不同的生成参数。

我更喜欢用一个大向量（即形状参数）来考虑，它包含根据协变量结构（即组成员资格）而变化的值。在线性模型上下文中，该向量可以等于线性预测变量（一旦通过链接函数适当转换）：设计矩阵的点积和回归系数的向量。

这是一个（非功能化）示例：

## 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)。

区间删失数据

Weibull 分布（R 中的pweibull）的累积分布函数 F(T) 给出了时间 T 之前的故障概率。所以， 如果在时间 T[0] 和 T[1] 之间对观察进行间隔删失，则对象在 T[0] 之间失败的概率T[1] 是 F(T[1]) - F(T[0]) ： 对象在 T[1] 之前失败的概率减去它在 T[0] 之前失败的概率（T[0] 之间的 PDF 的积分 和 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)

【讨论】：

@user36478 添加了关于区间删失数据的可能性

是否可以扩展区间删失模型以适应时变协变量？