cancel
Showing results for 
Show  only  | Search instead for 
Did you mean: 
  • JMP 19 is here! See the new features at jmp.com/new.
  • Register to attend Discovery Summit 2025 Online: Early Users Edition, Sept. 24-25.
Choose Language Hide Translation Bar
Shintaro1120
Level II

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

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.

1 ACCEPTED SOLUTION

Accepted Solutions

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

リプライありがとうございます。
 
ただ、修正ポアソン回帰はロバスト分散を使用することで正しい信頼区間を計算しているのですよね? 
はい。もうちょっと詳しく「正しい」を言うと、「標本サイズが大きい場合に正しい」と言えると思います。さらに詳しく言うと、「本当は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)) 

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
");

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

View solution in original post

11 REPLIES 11
Victor_G
Super User

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

Hi @Shintaro1120,

Welcome in the Community !
Poisson regression is available through Generalized Regression models in JMP Pro. 
You can check the available distributions for these models.

Could you please add more details (and if possible an anonimized or toy dataset) about what you're trying to achieve if you want a more detailed guidance/help from JMP users ?

Victor GUILLER

"It is not unusual for a well-designed experiment to analyze itself" (Box, Hunter and Hunter)

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

Hello, I am Yusuke Ono, Senior Tester at JMP Japan Group in SAS Institute Japan.

 

You can fit modified Poisson regression (and modified OLS regression) by Generalized Linear Mixed Model platform on JMP Pro 19, the next new major JMP release.

 

When you use this feature, the Y variable should be 0 and 1 (even if Y is originally binary categorical variable) in your data. You can convert the binary data to 0-1 indicator variable by Cols-Utilities-Indicator Columns. After preparing the 0-1 data, select Analyze - Fit Model. Select Generalized Linear Mixed Model for the Personality and select Poisson for Distribution (If you want to fit modified OLS, select Normal for Distribution). Set the 0-1 Y variable and Model Effects. Click Run button. Select Model Report - Empirical Standard Errors (Standard Errors by Sandwich Variance in Japanese) from the red triangle menu at the Poisson outline in the report window.  Now, all standard errors, tests, and confidence intervals are calculated by sandwich variances (a.k.a. empirical variances or robust standard variances). Note that the point estimates are not changed by this menu. Only standard errors, tests, and confidene intervals are changed. If you want Exp(b) where b is the point estimates or confidence limits, to calculated point estimates and confidence intervals for risk ratios, you need to output the estimation table to a JMP data table, and convert them to Exp(b) by yourself.

 

This "Sandwich variance" feature in JMP Pro 19's Generalized Linear Mixed Models (and Mixed Models) is really one of new features for JMP users in epidemiology and medical area. By this feature, you can also calculate standard errors for the Inverse Probability Weighting (IPW) (with assuming the weights are known).

 

On JMP 18 or before, there is no feature for fitting modified Poisson regression. You can fit Log-Binomial models (where the Distribution is binomial and the link function is Log function) on the Generalized Linear Models platform. But the fitting may fail and the predict probability values exceed zero or one in the Log-Binomial models. You can fit Genearlized Linear Models by selecting Analyze - Fit Model, and select Genearlized Linear Models for the Personality. The Generalized Linear Models is supported by standard JMP 18 or before. (The Genearalized Linear Mixed Models is supported only by JMP Pro.)

Yusuke Ono (Senior Tester at JMP Japan)

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

CORRECTION:

 

If you want Exp(b) where b is the point estimates or confidence limits, to calculated point estimates and confidence intervals for risk ratios, you need to output the estimation table to a JMP data table, and convert them to Exp(b) by yourself.

You can calculate Exp(b) in modified Poisson regression by Multiple Comparions. You can select this menu from the red triangle menu for the Poisson outline. Sorry for my misunderstanding.

Yusuke Ono (Senior Tester at JMP Japan)
Shintaro1120
Level II

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

Thank you for your reply.
Thank you also for the detailed explanation.
However, I am stuck halfway through.
I had been planning to use Generalized Linear Models to perform modified Poisson regression, but based on your advice, I will use Generalized Linear Mixed Model for the Personality to perform modified Poisson regression.
“Select Model Report - Empirical Standard Errors (Standard Errors by Sandwich Variance in Japanese) from the red triangle menu at the Poisson outline in the report window. Now, all standard errors, tests, and confidence intervals are calculated by sandwich variances (a.k.a. empirical variances or robust standard variances).” I don’t quite understand this. Could you please explain this part and your correction in more detail, even in Japanese?
I have taken a screenshot to illustrate the current situation and have attached it. I would appreciate it if you could review it.

Thank you very much for your assistance.

Shintaro1120
Level II

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

By the way, I am using JMP Pro 18. Could you please help me perform the modified Poisson regression?

Thank you very much for your assistance.
Best regards,

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

I am sorry for the confusion. The feature I explained is a new feature for JMP Pro 19. You cannot use the new feature on JMP Pro 18.

 

