使用 mle2 / 最大似然估計的刪失二項式模型預測的 95% 置信區間

我正在研究一個問題，其中我有多對目前在世的男性i，每對都有一個假定的父系祖先ni幾代前（基於家譜證據），並且我有關於他們的 Y 染色體基因型是否不匹配的信息（僅限父系遺傳，xi= 1 表示不匹配，0 表示匹配）。如果沒有錯配，他們確實有一個共同的父系祖先，但如果有的話，一定是由於一個或多個婚外情而導致的鏈條中的扭結（我只能檢測到是否沒有或至少發生了一個這樣的額外配對親子關係事件，即因變量被刪失）。我感興趣的是獲得最大似然估計（加上 95% 的置信限），而不僅僅是平均外對親子關係 (EPP) 率（每一代孩子從外對交配中獲得的概率），而是還試圖推斷額外配對親子關係率可能如何隨時間而變化（因為分離共同父系祖先的世代的 nr 應該有一些關於此的信息 - 當存在不匹配時我不知道' 雖然不知道 EPP 什麼時候會發生，因為它可能發生在那個假定的祖先的世代和現在之間的任何地方，但是當有匹配時，我們確信在任何前幾代中都沒有 EPP）。因此，我的因二項式變量和獨立協變量生成/時間都被審查了。基於發布的有點類似的問題在這裡，我已經弄清楚瞭如何對人口和時間平均額外配對親子關係率phat加上 R 中的 95% 輪廓似然置信區間進行最大似然估計，如下所示：

# Function to make overall ML estimate of EPP rate p plus 95% profile likelihood confidence intervals, 
# taking into account that for pairs with mismatches multiple EPP events could have occured
#
# input is 
# x=vector of booleans or 0 and 1s specifying if there was a mismatch or not (1 or TRUE = mismatch)
# n=vector with nr of generations that separate common ancestor
# output is mle2 maximum likelihood fit with best-fit param EPP rate p

estimateP = function(x, n) {
 if (!is.logical(x[[1]])) x = (x==1)
 neglogL = function(p, x, n)  -sum((log(1 - (1-p)^n))[x]) -sum((n*log(1-p))[!x]) # negative log likelihood, see see https://stats.stackexchange.com/questions/152111/censored-binomial-model-log-likelihood
 require(bbmle)
 fit = mle2(neglogL, start=list(p=0.01), data=list(x=x, n=n))
 return(fit)
}

一些試點數據的示例（來自Larmuseau 等人 ProcB 2010）：

n = c(7, 7, 7, 7, 7, 8, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13, 13, 15, 15, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 18, 18, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 22, 23, 23, 24, 24, 25, 27, 31) # number of generations/meioses that separate presumed common paternal ancestor
x = c(rep(0,6), 1, rep(0,7), 1, 1, 1, 0, 1, rep(0,20), 1, rep(0,13), 1, 1, rep(0,5)) # whether pair of individuals had non-matching Y chromosomal genotypes

人口的最大似然估計和時間平均外對親子關係加 95% 置信限：

fit = estimateP(x,n)
c(coef(fit),confint(fit))*100 # estimated p and profile likelihood confidence intervals
# p 2.5 % 97.5 % 
# 0.9415172 0.4306652 1.7458847

即，所有孩子中有 0.9% [0.43-1.75% 95% CLs] 來自與假設不同的父親。

然後，我想更進一步，並嘗試估計p作為世代函數的外對親子關係的可能時間趨勢ni（為簡單起見，假設觀察到外對親子關係事件的對數機率和代），考慮到如果發生不匹配，EPP 事件可能發生在共同祖先ni和現在（第 0 代）之間的任何地方，並且如果沒有不匹配，則不會發生 EPP 事件那對特定個體的任何前幾代人。

如果在我們假設一個孩子是從一對額外的交配中衍生出來的概率之前是恆定的，如果是一個隨機變量，等於當觀察到 Y 染色體錯配時（對應於 1 個或多個 EPP 事件）和否則，則觀察到不匹配的概率（即，) 父系祖先在世時幾代前（） 曾是，而觀察到 EPP 事件的機會是

在獨立觀察的數據集中的與祖輩一起生活幾代人以前，因此可能性是

導致對數可能性

考慮到在我想要的包含時間動態的更複雜的模型中成為一個函數現在，與

， IE

然後，我相應地更改了上述似然函數的定義，並使用mle2package中的函數將其最大化bbmle：

