感謝你的回覆。
但是，修改後的泊松迴歸是透過使用穩健變異數來計算正確的置信區間，對嗎？
是的。更具體地說，「正確」是指樣本量較大時的結果是正確的。更具體地說，即使反應變數實際上並不服從泊松分佈，使用假設泊松分佈的模型估計的結果在樣本量較大時也是正確的。用統計術語來說，我認為“無論反應變數的機率分佈如何，對於假設特定模型（在本例中為泊松分佈）的模型估計量的方差-協方差矩陣，夾心方差都是一個一致的估計量。”
我是否正確地理解穩健變異數和夾心變異數指的是同一件事？
是的。至少在這種情況下，「穩健方差」和「夾心方差」的含義完全相同。另一個術語是「經驗方差」。這種夾心方差可用於M估計（在某些條件下）。最大似然估計也包含在M估計中。
■如您所說，修改後的泊松迴歸首先使用最大似然估計泊松迴歸，然後使用穩健變異數計算標準誤差，然後使用夾心標準誤差計算檢定和信賴區間。
在泊松迴歸的最大似然估計中，假設的是泊松分佈，因此不假設二項分佈。透過使用穩健方差，即使反應變數的真實機率分佈不是泊松分佈，但模型假設的是泊松分佈，只要樣本量較大，估計仍然有效。
■特別是在流行病學中，我認為提出改進的泊松回歸的方法是基於以下動機的。
如果某事件發生的機率為p，常用的模型是邏輯斯蒂模型。為簡單起見，如果只有一個解釋變量，則模型公式如下：
若 x=1 時等式左邊為 log(p1/(1-p1))，x=0 時等式左邊為 log(p0/(1-p0))，則 β1 為
β1 = log(p1/(1-p1)) - log(p0/(1-p0))
根據對數的性質，這是
β1 = log({p1/(1-p1)}/{p0/(1-p0)})
透過對兩邊進行指數變換，我們得到 exp(log(a))=a，
exp(β1) = {p1/(1-p1)}/{p0/(1-p0)}
換句話說，β1 表示對數比值比（當 x 變化 1 時），exp(β1) 表示比值比（當 x 變化 1 時）。這沒問題，但有時比值比的解釋可能比較困難。如果 p1 和 p0 接近零，
{p1/(1-p1)}/{p0/(1-p0)}≒p1/p0
因此，比值比可以用風險比（比例之比）來近似。然而，這種近似只有在 p1 和 p0 接近零時才成立。邏輯迴歸的優點之一是，當模型右側 β0+β1x 在 -∞ 和 +∞ 之間移動時，p 大於 0 且小於 1。換句話說，無論模型右側如何變化，p 幾乎總是保留機率的屬性（介於 0 和 1 之間）。 （然而，在邏輯迴歸模型中，當機率恰好為 0 或恰好為 1 時，就會出現問題。）
由於優勢比難以解釋，因此風險比對應的模型對於分析師來說很容易解釋。
換句話說，左邊是比例的對數，而不是對數。
（諸如對數（log(p/(1-p)))和對數（log(p))之類的轉換公式稱為連結函數。）
在這種情況下，log(p1) - log(p0) = β1，所以
log(p1/p0) = β1
p1/p0 = exp(β1)
換句話說，在這個具有對數連結函數的模型中，β1 是對數風險比率，exp(β1) 是風險比率。
但是，該模型存在一個問題。例如，如果右邊變成β0 + β1x = 1.0，則log(p) = 1.0，或p = exp(1.0) ≒ 2.72，雖然是機率，但p有時可能會超過1。因此，如果基於二項分佈的似然函數進行最大似然估計，計算可能會中途失敗。
為了避免這個問題，可以考慮基於泊松分佈進行最大似然估計，這種方法在連結函數為對數函數時效果很好。然而，這樣做會導致使用泊松分佈計算標準誤差、檢定和置信區間（即使它們實際上被認為遵循二項分佈）。因此，有人提出使用夾心方差（穩健方差/經驗方差），該方差在指定反應變數的機率分佈時具有穩健性。當指定夾心變異數時，即使為反應變數指定了錯誤的機率分佈，如果樣本量較大（它們是真實標準誤差的一致性估計量），參數估計量的夾心標準誤差（和夾心協方差矩陣）也將是合理的。
順便說一句，比值比難以解釋的原因可能有兩個。首先，它們直觀上很難理解。 「a ÷ b」這個比率很容易理解，但「(a ÷ (1-a)) ÷ (b ÷ (1-b))」直覺上有點難以理解。此外，風險差異和風險比都具有一種稱為可折疊性的屬性，而比值比卻不具備這種屬性。
如上所述，改進的泊松回歸正受到關注，因為它彌補了比值比難以解釋的缺點，因此我個人認為，如果比值比不是特別難以解釋，那麼堅持使用傳統的邏輯回歸模型就可以了。 （例如，如果分析的目的是“預測”或“分類”，並且回歸係數的解釋並不重要，那麼我認為邏輯回歸比改進的泊松回歸更容易使用。）
■前幾天我提到的方法（調用 R 或 Python，以及使用 JMP 腳本語言從零開始編程）都需要編程，這就抵消了 JMP 的最大優勢（無需編碼即可享受交互式統計分析）。然而，即使是 JMP 用戶，尤其是那些進行日常統計處理的重度用戶，程式設計也是在 JMP 中完成的。
使用 JMP Pro 19，您將能夠透過滑鼠操作執行改進的泊松回歸。目前，如果您想在 JMP 中執行此操作，則需要進行一些程式設計。這比使用常規 JMP 更加困難。
下面是一個小的偽數據作為範例數據：
①如何從JMP調用R並運行改進的泊松回歸
②如何透過從 JMP 呼叫 Python 來運行改進的泊松回歸
3) 如何使用 JMP 腳本語言從零開始編寫改進的泊松回歸程序
首先，我希望您快速瀏覽一下，大致了解一下它的工作量。以下程式碼只是我個人創建的範例，並非我公司作為組織提供的產品（我對此承擔全部責任，但由於它只是示例，因此不提供任何保證）。
此外，以下解釋僅適用於此處提供的少量偽數據。
儘管我們在這裡為了解釋的目的使用了相當小的數據集，但夾心方差實際上是一種適用於相對大量數據的有效方法（它是一種基於漸近理論的方法）。
關於①和②，如果列名包含日文字符，可能無法正常運作。另外，我不熟悉 R 或 Python，所以可能有更優雅的套件或程式設計方法。關於下面的 Python 編程，我大量參考了 Copilot 的建議。
關於③，完全沒有考慮以下幾點。
a) 缺失值處理：資料中缺少的值將導致錯誤。
b) 設計矩陣秩虧的誤差處理：在線性迴歸分析、泊松迴歸等廣義線性模型中，如果解釋變數之間存在線性依賴關係，則參數估計值無法唯一確定。目前沒有針對此問題進行任何處理。
c) 分離誤差處理：邏輯迴歸和泊松迴歸可能會遇到稱為分離的問題（取決於資料），這會導致參數估計變為 +∞ 或 -∞。此問題未處理。關於③，我們重視可讀性，在編寫程式碼時並未考慮處理效率（尤其是記憶體消耗）。
我附加了一個包含結合步驟 1 和 3 的 JSL 程式碼的檔案。要運行它，您必須安裝 R 並能夠從 JMP 呼叫它。
我還附加了②的 Python 程式碼檔案。要執行此程式碼，您必須開啟 JMP 資料表，如上面 JSL 檔案的第 1-7 行所示。此外，首次執行此程式碼時，您必須執行第 2-4 行。
■0) 首先，建立下面①、②、③中所使用的樣本資料。
0-1) 選擇檔案 (F) > 新建 (N) > JSL 腳本 (Ctrl + T)。
0-2) 在呼叫的腳本編輯器中右鍵單擊，然後從彈出式選單中選擇[在視窗中顯示日誌]。這將在底部顯示日誌視窗。
0-3) 將以下程式碼複製並貼上到腳本編輯器的頂部。此程式碼用於建立範例資料。
New Table("Test",
New Column("X1", Nominal, Character,Values({"A","A","A","A","A","A","B","B","B","B","B","B","B","B","B"})),
New Column("X2", Values([5,4,3,2,1,4,5,4,3,3,3,2,1,5,4])),
New Column("Y", Values([0,1,0,1,0,1,1,0,1,1,0,1,1,0,1]));
);
0-4) 在腳本編輯器中按滑鼠右鍵，從彈出式選單中選擇[執行腳本]。
0-5) 將建立具有上述資料值的JMP資料表。
■①如何從JMP調用R並運行改進的泊松回歸
1）安裝 R 並檢查是否可以從 JMP 呼叫它
首先，要從 JMP 呼叫 R，必須在安裝 JMP（JMP、JMP Pro、JMP Student Edition）的機器上安裝 R。
如果已安裝 R，或安裝 R 後，請檢查是否可以從 JMP 呼叫 R。以下介紹如何使用 Windows 版本的選單操作來執行此操作。
1-1) 將以下程式碼複製並貼上到腳本編輯器中：
R Init();
R Submit("
R.Version();
print('Hello, World!')
");
1-2) 在腳本編輯器中按滑鼠右鍵，在彈出的選單中選擇[執行腳本]。
1-3) 如果腳本編輯器下方的日誌顯示R版本和「Hello, World!」之類的字串，則表示您已成功從JMP呼叫R。
2) 接下來，將 JMP 資料表傳送到 R，R 透過上面的 R Init() 函數連接到 JMP。您可以使用 R Send 函數將 JMP 資料表傳送到 R。
2-1) 開啟資料表，執行以下程式碼（如果已經連線到R，則不需要在開始時執行R Init()）。
R Init();
dt = Current Data Table();
dt << New Column("ID", Nominal, Formula(Row()));
Wait(0);
R Send(dt);
Current Data Table() 取得目前為目前資料表的 JMP 資料表的參考。
“dt << New Column("ID", Nominal, Formula(Row()));” 會為每一行建立一個包含連續數字的欄位。如果您已經有一個包含連續數字（ID 號碼）的列，則此行並非必要。
最後的R Send(dt)將資料表傳送給R。
3）運行R程式碼。
一旦將資料傳送到 R，下一步就是在 R 中估計修改後的泊松迴歸。下面，我們首先使用 geepack 函式庫中的 geeglm 函數展示一個範例，然後使用 rqlm 套件中的 rqlm 函數。
GEE（廣義方程式估計）最初是作為頻率數據等相關數據的估計方法提出的，並使用夾心方差計算標準誤差。即使對於每行獨立的數據，也可以使用此夾心方差計算函數。為此，請使用每行唯一的 ID 列（並使用獨立性作為相關結構）。基於上述原因，指定泊松分佈作為反應變數的機率分佈，並使用多重對數作為連結函數。
以下程式碼使用 geeglm 執行修正泊松迴歸，以 X1 和 X2 作為解釋變數（以 GEE 術語來說，即對數連結函數、泊松分佈和夾心方差，每行均視為獨立變數）。參數估計值和信賴區間均進行了指數變換。
結果輸出到腳本視窗底部的日誌。
R Submit("
if (!require('geepack')) {
install.packages('geepack')
}
library(geepack)
gee.result <- geeglm(Y~X1 + X2, data=dt, id = ID, family=poisson('log'))
summary(gee.result)
est <- coef(gee.result)
se <- diag(vcov(gee.result))
z <- qnorm(0.975)
print('Risk Ratios and the Confidence Interlvals')
print(cbind(exp(est), exp(est - se*z), exp(est + se*z)))
");
下面是一個使用 rqlm 套件的範例。此範例參考了以下外部資料編寫。我認為它比上面的 geepack 範例更容易指定。
野間尚 (2025)
rqlm：用於執行改進的泊松迴歸和改進的最小平方法迴歸分析的 R 套件
/* 以下のURLを参考にしました。https://www.ism.ac.jp/~noma/file/software/rqlm_JPN.pdf */
R Submit("
if (!require('rqlm')) {
install.packages('rqlm')
}
library(rqlm)
mpr.result <- rqlm(Y ~ X1 + X2, data=dt, family=poisson, eform=TRUE)
mpr.result
");
■①如何從JMP呼叫Python並執行改進的泊松迴歸
從 JMP 18 開始，安裝 JMP 也會安裝 Python，它僅在 JMP 內運行，因此無需單獨安裝 Python。
1-0) 以下內容將使用 Python 程式碼編寫。若要在 JMP 中建立 Python 程式碼，請選擇「檔案 (F)」>「新建 (N)」>「Python 程式碼」。您也可以在 Python 腳本編輯器中右鍵單擊，然後選擇“在視窗中顯示日誌”，以在底部顯示日誌視窗。
1-1) 安裝一些必要的 Python 套件。本例中，所需的套件包括 numpy、pandas 和 statsmodels。將以下程式碼複製並貼上到 Python 腳本編輯器中並運行。
from jmputils import jpip
jpip('install --upgrade', 'pip setuptools certifi')
jpip('install', 'numpy pandas statsmodels')
這段程式碼只需要執行一次，不需要再執行一次。
1-2）導入包或函數。
將以下程式碼複製並貼上到 Python 腳本編輯器中並運行它。
import jmp
import numpy as np
import pandas as pd
import statsmodels as sm
from statsmodels.genmod.generalized_estimating_equations import GEE
from statsmodels.genmod.families import Poisson
from statsmodels.genmod.families.links import log
from statsmodels.genmod.cov_struct import Independence
1-3）
將 JMP 資料表讀入 Python 並將其轉換為 pandas 資料框。
本例中讀取目前JMP資料表（目前開啟的JMP資料表），並將資料轉換為pandas資料框。
下面的程式碼需要根據您的資料進行修改。
dt = jmp.current()
data = pd.DataFrame(
{
'Y': dt['Y'][0:n],
'X1': dt['X1'][0:n],
'X2': dt['X2'][0:n],
'ID': dt['ID'][0:n]
}
)
要使用 Python 的 GEE 函數，似乎需要自行建立指示變數（虛擬變數）。
在此範例中，X1 是二元變數。這裡，我們建立一個虛擬變量，如果 X1=="B"，則為 1，否則為 0。
data['X1B'] = (data['X1'] == 'B').astype(int)
print(data)
1-4）
運行 GEE。指定泊松分佈作為反應變數的機率分佈，並對數作為連結函數。
此外，參數估計及其信賴區間也呈現指數變換。
model = GEE.from_formula( "Y ~ X1B + X2", groups="ID",
data=data,
family=Poisson(link=log()),
cov_struct=Independence())
model = GEE(Ydata, Xdata, groups=ID, cov_struct=cov_struct, family=Poisson(link=Log()))
result = model.fit()
print(result.summary())
print(result.cov_robust)
print(np.exp(result.params))
print(np.exp(result.conf_int()))
■③如何使用JMP腳本語言從零開始自行編寫改進的泊松回歸程序
改進的泊松迴歸首先使用泊松迴歸的最大似然估計來找到點估計。
在這種情況下，通常使用牛頓-拉夫森法。牛頓-拉夫森法是一種迭代計算，每次迭代都需要對數似然函數對參數（斜率、梯度）進行一次微分，對參數（Hessian矩陣）進行兩次微分。泊松迴歸屬於「廣義線性模型」系列，這使得計算斜率和Hessian矩陣更加容易。
此外，如果您自己從頭開始編程，您將需要為名義尺度變數創建自己的指示變數（虛擬變數）。
我不會深入細節，但如果你忽略錯誤處理（如何處理缺失值、如何處理秩不足、如何處理分離），你可以使用以下程式碼量執行修改後的泊松迴歸：
Names Default To Here(1);
dt = Current Data Table();
/* 自分で指示変数（ダミー変数）を作成しないといけない */
x1 = Design(Column(dt, "X1")<
<GET as="" matrix=""></GET>
<DIV>抱歉，解釋得有點長。要在 JMP Pro 18（JMP 學生版）中實現這一點，您需要進行上述程式設計工作。 JMP 腳本語言還允許您建立具有 GUI 的應用程式。 JMP 中一些不可用的功能在 JMP Marketplace 上可用（JMP Marketplace 提供的擴充功能按原樣提供）。遺憾的是，JMP Marketplace 上不提供改進的泊松回歸功能。如果您可以等待，我建議您等到 JMP Pro 19 發布。</DIV>
<DIV>
</DIV>
<DIV>某個流程能否在JMP中執行，可以分為以下幾個等級。本例中，我認為對應的是d級。</DIV>
<DIV>a) JMP 和 JMP Pro 中常用的功能 → 可立即完成</DIV>
<DIV>
<DIV>b) 如果使用JMP或JMP Pro的功能，只要稍加巧思就能很快完成。</DIV>
</DIV>
<DIV>
<DIV>c) JMP Marketplace 提供的功能 → 保持原樣，但可立即使用</DIV>
<DIV>d）上面未提供的簡單功能→可以透過程式設計（尤其是 R 或 Python）完成</DIV>
<DIV>e) 上述未提供的高階功能 → 只有具備相當的專業知識才能實現</DIV>
<DIV>
</DIV>
<DIV>最後，還有最後一點讓我感到困惑，所以我想補充一下我在開頭提到的「正確性」。迴歸分析本身，包括最簡單的最小平方法線性迴歸，有多種用途。例如，下面這篇論文將其分為描述、預測和控制/幹預。</DIV>
<DIV>岩崎，M.（2021）</DIV>
<DIV>統計因果推論視角下的多元迴歸分析</DIV>
<DIV>日本統計學會雜誌, 50(2), 363-379</DIV>
<DIV>
<A class="bluelink-style" href="https://doi.org/10.11329/jjssj.50.363" target="_blank" rel="nofollow noopener noreferrer">https://doi.org/10.11329/jjssj.50.363</A>
</DIV>
<DIV>以上解釋僅描述如何進行修正泊松回歸，並未以任何方式確立“正確性”，即在當前情況下進行回歸分析是否“正確”，進行修正泊松回歸是否“正確”，甚至未確立當前情況是否“正確”，即它是否滿足回歸分析的主要假設。至於你目前的統計分析是否符合你的實際情況（我自己沒有能力判斷），我建議你諮詢該領域的應用統計學家（我從未諮詢過，但我相信他們或許能夠提出分析本身的替代方案，並提供關於如何進行修正泊松回歸的詳細說明）。</DIV>
<DIV>
</DIV>
</DIV>