If you select Poisson distribution and Log link function in the Generalized Linear Models platform, the standard Poisson regression is fit. It assumes Poisson distribution, not binomial distribution, for the distribution of Y variable. The point estimates in the standard Poisson regression models are the same as the ones in the modified Poisson regression models, but standard errors are different. In modified Poisson regression, the standard errors are calculated by "sandwich variances", so the standard errors are valid also when the true distribution is binomial (if the sample size is large).

 

Even in JMP 18 or before, you can fit log-binomial model. If the response distribution is selected as Binomial and the link function is selected as Log in the Generalized Linear Models platform, then log-binomial model is estimated. Note that the optimization might fail (This is one of some reasons why the modified Poisson regression is preferred in practice). And, you need to be careful that in JMP's Generalized Linear Models platform, "effect-coding" is used for nominal X variables (like the Standard Least Squares and Nominal Logistic platform).

 

If you really want to do the modified Poisson regression on JMP Pro 18 or before, you need to program the algorithm from scratch with JMP Scripting Language(JSL), or program the code to invoke a routine to fit modified Poisson regression models in R or Python from JMP with JMP Scripting Language(JSL).

 

If you have more questions, I would like you to reply again.

 

 

混乱させてしまいすみません。先月に説明した機能はJMP Pro 19のものです。こちらの機能はJMP Pro 18以前ではまだ実装されていません。

 

もし「一般化線形モデル」プラットフォームにおいて、Poisson分布を応答変数の分布に、対数をリンク関数に指定すると、通常のPoisson回帰が実行されます。この分析では、二項分布ではなく、Poisson分布が仮定されます。通常のPoisson回帰の回帰係数に対する点推定値は、修正Poisson回帰のものと同じです。しかし、標準誤差は両者で異なります。修正Poisson回帰では、標準誤差は「サンドイッチ分散」によって計算されます。この「サンドイッチ分散」を用いることにより、真の分布が二項分布である場合でも、(大標本の場合には)標準誤差が妥当なものとなります。

 

JMP 18以前でも、対数-二項モデルならばあてはめることができます。「一般化線形モデル」プラットフォームにて、分布を二項分布に、リンク関数を対数に指定すると、この対数-二項モデルがあてはめられます。ただし、対数尤度を最大化するための最適化計算が失敗する場合もあります(対数-二項分布では最適化計算が失敗することが、修正Poisson回帰のほうが好まれる理由の1つです)。また、「一般化線形モデル」では、(「標準最小2乗」や「名義ロジスティック」と同様)名義尺度の説明変数に対して「効果コーディング」というコーディング方法を使っています。

 

どうしても、JMP Pro 18で修正Poisson回帰を行いたい場合には、(1) 自分自身でアルゴリズムをプログラミングするか、(2) RもしくはPythonで該当の機能を呼び出すプログラムを作成する、といったJMPスクリプト言語(JSL)によるプログラミングが必要となってきます。JMP Pro 19からは、先日申し上げた方法で簡単に実行できるようになります。

 

何かご不明な点などがありましたら、遠慮なくご返信ください。

Yusuke Ono (Senior Tester at JMP Japan)
Shintaro1120
Level II

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

I am a beginner when it comes to regression analysis, so I apologize for asking such a basic question.


I understand that the same value is derived for point estimates in Poisson regression and modified Poisson regression.

 

However, modified Poisson regression calculates the correct confidence interval by using robust variance, right?

Is it correct to understand that robust variance and sandwich variance refer to the same thing?

 

Thank you for suggesting a solution using JMP Pro 18. I am at a loss because I do not know how to use R or Python, but I would like to implement modified Poisson regression somehow.

 

回帰分析に関しまして、私は初学者なので初歩的な質問をしてしまい申し訳ございません。


点推定値に関しては、ポアソン回帰と修正ポアソン回帰とで同じ値が導かれるとのこと、承知いたしました。

 

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

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

 

JMP Pro 18を使用する場合の解決方法を提示してくださりありがとうございます。RやPythonの使い方がわからないので困り果てていますが、なんとかして修正ポアソン回帰を実装したいと思います。

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

リプライありがとうございます。
 
ただ、修正ポアソン回帰はロバスト分散を使用することで正しい信頼区間を計算しているのですよね? 
はい。もうちょっと詳しく「正しい」を言うと、「標本サイズが大きい場合に正しい」と言えると思います。さらに詳しく言うと、「本当は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)) 

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
");

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

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

いつも大変お世話になっております。

 

何度もすみません。

私は現在、MPH(Master of Public Health)課程の2年生です。

繰り返し申し上げています通り、私はJMPを使用して修正ポアソン回帰を実装したいです。

 

私のStudent editonのJMP Pro 18の有効期限は、28Feb2026です。

JMP Pro 19を購入しようと思いましたが、高額なため大学院生の私では手が出せませんでした。

 

その代わりに無料トライアルのJMP 19をインストールしましたが、修正ポアソン回帰は実装できないように思います。

やはり 様がおっしゃるとおり、

JMPのみで修正ポアソン回帰を実装するためにはJMP Pro 19が必要なのでしょうか?

 

可能でありましたら、解決策をご教示いただけますと幸いです。

 

ご指導ご鞭撻のほど何卒よろしくお願い申し上げます。

Recommended Articles