【深度學習】圖形模型基礎(5)：線性回歸模型第四部分：預測與貝葉斯推斷

1.引言

貝葉斯推斷超越了傳統估計方法，它包含三個關鍵步驟：結合數據和模型形成后驗分布，通過模擬傳播不確定性，以及利用先驗分布整合額外信息。本文將通過實際案例闡釋這些步驟，展示它們在預測和推斷中的挑戰和應用。

2 利用后驗模擬傳播推斷中的不確定性

在本文后續部分，我們將討論先驗分布的重要性。但首先，我們關注如何使用stan_glm函數進行模型擬合，并通過R中的模擬來表達推斷和預測中的不確定性。

**示例：**以之前文章的選舉預測模型為例，我們有以下估計結果摘要：

截距（Intercept）：中位數 46.3，MAD_SD 1.7
增長（growth）：中位數 3.0，MAD_SD 0.7
殘差標準差（sigma）：中位數 3.9，MAD_SD 0.7

這些數值來源于對參數向量 $（ a, b, σ ）$ 的模擬矩陣，代表截距、斜率和回歸的殘差標準差。我們得到的是一系列后驗模擬，而非單一的點估計，這反映了對參數的不確定性。

通過R語言提取模擬矩陣：

sims <- as.matrix(M1)

然后，我們可以直接從模擬中計算中位數和MAD_SD，以驗證摘要的準確性：

Median <- apply(sims, 2, median)
MAD_SD <- apply(sims, 2, mad)
print(cbind(Median, MAD_SD))

這里的MAD_SD是后驗模擬的縮放中位數絕對偏差，它比標準誤差或后驗標準差更穩定，因此我們有時稱之為“標準誤差”。

2.1 回歸系數的不確定性及回歸線的隱含不確定性

為了理解貝葉斯不確定性分布的概念，繪制系數的模擬抽樣圖是有幫助的。圖9.1展示了擬合的選舉預測模型的截距和斜率的后驗模擬直方圖。我們疊加了每個參數的估計值±1和2個MAD_SD，以展示數值摘要與模擬之間的關系。
在這里插入圖片描述
圖2.1 回歸模型y = a + bx + error中a和b的后驗模擬。對于每個參數，我們都展示了由stan_glm產生的模擬抽樣，以及后驗中位數±1和2倍的穩健標準誤差（mad sd）。

圖2.2 (a) 使用stan_glm擬合的選舉預測模型y = a + bx + error中，(a, b)的4000個后驗模擬結果。(b) 對應于這些(a, b)模擬結果中的100個的y = a + bx直線圖。

圖2.2a的散點圖顯示了截距和斜率的后驗分布。每個點 $（ a, b ）$ 可以映射到一條線 $y = a + b x$ ，圖2.2b展示了這些線中的100條，反映了擬合回歸的推斷不確定性。隨著樣本量的增加，系數的標準誤差減小，推斷的精度提高。

2.2 使用后驗模擬矩陣表達對參數估計或參數估計函數的不確定性

模擬總結推斷的優勢在于，我們可以直接使用這些模擬來傳播不確定性。例如，如果我們想要得到參數比率a/b的估計值和標準誤差，可以直接使用模擬矩陣：

a <- sims[, "Intercept"]
b <- sims[, "growth"]
z <- a / b
print(c(median(z), mad(z)))

這種方法允許我們直接從后驗模擬中提取信息，對參數的函數進行推斷，而不必依賴于傳統的點估計和標準誤差。

通過這種方式，貝葉斯推斷提供了一種靈活的方法來處理不確定性，并通過模擬使我們能夠對模型參數及其函數進行深入的分析和預測。

3. 預測與不確定性：點預測、后驗線性預測和后驗預測

在擬合了線性回歸模型 $\text{error}$ 后，我們可以利用它來預測一個或多個新數據點，其中預測變量為 $x_{\text{new}}$ 。我們可以進行三種不同級別的預測，這些預測反映了不確定性的不同程度：

點預測 $\hat{a} + b\hat{x}_{\text{new}}$ ：根據擬合模型，這是對于具有新 $x$ 值的新數據點 $y$ 的平均值的最佳點估計。這里使用 $\hat{a}$ 和 $\hat{b}$ 是因為點預測忽略了不確定性。
帶不確定性的線性預測器 $bx_{\text{new}}$ ，傳播了參數 $(a, b)$ 的推斷不確定性：這代表了對于具有預測變量 $x_{\text{new}}$ 的新數據點 $y$ 的預期或平均值的不確定性分布。
新觀察的預測分布 $bx_{\text{new}} + \text{error}$ ：這代表了具有預測變量 $x_{\text{new}}$ 的新觀察 $y$ 的不確定性。

