I would like to perform a modified Poisson regression using JMP.

Shintaro1120

I would like to perform a modified Poisson regression using JMP.

I would also like to display confidence intervals using robust variance.

Could you please explain in detail how to do this in JMP?

Thank you very much in advance.

yusuke_ono

リプライありがとうございます。

 

ただ、修正ポアソン回帰はロバスト分散を使用することで正しい信頼区間を計算しているのですよね？ 

はい。もうちょっと詳しく「正しい」を言うと、「標本サイズが大きい場合に正しい」と言えると思います。さらに詳しく言うと、「本当はPoisson分布なんかに従っていないのに、Poisson分布を仮定したモデルで推定した結果が、もしも標本サイズが大きければ正しくなる」と言えると思います。統計学のジャーゴンを使うと、「たとえ応答変数の確率分布が何であっても、サンドイッチ分散は、ある特定のモデル（今回はPoisson分布）を仮定したモデルの推定量の分散共分散行列に対する一致推定量になっている」と言えると思います。

 

ロバスト分散とサンドイッチ分散は同じものを指しているという理解でよろしいでしょうか？ 

はい。少なくとも今回の文脈に限れば、「ロバスト分散」と「サンドイッチ分散」で私はまったく同じものを指しています。他にも、「経験分散、経験的分散（empirical variance）」という呼び方もあります。このサンドイッチ分散は、（いくつかの条件のもとで）M推定（M estimation）であれば使えます。最尤推定もM推定に含まれます。

 

■おっしゃる通り、修正Poisson回帰（modified Poisson regression）は、いったんPoisson回帰を最尤推定した後、その後、ロバスト分散によって標準誤差を求め、そのサンドイッチ標準誤差から検定や信頼区間を求めています。

 

Poisson回帰の最尤推定では、Poisson分布を仮定するので、二項分布が仮定されていません。ロバスト分散を用いることにより、たとえ応答変数の真の確率分布がPoisson分布ではないのにモデルにPoisson分布を仮定したとしても、標本サイズが大きければ妥当な評価を行います。

 

■特に疫学において、修正Poisson回帰は次のようなモチベーションで提案された手法だと思います。

 

あるイベントが生じる確率をpとした場合、よく使われるモデルにロジスティックモデルがあります。簡単のために説明変数が1つの場合にすると、モデル式は次のようなものになります。

log(p/(1-p)) = β0 + β1 x

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) = β0 + β1 x

というように、左辺をロジットではなく、割合の対数にしたものです。

（ロジット（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を超えるときがあります。そのため、二項分布の尤度をもとに最尤推定すると、途中で計算が破綻することがあります。

 

この問題を回避するために、リンク関数が対数のときと相性がいいPoisson分布をもとに最尤推定をすることが考えられます。しかし、そうすると（本来は二項分布に従っていると考えられるのに）Poisson分布で標準誤差・検定・信頼区間が計算されることになります。そこで、応答変数の確率分布の指定に頑健であるサンドイッチ分散（ロバスト分散・経験分散）を用いることが提案されました。サンドイッチ分散を指定すると、たとえ応答変数の確率分布として間違ったものを指定しても、パラメータ推定量に対するサンドイッチ標準誤差（およびサンドイッチ共分散行列）は、標本サイズが大きければ、妥当なものとなります（真の標準誤差に対する一致推定量となっています）。

 

なお、「オッズ比の解釈が難しい」というのはたぶん2点あると思います。1点は、直観的に難しいということです。「a ÷ b」という比は分かりやすいですが、「(a ÷ (1-a)) ÷(b ÷ (1-b))」だと直感的には何が何だか少しわかりません。また、リスク差やリスク比は併合可能性（collapsibility）という性質があるのに対して、オッズ比は併合可能性が成立しません。

 

上記のように、修正Poisson回帰が注目されているのは、「オッズ比だと解釈が難しい」という欠点を補うものですので、オッズ比の解釈が別に難しくなければ、従来のロジスティック回帰モデルのままでもいいのではないかと個人的には思います。（たとえば、分析の目的が「予測」や「分類」であり、回帰係数の解釈が重視されないのであれば、ロジスティック回帰のほうが修正Poisson回帰よりも使い勝手がいいと思います。）

 

■先日、述べた方法（RやPythonを呼び出す、および、JMPスクリプト言語でスクラッチでプログラミングする）は、プログラミングが必要となりますので、JMPの最も良い点（ノーコードで対話的に楽しく統計分析ができる）が消えてしまいます。しかし、JMPのユーザーでも、特にルーチン的な統計処理を行っている方でヘビーユーザーの方は、プログラミングをJMP上で行っています。

 

JMP Pro 19になればマウス操作で修正Poisson回帰ができるようになります。現状でJMP上で行いたい場合にはどうしてもプログラミングが必要となります。通常のJMPの使い方に比べると、かなりのいばらの道になります。

 

以下に、小さめの疑似データをサンプルデータとして、

 