# ML estimation, assuming that EPP rate p shows a temporal trend
# where log(p/(1-p))=a+b*n
# input is 
# x=vector of booleans or 0 and 1s specifying if there was a mismatch or not (1 or TRUE = mismatch)
# n=vector with nr of generations that separate common ancestor
# output is mle2 maximum likelihood fit with best-fit params a and b

estimatePtemp = function(x, n) {
 if (!is.logical(x[[1]])) x = (x==1)
 pfun = function(a, b, n) exp(a+b*n)/(1+exp(a+b*n)) # we now write p as a function of a, b and n
 logL = function(a, b, x, n)  sum((log(1 - (1-pfun(a, b, n))^n))[x]) + 
   sum((n*log(1-pfun(a, b, n)))[!x]) # a and b are params to be estimated, modified from https://stats.stackexchange.com/questions/152111/censored-binomial-model-log-likelihood
 neglogL = function(a, b, x, n)  -logL(a, b, x, n) # negative log-likelihood
 require(bbmle)
 fit = mle2(neglogL, start=list(a=-3, b=-0.1), data=list(x=x, n=n))
 return(fit)
}

# fitted coefficients
estfit = estimatePtemp(x, n)
cbind(coef(estfit),confint(estfit)) # parameter estimates and profile likelihood confidence intervals
# 2.5 % 97.5 %
# a -3.09054167 -5.3191406 -1.12078519
# b -0.09870851 -0.2396262 0.02848305
summary(estfit)
# Coefficients:
# Estimate Std. Error z value Pr(z) 
# a -3.090542 1.057382 -2.9228 0.003469 **
# b -0.098709 0.067361 -1.4654 0.142819 

這讓我對 EPP 費率的演變有一個合理的歷史估計隨著時間的推移 ：

pfun = function(a, b, n) exp(a+b*n)/(1+exp(a+b*n)) 
xvals=1:max(n)
par(mfrow=c(1,1))
plot(xvals,sapply(xvals,function (n) pfun(coef(estfit)[1], coef(estfit)[2], n)), 
    type="l", xlab="Generations ago (n)", ylab="EPP rate (p)")

但是，對於如何計算該模型整體預測的 95% 置信區間，我仍然有些困惑。有人會知道如何做到這一點嗎？也許使用總體預測區間（通過根據多元正態分佈的擬合重新採樣參數）（或者 delta 方法也可以工作？）？有人可以評論我上面的邏輯是否正確嗎？我還想知道這種審查二項式模型是否以文獻中的某個標準名稱而聞名，以及是否有人知道在這種模型下進行此類 ML 計算的任何已發表的工作？（我有一種感覺，問題應該是相當標準的，並且與已經完成的事情相對應，但似乎找不到任何東西……）

[有關此主題/問題的更多背景的 PS 論文可在此處和[此處獲得]](https://bio.kuleuven.be/ento/pdfs/larmuseau_etal_tree_2016.pdf)

根據Ben Bolker 的章節和上面的評論，我最終認為 95% 的人口預測區間由下式給出

# predicted EPP rate p as a function of n (nr of generations ago)
# plus 95% population prediction intervals (cf chapter B. Bolker)
pfun = function(a, b, n) exp(a+b*n)/(1+exp(a+b*n)) # model prediction
xvals=1:max(n)
set.seed(1001)
library(MASS)
nresamp=100000
pars.picked = mvrnorm(nresamp, mu = coef(estfit), Sigma = vcov(estfit)) # pick new parameter values by sampling from multivariate normal distribution based on fit
yvals = matrix(0, nrow = nresamp, ncol = length(xvals))
for (i in 1:nresamp) {
 yvals[i,] = sapply(xvals,function (n) pfun(pars.picked[i,1], pars.picked[i,2], n))
}
quant = function(col) quantile(col, c(0.025,0.975)) # 95% percentiles
conflims = apply(yvals,2,quant) # 95% confidence intervals

par(mfrow=c(1,1))
plot(xvals, sapply(xvals,function (n) pfun(coef(estfit)[1], coef(estfit)[2], n)), 
    type="l", xlab="Generations ago (n)", ylab="EPP rate (p)", ylim=c(0,0.05))
lines(xvals, conflims[1,], col="red" , lty=2)
lines(xvals, conflims[2,], col="red" , lty=2)

引用自：https://stats.stackexchange.com/questions/221426