例如，假設有一項研究預測血壓 $y$ 基于藥物劑量 $x$ 。對于任何給定的 $x_{\text{new}}$ ，點預測是人群平均血壓的最佳估計，條件是劑量 $x_{\text{new}}$ ；線性預測器是人群中劑量為 $x_{\text{new}}$ 的人的平均血壓模型，其不確定性對應于參數 $a$ 和 $b$ 的推斷不確定性；預測分布代表了從這個人群中隨機抽取的一個人的血壓，在模型條件下指定的 $x_{\text{new}}$ 。

隨著樣本大小接近無窮大，參數 $a$ 和 $b$ 被估計得越來越精確，線性預測器的不確定性接近零，但新觀察的預測分布的不確定性不會接近零；它接近殘差標準差 $\sigma$ 。

3.1.使用 $p re d i c t$ 進行點預測

以預測選舉經濟為例，我們預測現任政黨的得票率，條件是經濟增長 2.0%。擬合的模型 $M 1$ 出現在第 113 頁。我們首先創建一個包含假設值 $x$ 的新數據集：

new <- data.frame(growth=2.0)

然后我們可以計算點預測：

y_point_pred <- predict(M1, newdata=new)

或者，我們可以“手動”執行這個計算：

new <- 2.0
a_hat <- coef(M1)[1]
b_hat <- coef(M1)[2]
y_point_pred <- a_hat + b_hat * new

無論哪種方式，點預測都是一個單一的數字，例如 $\times 2 = 52.3$ 。

3.2.帶不確定性的線性預測器

我們可以使用 posterior_linpred 來獲取擬合回歸線的不確定性值：

y_linpred <- posterior_linpred(M1, newdata=new)

這返回一個后驗模擬向量，其均值等于上面獲得的點預測，其標準差代表擬合模型的不確定性。

或者，我們可以使用函數 posterior_epred，它返回新數據點的預期預測。對于線性回歸，期望值 $E(y|x_{\text{new}})$ 與線性預測器 $X_{\text{new}}\beta$ 相同，但對于非線性模型，這兩個量是不同的。

3.3.新觀察構建預測分布

最后，我們可以構建一個向量，表示單個選舉的預測不確定性：

y_pred <- posterior_predict(M1, newdata=new)

要“手動”使用模擬獲得這個結果，我們只需在之前的計算中加上一個誤差項：

n_sims <- nrow(sims)
sigma <- sims[,3]
y_pred <- a + b * as.numeric(new) + rnorm(n_sims, 0, sigma)

在任何情況下，我們可以通過圖形和數字來總結模擬向量：

hist(y_pred)
y_pred_median <- median(y_pred)
y_pred_mad <- mad(y_pred)
win_prob <- mean(y_pred > 50)
cat("Predicted Clinton percentage of 2-party vote: ", round(y_pred_median, 1),", with s.e. ", round(y_pred_mad, 1), "\nPr (Clinton win) = ", round(win_prob, 2), sep="")

這是結果：

Predicted Clinton percentage of 2-party vote: 52.3, with s.e. 3.9
Pr (Clinton win) = 0.72

3.4.給定一系列輸入值的預測

我們還可以使用 predict、posterior_linpred 和 posterior_predict 來生成一系列預測值，例如預測經濟增長從 -2% 到 +4% 的可能值范圍的選舉結果：

new_grid <- data.frame(growth=seq(-2.0, 4.0, 0.5))
y_point_pred_grid <- predict(M1, newdata=new_grid)
y_linpred_grid <- posterior_linpred(M1, newdata=new_grid)
y_pred_grid <- posterior_predict(M1, newdata=new_grid)

結果是長度為 13 的向量(y_point_pred_grid)和兩個n_sims $\times$ 13 矩陣(y_linpred_grid) 和y_pred_grid，對應于為增長指定的每個13 個值的預測。

3.5.傳播不確定性

在上述計算中，我們表達了在各種預設的經濟增長值條件下的選舉結果的不確定性。但是，在選舉前，增長只能估計得相當近似，而且后來政府經常會修正這些數字。因此，在應用我們的模型時，考慮到這個預測變量的不確定性是有意義的。

假設在選舉前，我們對經濟增長的最佳估計是 2.0%，但有一些不確定性，我們將這種不確定性表示為標準差為 0.3% 的正態分布。然后我們可以傳播這個預測變量的不確定性，以獲得一個更完整地表達我們不確定性的預測分布。

我們只需要在我們的 R 代碼中添加一行來模擬預測的分布：

x_new <- rnorm(n_sims, 2.0, 0.3)
y_pred <- rnorm(n_sims, a + b * x_new, sigma)

然后像以前一樣總結，得到這個結果：

Predicted Clinton percentage of 2-party vote: 52.3, with s.e. 4.1
Pr (Clinton win) = 0.71.

點預測保持不變，為兩黨投票的 52.3%，但標準差略有增加，以反映這種額外的不確定性。

3.6.線性預測器和新觀察模擬的不確定性討論

示例：我們再次使用一個模型來預測體重（以磅為單位）從身高（以英寸為單位），來自一個調查，其回答已經被清理并保存在一個名為 earnings 的數據框中：