　①RをJMPから呼び出して修正Poisson回帰を実行する方法

　②PythonをJMPから呼び出して修正Poisson回帰を実行する方法

　③JMPスクリプト言語で自分でゼロからプログラミングして修正Poisson回帰を実行する方法

 

を紹介します。まずは、ざっと見ていただいて、どれぐらい面倒なのかの大まかな目安を体感していただければと思います。以下のコードは、あくまで私個人が作成したひとつのサンプルであり、弊社が組織として提供する製品ではまったくありません（責任は私個人だけに帰します。しかし、サンプルであり、何かしらの保証があるものではまったくありません）。

 

さらに、以下の説明は、ここで取り上げている小さめの疑似的なデータだけで動くような説明になっています。

ここでは説明のためにかなり小さめのデータを用いていますが、本来、サンドイッチ分散はある程度の大規模なデータに対して妥当な方法となっています（漸近論に基づく方法となっています）。

 

①，②については、列名に日本語が入っているとうまく動かないかもしれません。また、私自身、RやPythonに詳しくないため、よりエレガントなパッケージやプログラミング方法があるかもしれません。なお、以下のPythonでのプログラミングに関しては、Copilotの提案をかなり参考にしました。

 

③については、次のことを全く考慮していません。

 a) 欠測値の処理：データに欠測値があるとエラーとなります。

 b) 計画行列のランク落ちに対するエラー処理：　線形回帰分析やPoisson回帰をはじめとする一般化線形モデルでは、説明変数間に1次従属性があると、パラメータ推定値が一意に定まりません。それに対する処理をまったく行っていません。

 c) 分離に対するエラー処理：　ロジスティック回帰やPoisson回帰には、データによっては分離（separation）という問題が生じて、パラメータ推定値が+∞もしくは-∞になる場合があります。その処理を行っていません。また、③について、可読性を重視しており、処理の効率性（特にメモリ消費量）をまったく気にせずにプログラミングしています。

 

①と③とを一つにまとめたJSLコードのファイルを添付しておきます。これを実行するには、Rがインストールされていて、JMPから呼び出せるようになっていないといけません。

 

②のPythonコードのファイルも添付しておきます。これを実行するには、上記のJSLファイルの1～7行目部分のコードでJMPデータテーブルが開いていないといけません。また、一番、最初に実行するには、2～4行目を実行しなければいけません。

 

 

■0) まず、以下の①, ②、③で用いるサンプルデータを作成します。 

0-1) ファイル(F) > 新規作成(N) > JSLスクリプト（Cntl + 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データテーブルが作成されます。

 

■①RをJMPから呼び出して修正Poisson回帰を実行する方法

 

1) Rのインストール、および、JMPからRが呼び出されるかの確認

まず、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 Init()によってJMPと接続されたRに送ります。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にて修正Poisson回帰を推定します。以下では、まずgeepackライブラリのgeeglm関数を、その次に、rqlmパッケージのrqlm関数を用いた例を示します。

 

GEE（一般化方程式推定）は、もともとは度数データなどで相関がある場合の推定方法として提案されたものですが、そこで標準誤差を計算するときにサンドイッチ分散が使われていました。各行が独立なデータでも、このサンドイッチ分散を計算する機能を用いることができます。そのためには、ID列としては各行に一意なものを用います（また、相関構造には独立を用います）。上記で述べた理由から、応答変数の確率分布にはPoisson分布を、リンク関数には他対数を指定します。

 

以下のコードでは、geeglmを用いて、X1とX2を説明変数とした修正Poisson回帰をしています（GEEの言葉で言うと、対数リンク関数で、Poisson分布を仮定したモデルで、かつ、各行を独立とみなして、サンドイッチ分散を求めた）。パラメータ推定値や信頼区間は、指数（exp）で変換しています。

結果はスクリプトウィンドウ下部のログに出力されます。

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))&nbsp;

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パッケージ

https://www.ism.ac.jp/~noma/file/software/rqlm_JPN.pdf 

 

