感谢你的回复。
但是，修改后的泊松回归通过使用稳健方差来计算正确的置信区间，对吗？
是的。更具体地说，“正确”是指样本量较大时的结果是正确的。更具体地说，即使响应变量实际上并不服从泊松分布，使用假设泊松分布的模型估计的结果在样本量较大时也是正确的。用统计学术语来说，我认为可以这样说：“无论响应变量的概率分布如何，对于假设特定模型（在本例中为泊松分布）的模型估计量的方差-协方差矩阵，夹心方差都是一个一致的估计量。”
我是否正确地理解稳健方差和夹心方差指的是同一件事？
是的。至少在这种情况下，“稳健方差”和“夹心方差”的含义完全相同。另一个术语是“经验方差”。这种夹心方差可用于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>