fit_1 <- stan_glm(weight ~ height, data=earnings)
print(fit_1)

這產生了截距和高度的中位數、MAD_SD 值以及輔助參數 sigma 的中位數和 MAD_SD 值。

截距很難解釋 - 它代表了一個身高為零英寸的人的預測體重 - 所以我們進行轉換：

earnings$c_height <- earnings$height - 66
fit_2 <- stan_glm(weight ~ c_height, data=earnings)
print(fit_2)

這個回歸中的預測變量是相對于 66 英寸的身高。這種近似中心化的預測變量使得回歸系數的推斷更穩定和可解釋。

3.7.使用擬合的回歸進行推斷

在將擬合的回歸應用于新的數據點 $x_{\text{new}}$ 時，我們可以對線性預測器 $bx_{\text{new}}$ 或預測值 $y_{\text{new}} = a + bx_{\text{new}} + \epsilon$ 進行推斷。

例如，讓我們預測一個身高 70 英寸的人的體重，使得 $c_{\text{height}} = \text{height} - 66 = 4$ ：

線性預測器： $a + 4.0 b$
預測值： $\epsilon$ 。

有時我們對預測的期望或線性預測器感興趣，因為它代表這個身高的所有人在人群中的平均體重；在其他情況下，我們想要預測一個個體的體重。這取決于上下文，我們需要能夠評估每個的不確定性。

對于身高為70英寸的人，我們的模型預測其體重為153.2 + 4.0 * 4.9 = 172.8磅。如果此方程直接反映了真實世界的規律，而非基于數據的估計模型，那么我們可以將 $\hat{\sigma}$ = 29.1視為預測值標準差的一個合理估計。然而，在實際應用中，由于回歸參數（如斜率b和截距a）的估計存在不確定性，這種估計誤差的標準差往往會略高于 $\hat{\sigma}$ 。這種復雜性促使我們在一些回歸分析的文獻中見到了特定的預測標準差計算公式。

在標準的線性回歸模型 $y = a + b x + ε$ 中，線性預測器 $a + b x$ 的標準差可以通過以下公式計算：

$\hat{\sigma}_{\text{linpred}} = \sqrt[\hat{\sigma}]{\frac{1}{n} + \frac{(x_{\text{new}} - \bar{x})^2}{\sum_{i=1}^{n}(x_i - \bar{x})^2}}\quad \text{(3.1)}$

而包含隨機誤差項ε的完整預測值 $a + b x + ε$ 的標準差則為：

$\hat{\sigma}_{\text{prediction}} = \sqrt[\hat{\sigma}]{1 + \frac{1}{n} + \frac{(x_{\text{new}} - \bar{x})^2}{\sum_{i=1}^{n}(x_i - \bar{x})^2}} \quad \text{(3.2)}$

值得注意的是，盡管這些公式在簡單線性模型中非常有用，但它們并不總是適用于更復雜的模型或非線性預測場景。因此，一個更為靈活且強大的方法是利用模擬來評估預測的不確定性。

在R語言中，當我們使用stan_glm函數擬合模型后，可以通過posterior_linpred和posterior_predict函數分別獲得線性預測器和完整預測值的后驗模擬。這兩種預測方式在平均值上是一致的，即都等于點預測 $\hat{a} + \hat{b}x_{\text{new}}$ 。然而，它們在不確定性方面存在差異：隨著樣本量n的增加，線性預測器的標準差逐漸減小（與 $1/\sqrt{n}$ 成正比），而完整預測值的標準差則趨近于一個非零值 $σ$ ，這反映了即使模型參數完全已知，預測新數據點時仍會存在固有的不確定性。

3.8.后驗模擬

在使用 stan_glm 擬合模型后，我們可以使用 posterior_linpred 和 posterior_predict 獲取兩種預測的后驗模擬。

我們使用身高和體重示例進行演示，使用上述的貝葉斯回歸結果 fit_2，然后定義一個數據框，用于身高為 70 英寸的新人（因此中心化身高為 4.0 英寸）：

new <- data.frame(c_height=4.0)

我們可以計算點預測：

y_point_pred_2 <- predict(fit_2, newdata=new)

和線性預測器的模擬：

y_linpred_2 <- posterior_linpred(fit_2, newdata=new)

這產生了一個向量，表示可能的 $a + 4.0 b$ 值，變化來自系數的后驗不確定性。

我們可以為身高為 70 英寸的單個新人計算后驗預測模擬：

y_postpred_2 <- posterior_predict(fit_2, newdata=new)

產生了一個模擬向量，代表可能的 $\epsilon$ 值。

在任何情況下，我們都可以通過中位數、mad sd、直方圖等進行總結。

4. 先驗信息與貝葉斯綜合分析

傳統統計學方法根據單一數據集提供匯總和推斷。與之相對，貝葉斯方法將數據模型與先驗知識結合起來，以期獲得與這兩類信息源相一致的推斷結果。

