Как кодировать многопараметрическую функцию логарифмического правдоподобия в R

Question

Я хотел бы оценить силу следующей проблемы. Я заинтересован в сравнении двух групп, которые следуют распределению Вейбулла. Таким образом, группа A имеет два параметра (форма par = a1, масштаб par = b1), а два параметра имеют группу B (a2, b2). Имитируя случайные величины из интересующего распределения (например, предполагая, что различные параметры масштаба и формы, например, a1=1,5*a2 и b1=b2*0,5; или различия между группами имеют только параметры формы или масштаба), применяют Проверка отношения правдоподобия для проверки, если a1=a2 и b1=b2 (или, например, a1=a1, когда мы знаем, что b1=b2), и оценка мощности теста.

Возникают вопросы: какова логарифмическая вероятность для полных моделей и как ее кодировать в 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

Answer 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)

Чтобы написать логарифмическую функцию правдоподобия, чтобы найти 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).

Данные с интервалом цензуры

Кумулятивная функция распределения F(T) распределения Вейбулла (pweibull в R) дает вероятность отказа до времени T. Таким образом, если наблюдение представляет собой интервал цензуры между моментами времени T[0] и T[1], вероятность того, что объект потерпит неудачу между T[0] и T[1], равна F(T[1]) - F(T[0] ]): вероятность того, что объект потерпел неудачу до T[1], минус вероятность того, что он потерпел неудачу до T[0] (интеграл PDF от T[0] и T[1]). Поскольку 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)