次の問題の検出力を計算したいと思います。ワイブル分布に従う 2 つのグループを比較することに興味があります。したがって、グループ A には 2 つのパラメータ (形状 par = a1、スケール par = b1) があり、2 つのパラメータにはグループ 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 の場合、2 つの異なる観測スキームを考慮した場合、つまり 4 つのパラメーターを推定する必要がある場合 (または、形状パラメーターの違いを調べたい場合は、 3 つのパラメータを推定する必要があります)?
別々のグループに対して 2 つの対数尤度を構築し、それを合計して推定することは可能ですか (つまり、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