要全面理解統計推斷，無論是貝葉斯還是傳統方法，所需的概率論知識往往超出了本書的覆蓋范圍。但是，即便不完全理解概率論的數學基礎，統計學的用戶也能取得顯著進展。本書的核心理念是直接深入探討方法本身以及理解這些方法的方式，希望那些希望深入研究的讀者可以在之后進一步學習相關理論和推導。本文末尾的參考文獻注釋為有興趣深入了解的讀者提供了一些指導。

4.1.數據和先驗信息的統一度量

我們首先介紹正態分布的貝葉斯推斷公式，這里用 $\theta$ 表示我們想要估計的連續參數。我們的目標是將先驗估計 $\theta_{\text{prior}}$ （帶有先驗標準誤差 $se_{\text{prior}}$ ）與數據估計 $\theta_{\text{data}}$ （帶有標準誤差 $se_{\text{data}}$ ）結合起來。

最終，我們得到的貝葉斯估計是 $\theta_{\text{Bayes}}$ ，其標準誤差為 $se_{\text{Bayes}}$ ，具體公式如下：

$\hat{\theta}_{\text{Bayes}} = \frac{\frac{\hat{\theta}_{\text{prior}}}{se_{\text{prior}}^2} + \frac{\hat{\theta}_{\text{data}}}{se_{\text{data}}^2}}{\frac{1}{se_{\text{prior}}^2}+\frac{1}{se_{\text{data}}^2}}\tag{4.1}$

$se_{\text{Bayes}} = \frac{1}{\sqrt{\frac{1}{se_{\text{prior}}^2} + \frac{1}{se_{\text{data}}^2}}}\tag{4.2}$

這些公式雖然初看可能有些復雜，但它們實際上揭示了在這種簡單情況下貝葉斯推斷是如何在先驗知識和數據之間做出權衡的。特別是，公式 (4.1) 實質上是先驗估計和數據估計的加權平均，每個估計的權重與其標準誤差的平方成反比。公式 (4.1) 可以理解為將 $1/se^2$ 視為精度，并認識到后驗精度等于來自先驗的精度加上來自數據的精度。

分享這些公式的原因是，這種加權平均的概念是理解貝葉斯邏輯的有用途徑。權衡的具體細節取決于先驗和數據的標準誤差：如果兩個標準誤差相等，貝葉斯估計將正好位于先驗和數據估計的中間；如果不相等，則貝葉斯估計將更接近于具有較低標準誤差的信息源。

公式 (4.2) 可能不太直觀，但可以理解為先驗和數據的逆方差相加，因此 $se_{\text{Bayes}}$ 必然低于 $se_{\text{prior}}$ 和 $se_{\text{data}}$ 中的任何一個：信息的結合提高了精確度。

4.2.貝葉斯信息聚合實例

我們通過一個簡單的例子來演示如何使用公式 (4.1) 和 (4.2) 來結合先驗信息和數據，進行后驗估計和標準誤差的計算。假設即將舉行一場選舉，基于經濟和政治情況的先前模型預測民主黨候選人將獲得兩黨票數的 52.4%，預測不確定性為 4.1%。按照上面的符號，我們有 $\theta_{\text{prior}} = 0.524$ 和 $se_{\text{prior}} = 0.041$ 。

現在，我們進行了一項包含 400 人的調查，其中 190 人表示將投票給民主黨候選人，210 人支持共和黨。如果這項調查是選民的簡單隨機樣本，并且在調查和選舉之間沒有選民改變主意，那么數據估計為 $\theta_{\text{data}} = \frac{190}{400} = 0.475$ ，標準誤差為 $se_{\text{data}} = \sqrt{\frac{0.475(1 - 0.475)}{400}} = 0.025$ 。如果考慮到可能的非抽樣誤差，更現實的分析可能會增加標準誤差，也可能改變點估計，但這里我們不討論這種情況。

圖 4.1a 展示了上述兩段中描述的先驗分布和可能性（數據中的信息）。

在這里插入圖片描述
圖4.1 (a) 似然性（代表數據估計和標準誤差的分布）和先驗分布，這個例子結合了民意調查（數據）的估計值和擬合模型（先驗）的預測值。在這種情況下，數據估計值為0.475，標準誤差為0.025，先驗估計值為0.524，先驗標準誤差為0.041。(b) 貝葉斯推斷得出后驗分布，它是似然性和先驗之間的折衷。在這個例子中，數據比先驗更具信息量，因此后驗分布更接近似然性分布而不是先驗分布。
現在，我們結合先驗和數據估計，這里的關鍵假設是，先驗和數據代表兩種不同的信息來源，這是完全合理的。更準確地說，我們假設兩個估計的不確定性在統計上是獨立的。

這里，先驗分布是基于過往選舉的預測以及當前經濟表現的衡量，預測不確定性來自于擬合模型中的非零殘差。數據來自于新選舉中潛在選民的調查，數據估計的不確定性來自于抽樣變異性，以及分析中考慮的任何非抽樣誤差。

