Comment on “Bayesian additional evidence for decision making under small sample uncertainty”

 ## reproduce BAE calculations from Sondhi et al. (2021)

## -----------------------------------------------------------------------------

## data from Sondhi et al. (2021)
hr <- 0.42
hrCI <- c(0.14, 1.23)
x <- log(hr)
za <- qnorm(p = 0.975)
se <- (log(hrCI)[2] - log(hrCI)[1])/(2*za)
p <- 2*pnorm(q = abs(x)/se, lower.tail = FALSE)

## prior evidence from Innocenti et al. (2019)
hr2 <- 0.13
hr2CI <- c(0.06, 0.30)
x2 <- log(0.13)
se2 <- (log(hr2CI[2]) - log(hr2CI[1]))/(2*za)

## Bayesian additional evidence based on effect size and standard error
baeES <- function(x, se, alpha) {
    bae <- sign(x)*sqrt(2)*qnorm(p = 1 - alpha/2)*se - x
    return(bae)
}

## Bayesian additional evidence based on confidence interval limits
baeCI <- function(L, U) {
    x <- mean(c(L, U))
    bae <- (sign(x)*sqrt(2)*(U - L) - (U + L))/2
    return(bae)
}

## compute BAE for data from Sondhi et al. (2021)
bae <- baeCI(L = log(hrCI)[1], U = log(hrCI)[2])


## functions to compute prior parameters of the 3 reverse-Bayes approaches
## -----------------------------------------------------------------------------

## determine advocacy prior such that (1 - alpha) posterior credible interval
## fixed to zero and prior coefficient of variation fixed to 1/z_(alpha/2)
advPriorCV <- function(x, se, alpha) {
    ## x: data
    ## se: standard error of data
    ## alpha: significance level
    z <- x/se
    za <- qnorm(p = 1 - alpha/2)
    if (abs(z) >= za) {
        out <- out <- c("mean" = NaN, "var" = NaN)
    } else {
        meanPrior <- 2*x/(1 - z^2/za^2)
        varPrior <- meanPrior^2/za^2
        out <- c("mean" = meanPrior, "var" = varPrior)
    }
    return(out)
}

## determine advocacy prior such that (1 - alpha) posterior credible interval
## fixed to zero and prior variance fixed to variance of data
advPriorVar <- function(x, se, alpha) {
    ## x: data
    ## se: standard error of data
    ## alpha: significance level
    z <- x/se
    za <- qnorm(p = 1 - alpha/2)
    if (abs(z) >= za) {
        out <- out <- c("mean" = NaN, "var" = NaN)
    } else {
        meanPrior <- sign(x)*sqrt(2)*za*se - x
        varPrior <- se^2
        out <- c("mean" = meanPrior, "var" = varPrior)
    }
    return(out)
}

## determine advocacy prior such that (1 - alpha) posterior credible interval
## fixed to zero and prior mean fixed to value of data
advPriorMean <- function(x, se, alpha) {
    ## x: data
    ## se: standard error of data
    ## alpha: significance level
    z <- x/se
    za <- qnorm(p = 1 - alpha/2)
    if (abs(z) >= za) {
        out <- out <- c("mean" = NaN, "var" = NaN)
    } else {
        meanPrior <- x
        varPrior <- se^2/(za^2/z^2 - 1)
        out <- c("mean" = meanPrior, "var" = varPrior)
    }
    return(out)
}

## compute posterior credible interval to check correctness of functions
postCrI <- function(x, se, meanPrior, varPrior, alpha) {
    postVar <- 1/(1/se^2 + 1/varPrior)
    postMean <- (x/se^2 + meanPrior/varPrior)*postVar
    za <- qnorm(p = 1 - alpha/2)
    crI <- postMean + c(-1, 1)*za*sqrt(postVar)
    names(crI) <- c("lower", "upper")
    return(crI)
}


## create similar figure as in Sondhi et al. (2021)
## -----------------------------------------------------------------------------

## significance level
alpha <- 0.05

## bae
bae <- baeES(x = x, se = se, alpha = alpha)
bae <- baeCI(L = log(hrCI[1]), U = log(hrCI[2]))
baePost <- postCrI(x = x, se = se, meanPrior = bae, varPrior = se^2,
                   alpha = alpha)

