Thank you for your reply.
However, modified Poisson regression calculates the correct confidence intervals by using robust variance, right?
Yes. To be more specific about "correct," I think I would say it's correct when the sample size is large. Even more specific, I think I would say that the results estimated using a model that assumes a Poisson distribution, even though the response variable doesn't actually follow a Poisson distribution, will be correct if the sample size is large. Using statistical jargon, I think I could say that "no matter what the probability distribution of the response variable is, the sandwich variance is a consistent estimator for the variance-covariance matrix of the estimator of a model that assumes a specific model (in this case, a Poisson distribution)."
Am I correct in understanding that robust variance and sandwich variance refer to the same thing?
Yes. At least in this context, I mean exactly the same thing by "robust variance" and "sandwich variance." Another term is "empirical variance." This sandwich variance can be used for M estimation (under certain conditions). Maximum likelihood estimation is also included in M estimation.
■As you said, modified Poisson regression first performs maximum likelihood estimation of Poisson regression, then calculates the standard error using robust variance, and then calculates the test and confidence interval from the sandwiched standard error.
In maximum likelihood estimation of Poisson regression, a Poisson distribution is assumed, and therefore a binomial distribution is not assumed. By using robust variance, even if the true probability distribution of the response variable is not a Poisson distribution but the model assumes a Poisson distribution, the estimation is valid as long as the sample size is large.
■In epidemiology in particular, I believe that modified Poisson regression is a method proposed with the following motivations.
If the probability of an event occurring is p, a commonly used model is the logistic model. For simplicity, if there is only one explanatory variable, the model formula is as follows:
If the left side of the equation when x=1 is log(p1/(1-p1)), and the left side when x=0 is log(p0/(1-p0)), then β1 is
β1 = log(p1/(1-p1)) - log(p0/(1-p0))
By the properties of logarithms, this is
β1 = log({p1/(1-p1)}/{p0/(1-p0)})
By exponentially transforming both sides, we get exp(log(a))=a,
exp(β1) = {p1/(1-p1)}/{p0/(1-p0)}
In other words, β1 represents the log odds ratio (when x changes by 1), and exp(β1) represents the odds ratio (when x changes by 1). This is fine, but sometimes it can be difficult to interpret the odds ratio. If p1 and p0 are close to zero,
{p1/(1-p1)}/{p0/(1-p0)}≒p1/p0
Thus, the odds ratio can be approximated by the risk ratio (ratio of proportions). However, this approximation only holds true if p1 and p0 are close to zero. One of the advantages of logistic regression is that p is greater than 0 and less than 1 when the right-hand side of the model, β0+β1x, moves between -∞ and +∞. In other words, no matter what the right-hand side of the model looks like, p will almost always retain the properties of a probability (between 0 and 1). (However, in a logistic regression model, problems arise when the probability is exactly 0 or exactly 1.)
Since odds ratios are difficult to interpret, the model corresponding to risk ratios is easy for analysts to interpret.
In other words, the left side is the logarithm of the proportion rather than the logit.
(Transformation formulas such as logit (log(p/(1-p))) and logarithm (log(p)) are called link functions.)
In this case, log(p1) - log(p0) = β1, so
log(p1/p0) = β1
p1/p0 = exp(β1)
In other words, in this model with a log link function, β1 is the log risk ratio and exp(β1) is the risk ratio.
However, this model has a problem. If the right-hand side becomes, for example, β0 + β1x = 1.0, then log(p) = 1.0, or p = exp(1.0) ≒ 2.72, and even though it is a probability, p can sometimes exceed 1. As a result, if you perform maximum likelihood estimation based on the likelihood of a binomial distribution, the calculation may fail midway.
To avoid this problem, one could consider performing maximum likelihood estimation based on the Poisson distribution, which works well when the link function is logarithmic. However, doing so would result in standard errors, tests, and confidence intervals being calculated using the Poisson distribution (even though they are actually thought to follow a binomial distribution). Therefore, it was proposed to use sandwich variance (robust variance/empirical variance), which is robust when specifying the probability distribution of the response variable. When sandwich variance is specified, even if the wrong probability distribution is specified for the response variable, the sandwich standard errors (and sandwich covariance matrix) for the parameter estimators will be reasonable if the sample size is large (they are consistent estimators of the true standard errors).
Incidentally, there are probably two reasons why odds ratios are difficult to interpret. The first is that they are intuitively difficult. The ratio "a ÷ b" is easy to understand, but "(a ÷ (1-a)) ÷ (b ÷ (1-b))" is a little hard to understand intuitively. Also, while risk differences and risk ratios have a property called collapsibility, odds ratios do not have this property.
As mentioned above, modified Poisson regression is gaining attention because it compensates for the drawback that odds ratios are difficult to interpret, so I personally think that if odds ratios are not particularly difficult to interpret, it would be fine to stick with the traditional logistic regression model. (For example, if the purpose of the analysis is "prediction" or "classification" and interpretation of the regression coefficients is not important, I think logistic regression is easier to use than modified Poisson regression.)
■The methods I mentioned the other day (calling R or Python, and programming from scratch using the JMP scripting language) require programming, which eliminates the best thing about JMP (enjoying interactive statistical analysis without coding). However, even among JMP users, especially heavy users who perform routine statistical processing, programming is done in JMP.
With JMP Pro 19, you'll be able to perform modified Poisson regression with mouse operations. Currently, if you want to do this in JMP, you'll need to do some programming. This is a much more difficult path than using regular JMP.
Below is a small pseudo-data as sample data:
①How to call R from JMP and run a modified Poisson regression
②How to run a modified Poisson regression by calling Python from JMP
3) How to program modified Poisson regression yourself from scratch using the JMP scripting language
First, I'd like you to take a quick look at it to get a rough idea of how much work it is. The code below is merely a sample created by me personally, and is in no way a product provided by our company as an organization (I take sole responsibility for it, but as it is a sample, it does not come with any guarantees whatsoever).
Furthermore, the following explanation is designed to work only with the small, pseudo-data presented here.
Although we are using fairly small data sets here for the purposes of explanation, sandwich variance is actually a valid method for relatively large amounts of data (it is a method based on asymptotic theory).
Regarding ① and ②, if the column names contain Japanese characters, they may not work properly. Also, I'm not familiar with R or Python, so there may be more elegant packages or programming methods. Regarding the Python programming below, I heavily referenced Copilot's suggestions.
Regarding ③, the following points are not taken into consideration at all.
a) Missing Value Handling: Missing values in the data will result in an error.
b) Error handling for rank deficiency in the design matrix: In generalized linear models such as linear regression analysis and Poisson regression, if there is linear dependence between explanatory variables, parameter estimates cannot be uniquely determined. No processing is done to address this issue.
c) Error handling for separation: Logistic regression and Poisson regression can encounter a problem called separation depending on the data, which can cause parameter estimates to become +∞ or -∞. This issue is not handled. Regarding ③, we place importance on readability and have programmed without any concern for processing efficiency (especially memory consumption).
I have attached a file with the JSL code that combines steps 1 and 3. To run this, you must have R installed and be able to call it from JMP.
I've also attached the Python code file for step ②. To run this, you must have a JMP data table open, as per lines 1-7 of the JSL file above. Also, to run it for the first time, you must run lines 2-4.
■0) First, create the sample data to be used in ①, ②, and ③ below.
0-1) Select File (F) > New (N) > JSL Script (Ctrl + T).
0-2) Right-click within the called script editor and select [Show log in window] from the popup menu. This will display the log window at the bottom.
0-3) Copy and paste the following code into the top of the script editor. This code creates sample data.
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) Right-click in the script editor and select [Run Script] from the pop-up menu.
0-5) A JMP data table with the above data values will be created.
■①How to call R from JMP and run a modified Poisson regression
1) Install R and check that it can be called from JMP
First, to call R from JMP, R must be installed on the machine where JMP (JMP, JMP Pro, JMP Student Edition) is installed.
If R is already installed, or after installing R, check whether R can be called from JMP. The following explains how to do this using the menu operations in the Windows version.
1-1) Copy and paste the following code into the script editor:
R Init();
R Submit("
R.Version();
print('Hello, World!')
");
1-2) Right-click in the script editor and select [Run Script] from the pop-up menu.
1-3) If the log below the script editor displays the R version and a string such as "Hello, World!", then you have successfully called R from JMP.
2) Next, send the JMP data table to R, which is connected to JMP by the R Init() function above. You can send a JMP data table to R using the R Send function.
2-1) With the data table open, run the following code (if you are already connected to R, you do not need to run the R Init() at the beginning).
R Init();
dt = Current Data Table();
dt << New Column("ID", Nominal, Formula(Row()));
Wait(0);
R Send(dt);
Current Data Table() obtains a reference to the JMP data table that is currently the current data table.
"dt << New Column("ID", Nominal, Formula(Row()));" creates a column of sequential numbers for each row. If you already have a column of sequential numbers (ID numbers), this line is not necessary.
The last R Send(dt) sends the data table to R.
3) Run the R code.
Once the data is sent to R, the next step is to estimate a modified Poisson regression in R. Below, we first show an example using the geeglm function from the geepack library, and then the rqlm function from the rqlm package.
GEE (Generalized Equation Estimation) was originally proposed as an estimation method for correlated data such as frequency data, and sandwich variance was used to calculate standard errors. This sandwich variance calculation function can be used even for data where each row is independent. To do this, use an ID column that is unique for each row (and use independence as the correlation structure). For the reasons stated above, specify the Poisson distribution as the probability distribution of the response variable and the multilogarithm as the link function.
The following code uses geeglm to perform a modified Poisson regression with X1 and X2 as explanatory variables (in GEE terms, a log link function, a Poisson distribution, and sandwich variances, with each row considered independent). Parameter estimates and confidence intervals are exponentially transformed.
The results are output to the log at the bottom of the script window.
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)))
");
Below is an example using the rqlm package. This was programmed with reference to the following external materials. I think it is easier to specify than the geepack example above.
Noma, Hisashi (2025)
rqlm: An R package for performing modified Poisson regression and modified least squares regression analyses
/* 以下の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
");
■①How to call Python from JMP and run a modified Poisson regression
Starting with JMP 18, installing JMP also installs Python, which runs only within JMP, so there is no need to install Python separately.
1-0) The following will be written using Python code. To create Python code in JMP, select File (F) > New (N) > Python Code. You can also right-click in the Python script editor and select "Show log in window" to display the log window at the bottom.
1-1) Install some necessary Python packages. In this example, the required packages are numpy, pandas, and statsmodels. Copy and paste the following code into the Python script editor and run it.
from jmputils import jpip
jpip('install --upgrade', 'pip setuptools certifi')
jpip('install', 'numpy pandas statsmodels')
This code only needs to be executed once and there is no need to execute it again.
1-2) Import the package or function.
Copy and paste the following code into your Python script editor and run it.
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)
Read a JMP data table into Python and convert it to a pandas data frame.
In this example, the current JMP data table (the JMP data table that is currently open) is read and the data is converted to a pandas data frame.
The code below needs to be modified depending on your data.
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]
}
)
To use Python's GEE function, it seems necessary to create indicator variables (dummy variables) yourself.
In this example, X1 is a binary variable. Here, we create a dummy variable that is 1 if X1=="B" and 0 otherwise.
data['X1B'] = (data['X1'] == 'B').astype(int)
print(data)
1-4)
Run GEE. Specify Poisson distribution as the probability distribution of the response variable and logarithm as the link function.
Additionally, parameter estimates and their confidence intervals were exponentially transformed.
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()))
■③How to program modified Poisson regression yourself from scratch using the JMP scripting language
Modified Poisson regression first finds point estimates using maximum likelihood estimation of Poisson regression.
In this case, the Newton-Raphson method is often used. The Newton-Raphson method is an iterative calculation, and for each iteration, the log-likelihood function needs to be differentiated once with respect to the parameters (slope, gradient) and twice (Hessian matrix). Poisson regression belongs to a family of models called "generalized linear models," which make it easier to calculate the slope and Hessian matrix.
Also, if you are programming from scratch yourself, you will need to create your own indicator variables (dummy variables) for nominal scale variables.
I won’t go into the details, but if you ignore error handling (how to handle missing values, how to handle rank deficiency, how to handle separation), you can perform modified Poisson regression with the following amount of code:
Names Default To Here(1);
dt = Current Data Table();
/* 自分で指示変数（ダミー変数）を作成しないといけない */
x1 = Design(Column(dt, "X1")<
<GET as="" matrix=""></GET>
<DIV>I apologize for the lengthy explanation. To do this in JMP Pro 18 (JMP Student Edition), you'll need to do the above amount of programming. The JMP scripting language also allows you to create GUI-equipped applications. Some features not available in JMP are available on the JMP Marketplace (extensions offered there are provided as-is). Unfortunately, the modified Poisson regression feature is not available on the JMP Marketplace. If you can wait, I recommend waiting until JMP Pro 19 is released.</DIV>
<DIV>
</DIV>
<DIV>Whether or not a certain process can be performed in JMP can be divided into the following levels. In this case, I think it corresponds to level d).</DIV>
<DIV>a) Functions commonly available in JMP and JMP Pro → Can be done in no time</DIV>
<DIV>
<DIV>b) If you use the functions of JMP or JMP Pro, you can do it quickly with some ingenuity.</DIV>
</DIV>
<DIV>
<DIV>c) Features provided on the JMP Marketplace → As is, but ready to use</DIV>
<DIV>d) A simple function that is not provided above → It can be done by programming (especially in R or Python)</DIV>
<DIV>e) Highly advanced functions not provided above → Only possible with considerable expertise</DIV>
<DIV>
</DIV>
<DIV>Finally, there's one thing that bothered me, so I'd like to add a comment about the "correctness" I mentioned at the beginning. Regression analysis itself, including the simplest least-squares linear regression, has a variety of uses. For example, the following paper divides it into description, prediction, and control/intervention.</DIV>
<DIV>Iwasaki, M. (2021)</DIV>
<DIV>Multiple regression analysis from the perspective of statistical causal inference</DIV>
<DIV>Journal of the Japan Statistical Society, 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>The explanation above merely describes how to perform modified Poisson regression, and does not in any way establish "correctness" in the sense of whether it is "correct" to perform regression analysis in the current situation, whether performing modified Poisson regression is "correct," or even whether the current situation is "correct" in the sense that it satisfies the major assumptions made by regression analysis. Whether your current statistical analysis is "correct" for your actual situation (as I myself have no ability to judge) is something I would recommend consulting an applied statistician in that field (I have never consulted with one, but I believe they may be able to suggest alternatives to the analysis itself and also provide detailed instructions on how to perform modified Poisson regression).</DIV>
<DIV>
</DIV>
</DIV>