在這個例子中，先驗標準誤差為 4.1%，數據標準誤差為 2.5%，因此數據比先驗提供了更多的信息，貝葉斯估計將更接近數據。

我們可以輕松地進行如下計算：

theta_hat_prior <- 0.524
se_prior <- 0.041
n <- 400
y <- 190
theta_hat_data <- y / n
se_data <- sqrt((y/n) * (1 - y/n) / n)
theta_hat_bayes <- (theta_hat_prior / (se_prior^2) + theta_hat_data / (se_data^2)) / (1 / (se_prior^2) + 1 / (se_data^2))
se_bayes <- sqrt(1 / (1 / (se_prior^2) + 1 / (se_data^2)))

計算結果為 $\theta_{\text{Bayes}} = 0.488$ 和 $se_{\text{Bayes}} = 0.021$ 。估計值確實位于先驗和數據估計之間，更接近數據，如圖 4.1b 所示。

現在考慮相同的情況，但數據不確定性更大。只需取相同的點估計，并將數據標準誤差 $se_{\text{data}}$ 從 0.025 增加到 0.075，貝葉斯估計值變為 0.512；現在它更接近先驗。0

4.3.分配先驗分布和執行貝葉斯計算的不同方式

統計學中可能令人困惑的一點是，相似的分析可以用不同的方法執行。例如，在第 7.3 節中我們討論了如何將兩個組的比較表達為對指示變量的回歸，相應的標準誤差可以使用公式 (4.1) 計算，或者使用最小二乘回歸的代數方法。

同樣，貝葉斯推斷，至少在具有單一先驗信息源和單一數據源的簡單例子中，可以直接使用方程 (4.1) 和 (4.2) 來執行，或者使用回歸分析的線性代數，如 stan_glm 函數實現的那樣。

另一個復雜性在于，貝葉斯推斷可以應用于任何不確定的量。在上述例子中，不確定的量是可觀察的結果，即下次選舉中現任政黨的得票率。在后續文章中，不確定的量是不同類別父母生女孩概率的差異。在第6節中，我們討論了更一般的問題，其中不確定的量是回歸系數或一組回歸系數。

在這些不同情況下，先驗分布的設置將不同。對于投票的例子，基于過去數據的回歸模型提供了單個可觀察結果的預測推斷，這可以作為與當前選舉民意調查結合時的先驗信息。在出生的例子中，我們需要有關人口差異的先驗信息。在一般回歸建模的例子中，我們必須為所有系數指定先驗信息，這在實踐中通常意味著對那些我們知之甚少或不想做出強假設的系數設置弱先驗。

5. 貝葉斯推斷實例：美貌與性別比

5.1.美貌與性別比研究案例

我們能夠利用先驗信息來提高從嘈雜研究中得出的估計的精確度。舉例來說，若干年前，一位研究者通過對3000名美國人的調查數據分析，發現了父母外貌吸引力與他們孩子性別之間的某種相關性。這項調查將成年人分為五個吸引力等級，并發現在最高吸引力等級的父母中，有56%的子女是女性，而在其他等級中這一比例為48%。觀察到的8個百分點的差異，其標準誤差（基于比例差異的標準公式）為3個百分點。

這種觀察差異超過零的兩個標準誤差，滿足了統計顯著性的標準。事實上，關于外貌更有吸引力的父母更有可能擁有女兒的說法，已經發表在科學期刊上，并受到了廣泛關注。

5.2.先驗信息的應用

然而，我們還有更多已知信息。眾所周知，人類性別比的變化范圍非常有限。例如，美國最近的一項統計報告顯示，在白人中女孩的比例為48.7%，在黑人中為49.2%。在基于出生順序、母親年齡或出生季節等因素的比較中，也發現了少于0.5個百分點的類似微小差異。鑒于外貌吸引力本身只是一個主觀測量，我們很難相信，更有吸引力和較無吸引力父母之間的差異會達到0.5%。現在，我們使用貝葉斯推斷方法來進一步分析。這里我們關注的參數 $\theta$ 是指美貌父母生女兒的概率與其他父母相比，在全美人群中的差異。通常，我們關注的是 $\theta$ 與零的比較：是否有充分證據表明 $\theta > 0$ ，這將意味著更有吸引力的父母更有可能生女兒？

5.3.先驗估計及標準誤差

我們可以將我們的科學認知表達為 $\theta$ 的先驗分布，平均值為0%，標準偏差為0.25%。先驗平均值的零表示，在看到數據之前，我們沒有理由預期外貌更有吸引力的父母生女兒的比率會有任何提高或降低。先驗標準偏差的0.25%表示，我們認為 $\theta$ 的真實值超過0.5%或低于-0.5%的可能性非常低。

為了方便，我們在百分比尺度上表達我們的估計和不確定性，避免使用如0.0025這樣的表達式，以免在計算中丟失小數點。