## AnCred
anCredPars <- advPriorCV(x = x, se = se, alpha = alpha)
anCredPost <- postCrI(x = x, se = se, meanPrior = anCredPars[1],
                      varPrior = anCredPars[2], alpha = alpha)

## fixed mean
fmPars <- advPriorMean(x = x, se = se, alpha = alpha)
fmPost <- postCrI(x = x, se = se, meanPrior = fmPars[1],
                  varPrior = fmPars[2], alpha = alpha)

## CIs
lower <- c(x - za*se,
           bae - za*se, anCredPars[1] - za*sqrt(anCredPars[2]),
           fmPars[1] - za*sqrt(fmPars[2]),
           baePost[1], anCredPost[1], fmPost[1],
           x2 - za*se2)
est <- c(x,
         bae, anCredPars[1], fmPars[1],
         mean(baePost), mean(anCredPost), mean(fmPost),
         x2)
upper <- c(x + za*se,
           bae + za*se, anCredPars[1] + za*sqrt(anCredPars[2]),
           fmPars[1] + za*sqrt(fmPars[2]),
           baePost[2], anCredPost[2], fmPost[2],
           x2 + za*se2)
colpalette <- c("#000000", "#1B9E77", "#D95F02", "#7570B3", "#000000")
cols <- colpalette[c(1, 2, 3, 4, 2, 3, 4, 5)]


xseq <- c(1, 2.8, 3, 3.2, 1.8, 2, 2.2, 4)
xbks <- c(1, 2, 3, 4)
xlab <- c("Data", "Posterior", "Advocacy Prior", "Additional data")
ybks <- c(log(10), log(1), log(1/10), log(1/100), log(1/1000))
plot(x = xseq, y = est, type = "n", ylim = c(min(lower), log(5)),
     xlim = c(0.8*min(xbks), 1.1*max(xbks)), xaxt = "n", yaxt = "n",
     xlab = "", ylab = "", las = 2)
mtext(side = 2, text = "Hazard Ratio", line = 3.3)
abline(h = 0, lty = 2, col = adjustcolor(col = 1, alpha = 0.5), lwd = 1.5)
axis(side = 2, at = ybks, labels = pCalibrate::formatBF(exp(ybks)), las = 2)
axis(side = 1, at = xbks, labels = xlab)
abline(h = ybks, lty = 1, col = adjustcolor(col = 1, alpha = 0.1))
box()
arrows(x0 = xseq, y0 = lower, y1 = upper, angle = 90, code = 3, length = 0.1,
       lwd = 1.2, col = cols)
points(x = xseq, y = est, pch = 20, cex = 1.2, col = cols)
legend("bottomright", c("BAE", "AnCred", "fixed mean"), pch = 20,
       col = colpalette[2:4], bg = "white")


## plot partition of the prior parameter space
## -----------------------------------------------------------------------------

## BAE for different prior variances than the one observed in the data
bae2 <- function(x, se, alpha, g) {
    ## g: the relative prior variance
    bae <- (sign(x)*qnorm(p = 1 - alpha/2)*se*sqrt(g/(1 + g)) - x*g/(1 + g))*(1 + g)
    return(bae)
}

## prior parameter space data
priorsDF <- data.frame(mean = c(anCredPars[1], bae, fmPars[1]),
                       var = c(anCredPars[2], se^2, fmPars[2]),
                       type = c("AnCred", "BAE", "fixed mean"),
                       color = colpalette[c(3, 2, 4)])

## plot prior parameter space
loggseq <- seq(log(0.01), log(25), length.out = 500)
gbks <- c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20)
baeseq <- bae2(x = x, se = se, alpha = alpha, g = exp(loggseq))
plot(x = loggseq, y = baeseq/x, type = "l", las = 1, xaxt = "n",
     xlab = bquote("relative prior variance" ~ italic(g) == tau^2/sigma^2),
     ylab = "", xlim = log(c(0.1, 20)), ylim = c(0, 6), lwd = 1.8)
axis(side = 1, at = log(gbks), labels = gbks)
mtext(side = 2, text =bquote("relative prior mean" ~ mu/hat(theta)),
      line = 2)