/* 以下の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
");

 

■①PythonをJMPから呼び出して修正Poisson回帰を実行する方法

JMP18から、JMPをインストールすると、JMP内部からだけで動作するPythonが一緒にインストールされるようになりました。よって、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) パッケージや関数をimportします。

以下のコードを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を実行します。応答変数の確率分布にPoisson分布を、リンク関数にログを指定します。

また、パラメータ推定値とその信頼区間を指数変換しています。

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スクリプト言語で自分でゼロからプログラミングして修正Poisson回帰を実行する方法

 

修正Poisson回帰は、まずPoisson回帰の最尤推定で点推定値を求めます。

このときには、Newton-Raphson法という方法がよく使わています。Newton-Raphson法は反復計算なのですが、各反復ごとに対数尤度関数をパラメータで1回微分したもの(傾き、グラディエント)と、2回微分したもの（ヘッセ行列）が必要です。Poisson回帰は「一般化線形モデル（generalized linear models）」というモデル族に属しているのですが、一般化線形モデルではこの傾きとヘッセ行列の計算が簡単になります。

 

また、自分でゼロからプログラミングする場合には、名義尺度の変数は自分で指示変数（ダミー変数）を作成する必要があります。

 

細かい点の説明は省略しますが、エラー処理（欠測値をどうするか、ランク落ちをどうするか、分離をどうするか）を無視すれば、以下ぐらいのコード量で修正Poisson回帰は行えます。

 

Names Default To Here(1);
dt = Current Data Table();
/* 自分で指示変数（ダミー変数）を作成しないといけない */
x1 = Design(Column(dt, "X1")<<Get As Matrix);
//x1 = x1[0, 1::(NCol(x1)-1)];
x1 = x1[0, 2];
/* x2は連続変数なので、そのまま使える */
x2 = Column(dt, "X2")<<Get As Matrix;
/* この例では、yはすでに0,1になっているものとする */
y = Column(dt,"y")<<Get As Matrix;

/* 計画行列（切片を付ける) */
x = J(NRow(x1),1,1)||x1||x2;

/* Newton-Raphson法スタート */
b = J(NCol(x),1,0);
maxit = 1000;
check = 0;
For(i = 1, i <= maxit, i++,
eta = x*b;
pred&nbsp; =Exp(eta);
g = x`*(y - pred);
h = x` * Diag(pred) * x;
If(check == 1, Break());
new b = b + Inv(h)*g;
If(Max(Abs(new b - b))<1e-18, check=1);
b = new b;
);

/* サンドイッチ分散 */
/* Inv(h) * X` * Diag((y - pred)^2) * X * Inv(h)が、A*B*Aの形になっているので「サンドイッチ分散」*/
se = Sqrt(VecDiag(Inv(h) * X` * Diag((y - pred)^2) * X * Inv(h)));
p = 2*Normal Distribution(-Abs(b:/se));

/* 結果の表示 */
New Window("Test",
Outline Box("Result",
Table Box(
Number Col Box("exp(Est)", Exp(b)),
Number Col Box("Stderr for Est", se),
Number Col Box("Low CL exp(Est)", Exp(b - Normal Quantile(0.975)*se)),
Number Col Box("Up CL exp(Est)", Exp(b + Normal Quantile(0.975)*se)),
Number Col Box("p-value", p)
)
)
);

説明が長くなってすみません。JMP Pro 18（JMP Student Edition）上で行おうとすると、上記ぐらいのプログラミングする必要があります。なお、JMPスクリプト言語では、GUIを備えたアプリケーションも作成できます。いくつかJMPに備わっていない機能は、JMP Marketplaceというサイトで提供されています（こちらで提供されている拡張機能は、現状有姿(as-is)で提供されているものです）。残念ながら、修正Poisson回帰の機能はJMP Marketplaceでも提供されていません。もし待てるのであれば、JMP Pro 19がリリースされるまで待っていただければと思います。

 

JMPである処理ができるか、できないかについて、以下のようなレベルに分けることができるかもしれません。今回は、レベルd)に相当するものだと思います。

a) JMPやJMP Proで普通に用意されている機能　→　あっという間にできる

b) JMPやJMP Proの機能を利用する場合　→　工夫すればすぐにできる

c) JMP Marketplaceで提供されている機能　→　あくまで現状有姿だがすぐにできる

d) 上記で提供されていない機能でやや簡単なもの　→　プログラミングすればできる（特にRやPythonでできれば）

e) 上記で提供されていない機能でかなり高度なもの → かなりの専門知識がある場合にかぎりできる

 

最後に一点だけ、気になったので、冒頭の「正しさ」について補足しておきます。そもそも、回帰分析自体、一番単純な最小2乗の線形回帰も含めて、様々な使われ方があります。例えば、以下の論文では、記述、予測、制御・介入と分けています。

岩崎学（2021）

統計的因果推論の視点による重回帰分析

日本統計学会誌，50(2)，363-379

https://doi.org/10.11329/jjssj.50.363

上記のここまでの説明はあくまで修正Poisson回帰の実行方法を述べただけであり、いま置かれた状況で、回帰分析を行うことが「正しい」のかや、修正Poisson回帰を行うのが「正しい」のか、さらに、現在の状況が回帰分析が仮定している大きな前提を満たしているという意味で「正しい」のか、などの意味での「正しさ」を実現するものではまったくありません。現実の状況に対して現在の統計分析が「正しい」かどうかは、（私自身は判断することがまったくできないため）その分野の応用統計家の方にご相談することをお勧めしたいです（私自身、コンサルを受けたことはないのですが、分析そのものの代替案を提示することもあるでしょうし、また、修正Poisson回帰の実行方法も詳しく教えてくれるのではないかと思います）。

 

Modified Poisson Regression.jsl

Modified Poisson Regression.py