按照百分比尺度，這項調查提供了 $\theta_{\text{data}} = 8\%$ 的估計值，以及 $se_{\text{data}} = 3\%$ 的標準誤差。我們現在可以看到，與數據相比，先驗信息具有更高的信息量：數據的標準誤差是先驗不確定性的10倍以上。

5.4.貝葉斯估計

使用公式 (4.1) 和 (4.2) 結合先驗和數據，我們得到 $\theta_{\text{Bayes}} = 0.6\%$ 和 $se_{\text{Bayes}} = 0.25\%$ ：根據貝葉斯方法，美貌和非美貌父母生女孩比例的差異估計為0.06%（即不到0.1個百分點），標準誤差為0.25%。這一估計值非常小，遠低于標準誤差，表明這項調查的數據實際上并沒有提供有關性別比變化的有價值信息。

怎么會這樣？3000的樣本量聽起來相當不錯，但為何它在這里幾乎未提供任何信息？原因是我們在這個研究中需要比一般調查設置更高的精確度。以下是一些快速計算來說明這一點：

來自大小為3000的簡單隨機樣本的估計比例的標準誤差 $\sqrt\frac{p(1-p)}{3000}$ ，如果 $p$ 接近0.5，大約是 $\sqrt{\frac{0.5 \times 0.5}{3000}} = 0.009$ 。3000人的調查讓我們能夠將比例估計在大約1%的精度內。
現在考慮比較兩個各為1500人的群體的比例差異。這些比例差異的標準誤差 $\sqrt{\frac{p_1(1-p_1)}{1500} + \frac{p_2(1-p_2)}{1500}}$ ，如果 $p_1$ 和 $p_2$ 接近0.5，大約是 $\sqrt{\frac{2 \times 0.5 \times 0.5}{1500}} = 0.018$ 。因此，在最佳條件下，我們可以將比例差異估計在大約2個百分點內。
來自調查的實際標準誤差稍高，為3.3個百分點，因為兩個比較組的人數并不相等；大約只有10%的受訪者被標記為“非常有魅力”。例如， $\sqrt{\frac{0.5 \times 0.5}{300} + \frac{0.5 \times 0.5}{2700}} = 0.03$ 。
對比例差異的估計達到2或3個百分點的精度已經相當不錯，對于許多目的來說已經足夠。但是，對于研究性別比的差異，我們需要更高的精確度，更接近0.1個百分點的尺度。

圖4.1 (a) ：線性回歸y = a + bx + error中參數a和b的似然函數，用于描述經濟增長xi上的選舉結果yi。(b) ：似然函數在眾數處的模式，橢圓總結了眾數處對數似然函數的逆二階導數矩陣。? 從對應于平坦先驗的后驗分布中隨機抽取的4000個(a, b)模擬值。

6. 回歸中的均勻先驗、弱信息先驗和信息先驗

在前面的章節中，我們通過模擬執行了貝葉斯分析，默認使用了均勻先驗來捕捉估計和預測中的不確定性。在那些例子中，現有數據相對于其他容易獲得的信息來說較為強大，因此我們并不過分擔憂先驗分布的選擇。然而，當數據提供的信息較少時，我們就需要更慎重地考慮分析中應包含哪些先驗信息。在本節中，我們將更正式地討論這個話題，并探討如何使用stan_glm來處理均勻先驗、弱信息先驗和信息先驗。

6.1.均勻先驗分布

均勻先驗，也稱為非信息性先驗，意味著我們對參數的先驗知識是中立的。在這種情況下，后驗分布僅僅是似然函數乘以一個常數，最大似然估計也成為了后驗分布的眾數。stan_glm可以通過設置prior參數為NULL來使用均勻先驗：

M3 <- stan_glm(vote ~ growth, data=hibbs,prior_intercept=NULL, prior=NULL, prior_aux=NULL)sims <- as.data.frame(M3)a <- sims[,1]b <- sims[,2]plot(a, b)

為了理解似然與先驗的關系，我們可以提取系數的后驗模擬并將其可視化。

6.2.默認先驗分布

默認情況下，stan_glm使用的是弱信息先驗，這種先驗足以在數據較少時保持推斷結果在合理范圍內。對于一個具有形式 $y = a + b 1 x 1 + b 2 x 2 + ... + b K x K + error$ 的模型，每個系數 $b_k$ 都有一個均值為0、標準差為 $2.5\times \frac{sd(y)} {sd(x_k)}$ 的正態先驗分布。截距 $a$ 沒有直接賦予先驗分布，而是將先驗賦予了在數據中預測變量 $x$ 的平均值處的 $y$ 的預期值，即：
$\bar{x}_1 + b2 \bar{x}_2 + ... + bK \bar{x}_K$
這個居中的截距被賦予了一個均值為 $\mu_y$ 、標準差為 $2.5\times{sd(y)}$ 的正態先驗分布。殘差標準偏差 $σ$ 被賦予了一個速率參數為 $1/ s d (y)$ 的指數先驗分布。