## draw non-credibiliy region
polygon(x = c(head(loggseq)[1], loggseq, tail(loggseq)[1]),
        y = c(10, baeseq, 10)/x, density = 10,
        angle = 135, border = NA, col = adjustcolor(col = "black", alpha = 0.3),
        lty = 2, lwd = 1)
## draw credibiliy region
polygon(x = c(head(loggseq)[1], loggseq, tail(loggseq)[1]),
        y = c(-100, baeseq, -100)/x, density = 10,
        angle = 45, border = NA, col = adjustcolor(col = "black", alpha = 0.3),
        lty = 3, lwd = 1.5)
## draw significance region
lines(x = loggseq, y = -qnorm(p = 1 - alpha/2)*sqrt(exp(loggseq)*se^2)/x,
      col = adjustcolor(colpalette[3], alpha = 0.5))
## draw fixed mean/variance
abline(h = 1, col = adjustcolor(colpalette[4], alpha = 0.5))
abline(v = 0, col = adjustcolor(colpalette[2], alpha = 0.5))
## draw tipping points
points(x = log(priorsDF[,2]/se^2), y = priorsDF[,1]/x, col = priorsDF[,4],
       pch = 20, cex = 1.6)
legend("topleft", c("BAE", "AnCred", "fixed mean"), pch = 20,
       col = colpalette[2:4], bg = "white")
text(x = log(0.28), y = 1.5,
     labels = bquote("posterior credible at" ~ alpha == .(alpha)), cex = 1.1)
text(x = log(5), y = 0.5,
     labels = bquote("posterior not credible at" ~ alpha == .(alpha)), cex = 1.1)

留言

張貼留言

這個網誌中的熱門文章

可轉移性、普遍性、代表性和外部有效性

 ### 可轉移性（Transportability） - **定義：** 可轉移性指的是將一個研究中得到的因果效應估計從一個人群（研究人群）應用到另一個不同的人群（目標人群）的能力。它涉及特定的統計方法來調整人群之間的差異，並特別關注為了合理和實施此類應用所需的假設和方法。 - **用途：** 當研究人員想知道在一個設置或人群中觀察到的干預效果是否在不同的設置或人群中保持時使用。 ### 普遍性（Generalizability） - **定義：** 普遍性是指研究發現可以應用到更廣泛人群、設置、時間或未直接研究的測量指標的程度。它是一個比可轉移性更廣泛的概念，更多地關注研究結果在特定研究條件之外的適用性。 - **用途：** 在評估研究結果是否適用於明確研究之外的更廣泛背景時考慮，而不必然採用特定的統計方法來進行此類應用。 ### 代表性（Representativeness） - **定義：** 代表性指的是研究樣本如何緊密地反映其所來自的總體或目標人群的特徵。它涉及抽樣過程和樣本在多大程度上體現了更廣泛人群的特徵。 - **用途：** 這一概念在研究設計階段至關重要，影響研究結果能夠被泛化到總體的程度。 ### 外部有效性（External Validity） - **定義：** 外部有效性是指研究結果可以被泛化或應用到原始研究設置、人群、治療變量和結果之外的背景的程度。它包括普遍性和可轉移性的概念，但範圍更廣，包括隨時間和在不同條件下應用研究發現的適用性。 - **用途：** 它涉及到研究發現超出特定研究範圍的整體適用性和相關性，考慮諸如時間、設置和人群等因素。 **關鍵差異和關係：** - **可轉移性與普遍性：** 可轉移性是普遍性的一個子集，特別關注使用統計調整將結果應用到不同人群的能力。普遍性是一個更廣泛的概念，不必然涉及此類調整。 - **代表性與外部有效性：** 代表性樣本可以增強研究的外部有效性，但外部有效性也取決於其他因素，如研究設計和執行。 - **外部有效性包含：** 普遍性和可轉移性，以及其他考慮因素，以便將研究發現應用於原始背景之外。 總之，雖然所有這些概念都涉及到將研究發現應用於原始研究條件之外，但它們每一個都關注此應用的不同方面，從用於調整不同人群之間差異的統計方法（可轉移性）到研究發現的整體適用性（外部有效性）。 - ...
閱讀完整內容

頻率學派 vs 貝氏學派