在stan_glm中使用先驗分布時，我們需要注意，默認的先驗設置實際上并不是純粹的貝葉斯先驗，因為它們是依賴于數據來確定其具體形態的。按照貝葉斯理論，先驗分布應當在我們觀察到數據之前就已經確定。技術上，我們可以把stan_glm的默認先驗看作是對理想貝葉斯先驗的一種近似處理，這種處理方式本意是根據變量( x )和( y )在總體中的均值和標準差來進行縮放，而在實際操作中，我們通常使用樣本的均值和標準差作為這些參數的估計。

若要展示通過stan_glm擬合的模型中所使用的先驗分布，我們可以利用prior_summary()函數來達成這一目的。

在R語言環境中，通過輸入vignette("priors", package="rstanarm")可以調閱有關rstanarm包中先驗工作機制的更詳盡的說明。

為了深入理解默認先驗的含義，我們應該回憶一下，回歸系數本質上表示了在預測變量 $x$ 上相差一個單位的兩個人在結果 $y$ 上的預期差異。如果我們已經對 $x$ 和 $y$ 進行了標準化處理，那么這個系數就代表了 $x$ 每變化一個標準差， $y$ 中標準差的預期變化量。通常，我們預期這種差異的絕對值會小于1。因此，一個均值為0、標準差為2.5的正態先驗分布能夠幫助我們將不穩定的系數估計值向這個范圍內聚合。選擇2.5作為縮放因子可能看起來有些隨意，但這是為了在估計過程中引入一定的穩定性，同時確保當數據本身提供的信息量達到中等水平時，對系數估計的影響降到最低。

默認先驗分布可以隱式指定：

M1 <- stan_glm(vote ~ growth, data=hibbs)

或者顯式指定：

sd_x <- sd(hibbs$growth)
sd_y <- sd(hibbs$vote)
mean_y <- mean(hibbs$vote)
M1a <- stan_glm(vote ~ growth, data=hibbs, prior=normal(0, 2.5*sd_y/sd_x),prior_intercept=normal(mean_y, 2.5*sd_y), prior_aux=exponential(1/sd_y))

這兩個模型本質上是相同的，唯一的區別在于對于M1a我們直接指定了先驗分布。

6.3.弱信息先驗的設定

默認的先驗設置旨在確保推斷結果的穩健性。然而，在許多情況下，通過引入針對具體問題量身定制的先驗信息，我們可以進一步提高分析的準確性。以選舉預測為例，即使不依賴現有數據，我們也能對模型中的系數提出一些假設，該模型表達為： $\text{vote} = a + b \times \text{growth} + \text{error}$ 。

先看截距項。在stan_glm函數中，prior_intercept參數默認根據預測變量的平均值來設定，即 $\bar{x}$ 。在此背景下，它反映了當經濟增長處于歷史平均水平時，現任政黨所獲得的兩黨票數的平均比例。即使不參考現有數據，我們也知道投票百分比 $y$ 應在0到100之間，很可能接近50%。因此，我們可能會選擇為 $\bar{x}$ 設定一個均值為50、標準差為10的正態先驗分布。

接下來考慮經濟增長對現任政黨票數百分比的影響系數 $b$ 。直覺上，我們可能期望 $b$ 為正，但僅憑先驗知識過于確定是不明智的；我們需要保留對數據可能揭示的意外情況的開放性。那么， $b$ 的大小可能在什么范圍內呢？過去一個世紀的經濟增長通常在0到4%之間變動，而經濟增長變量 $x$ 以百分比形式編碼。1%的經濟增長差異是否可能導致預測的票數差異達到10%？這是有可能的，但很難想象 $b$ 的值會更大。因此，我們為 $b$ 指定一個均值為5、標準差為5的正態先驗分布，這表明 $b$ 很可能是正的，并且很可能小于10。

這些先驗選擇并非唯一合理的，但它們展示了如何利用對問題結構的了解來設定一個弱信息先驗，為預期結果提供一個大致范圍。

現在，我們可以整合這些信息來擬合模型：

M4 <- stan_glm(vote ~ growth, data=hibbs,prior=normal(5, 5), prior_intercept=normal(50, 10))

擬合模型后，我們打印出結果：

Median MAD_SD
(Intercept) 46.2 1.7
growth 3.1 0.7
Auxiliary parameter(s):
Median MAD_SD
sigma 3.9 0.7

這些推斷結果與之前使用非信息性先驗和默認先驗得到的結果基本一致，這在當前例子中是合理的，因為我們的先驗與數據相比包含的信息量較少（可以從 $b$ 的后驗標準誤差0.7遠小于先驗標準差5這一事實中看出）。

6.4.信息先驗的應用示例：美貌與性別比

在選舉預測的案例中，加入信息先驗雖無傷大雅，但并非必要。然而，有些情況下，由于數據的噪聲較大，先驗信息在維持合理的估計結果方面顯得尤為重要。

作為示例，我們再次審視第9.4節中提到的美貌與性別比的關系，并將其放入回歸分析的框架中考察。

相關數據來源于對美國青少年的吸引力進行評估的一項研究，這些評估是通過面對面調查，由訪問員根據五分制量表進行的。多年后，這些受訪者中的許多人已育有子女，圖9.5a展示了不同吸引力等級的父母所生女孩的百分比。原先對這些數據的分析是將“極具吸引力”的父母的孩子性別比與所有其他父母的孩子性別比進行比較（這也是第9.4節中討論的8個百分點估計值的來源），但在這里，我們將美貌視為一個連續的預測變量，使用lm進行最小二乘回歸分析，得到以下結果：

coef.est coef.se
(Intercept) 49.4 1.9
x 1.5 1.4

圖6.1b展示了最小二乘回歸線。從斜率的標準誤差可以看出，斜率非常不確定，這表明僅憑數據本身，我們對x和y之間關系的了解非常有限。

在這里插入圖片描述
圖6.1 (a) 一項調查顯示，來自五個不同吸引力類別的父母所生的女孩出生百分比數據；(b) 帶有擬合回歸線的數據。圖9.6展示了擬合中的不同不確定性表達方式。

接下來，我們使用默認先驗進行貝葉斯推斷：

fit_default <- stan_glm(y ~ x, data=sexratio)
print(fit_default)

得到的結果如下：

Median MAD_SD
(Intercept) 49.4 2.0
x 1.4 1.4
Auxiliary parameter(s):
Median MAD_SD
sigma 4.6 1.7

這與最小二乘估計幾乎相同，這是合理的，因為默認先驗設計得較弱。

然后，我們考慮一個信息先驗。如第9.4節所討論的，女孩出生的百分比通常非常穩定，在48.5%到49%之間。我們可以利用這些信息為回歸系數設定弱信息先驗：

對于截距 $a$ ，stan_glm中的prior_intercept是針對數據中預測變量的平均值設定的；因此我們需要一個針對平均美貌水平父母的女孩百分比的先驗。對于這個未知量，我們選擇一個以48.8為中心，標準差為0.5的正態分布，表明在平均吸引力水平下y的預期值大致在48.8 ± 0.5的范圍內。
對于斜率 $b$ ，我們選擇一個以0為中心，標準差為0.2的正態分布，表明我們沒有先驗期望女孩的出生與父母的美貌呈正相關或負相關，我們期望該系數大致在-0.2到0.2之間。數據中的預測變量x的范圍為4，因此這相當于說，我們期望在人群中，最高和最低吸引力類別之間女孩出生百分比的差異不超過0.8。

然后我們擬合一個信息先驗回歸：

fit_post <- stan_glm(y ~ x, data=sexratio,prior=normal(0, 0.2), prior_intercept=normal(48.8, 0.5))
print(fit_post)

擬合模型的摘要如下：

Median MAD_SD
(Intercept) 48.8 0.5
x 0.0 0.2
Auxiliary parameter(s):
Median MAD_SD
sigma 4.3 1.3

x的估計系數并不完全是零；如果我們感興趣，我們可以將結果顯示到更多的有效數字，但無論如何，顯然任何估計都被不確定性所淹沒。在這種情況下，數據提供的信息如此之少，以至于后驗基本上與先驗相同。

圖6.2展示了使用默認先驗（上圖）和信息先驗（下圖）擬合的模型的 $(a, b)$ 的后驗模擬以及相應的回歸線 $y = a + b x$ 。在這種情況下，先驗比數據提供了更多的信息。
在這里插入圖片描述

圖6.2 系數和回歸線的100次后驗模擬，針對圖6.1中的美丑與性別比例數據，在給定兩種不同的先驗下擬合。頂部圖展示了與數據和默認先驗一致的參數值及其對應的回歸線的廣泛范圍。底部圖展示了使用信息先驗擬合模型的結果，這極大地縮小了可能性的范圍。在給定這個先驗的情況下，數據幾乎沒有為a和b提供任何證據。

7.文獻注釋

想要深入了解貝葉斯方法論、相關技術和應用案例的讀者，我們有以下推薦。首先，可以閱讀 Gelman 等人在 2013 年出版的書籍，以及 McElreath 在 2020 年的著作。在當代，貝葉斯推斷大多通過模擬手段而非解析公式來執行，盡管如此，我們認為像公式（4.1）和公式（4.2）這樣的解析表達式對于理解先驗信息與局部數據之間的相對權重具有啟發性作用。對于那些對弱信息先驗感興趣的讀者，我們建議閱讀 Gelman、Simpson 和 Betancourt 在 2017 年的研究成果。

此外，關于美貌與性別比的案例分析，Gelman 和 Weakliem 在 2009 年的論文中進行了深入探討。這篇論文不僅展示了在統計分析中如何運用貝葉斯方法和非貝葉斯方法來利用先驗信息，而且對于理解貝葉斯分析在實際問題中的應用提供了寶貴的視角。