機率是重複事件長期發生的頻率 vs 主觀上某事件發生的不確定性（比較像是一般人的觀念）。 傳統 vs 新 簡單vs 複雜 人工計算可行 vs 不可行 母數（參數）是未知的常數（沒有機率可言） vs 未知的隨機變數（有機率） 樣本是隨機的變數 vs 已經發生的固定數 客觀（只依據本次試驗的結果）vs 主觀（依據先驗信念及本次試驗的結果而更新成後驗信念）。 H: 假設，D: 資料。H0: 虛無假設（參數或母數 θ 等於 0），H1: 對立假設，p(H1)=1-p(H0)。 一個棒球比賽的觀眾關心的是「機率」p(D|H)，亦即預測「這個左打者擊中球的機率是多少？」；當他看見一名打者擊中了球以後，計算這一名打者是一個左打者的機率是多少則是屬於「逆機率」p(H|D)。 醫生關心的是「逆機率」p(H|D)：當看見某症候、檢驗、檢查 D時，用貝氏定理計算病人有某疾病 H 的機率。 估計 p 值（H0 為真時，看見 D 或比它更極端的機率 p(D|H0)) vs 估計後驗機率 p(H1|D)：p(D|H) ≠ p(H|D)。 95% 信賴區間包含母數的機率是 0% 或 100% vs 95% 置信區間包含母數的機率是 95%。 頻率學派：當 p < α（例如 0.05）時拒絕 H0，當 p ≥ α 時無法拒絕 H0。但是我們無法知道 p(H0|D)、p(H1|D)，而且無論 p 值多大都無法接受 H0。 貝氏定理： 1. H1 的後驗勝算 p(H1|D)/(p(H0|D))= 先驗勝算 x 似然比或貝氏因子 p(D|H1)/p(D|H0)，勝算=機率/(1-機率），機率=勝算/(勝算+1)。 2. 後驗機率 p(H|D）= 先驗機率 p(H) x p(D/H)/p(D) p(D): p(D/H0) + p(D/H1) RCT 時，頻率學派的 α 是 H0 為真（無效）時，你卻斷言為假（有效）的機率：p(斷言有效|無效)。傳統上它被稱為「第一型錯誤率」，其實它應該是「第一型斷言率」，而不是治療無效你卻以為有效的錯誤機率：貝氏學派的 1 - p(有效|資料，先驗機率）。後者才是我們所關心的。 ESP and the significance of significance 22 https://understandinguncertainty.org/node/1286 爲什麼我...
閱讀完整內容

貝氏分析計算器

貝氏分析計算器 ：20250130 by 顧進裕, Claude-3.5-Sonnet 平均值與標準差計算器 貝氏定理：後驗分布 p(H|D) = 先驗分布 p(H) × 似然分布 p(D|H) / 邊際分布 p(D) 似然分布：觀察到的平均值與標準誤 se 後驗分布：95% 可信區間（credible interval） 1. 等尾可信區間：適用於對稱性分布，可以跟 95% CI 比較 2. 最高後驗密度：適用於偏斜或多峰分布 3. 分位數可信區間：適用於單峰的對稱性分布 後驗預測分布：95% 預測區間，p 值靠近 0.5 時表示模型是可靠的 假設： • log(OR), log(RR), log(HR) 都是常態分布 • 先驗 log(OR) = log(RR) = log(HR) = 0（即先驗 OR = RR = HR = 1） • log(先驗反應比率)=log(mean 1/mean2) • 大部分的 OR, RR, HR, log(先驗反應比率) 介於 0.5～2 之內  • 先驗標準平均差 SMD =平均差/標準差 sd= 0 • 大部分的 SMD 介於 -1～1（亦即二組相差一個標準差）之內 多重先驗敏感性分析： • 強先驗: gaussian(0, 0.2)，95% confidence interval（信賴區間） -0.4, 0.4, 亦即 OR, RR, HR, 先驗反應比率 0.7, 1.5 • 內定先驗: gaussian(0, 0.42)，95% CI（信賴區間） -0.84, 0.84, 亦即 OR, RR, HR, 先驗反應比率 0.4, 2.7 • 弱先驗: gaussian(0, 1)，95% CI（信賴區間） -4, 4, 亦即 OR, RR, HR, 先驗反應比率 0.14, 7.4
閱讀完整內容