Gracias por tu respuesta.
Sin embargo, la regresión de Poisson modificada calcula los intervalos de confianza correctos utilizando una varianza robusta, ¿verdad?
Sí. Para ser más específico sobre "correcto", creo que diría que es correcto cuando el tamaño de la muestra es grande. Aún más específico, creo que diría que los resultados estimados utilizando un modelo que asume una distribución de Poisson, aunque la variable de respuesta no siga realmente una distribución de Poisson, serán correctos si el tamaño de la muestra es grande. En jerga estadística, creo que podría decir que "independientemente de cuál sea la distribución de probabilidad de la variable de respuesta, la varianza sándwich es un estimador consistente para la matriz de varianza-covarianza del estimador de un modelo que asume un modelo específico (en este caso, una distribución de Poisson)".
¿Estoy en lo cierto al entender que la varianza robusta y la varianza sándwich se refieren a lo mismo?
Sí. Al menos en este contexto, me refiero exactamente a lo mismo con "varianza robusta" y "varianza sándwich". Otro término es "varianza empírica". Esta varianza sándwich puede utilizarse para la estimación M (en ciertas condiciones). La estimación de máxima verosimilitud también se incluye en la estimación M.
■Como usted dijo, la regresión de Poisson modificada primero estima la regresión de Poisson usando máxima verosimilitud, luego calcula el error estándar usando varianza robusta y luego usa el error estándar intercalado para calcular pruebas e intervalos de confianza.
En la estimación de máxima verosimilitud de la regresión de Poisson, se asume una distribución de Poisson y, por lo tanto, no una distribución binomial. Al utilizar la varianza robusta, incluso si la distribución de probabilidad real de la variable de respuesta no es una distribución de Poisson, pero el modelo la asume, la estimación es válida siempre que el tamaño de la muestra sea grande.
■En epidemiología en particular, creo que la regresión de Poisson modificada es un método propuesto con las siguientes motivaciones.
Si la probabilidad de que un evento ocurra es p, un modelo comúnmente utilizado es el modelo logístico. Para simplificar, si solo hay una variable explicativa, la fórmula del modelo es la siguiente:
Si el lado izquierdo de la ecuación cuando x=1 es log(p1/(1-p1)), y el lado izquierdo cuando x=0 es log(p0/(1-p0)), entonces β1 es
β1 = log(p1/(1-p1)) - log(p0/(1-p0))
Por las propiedades de los logaritmos, esto es
β1 = log({p1/(1-p1)}/{p0/(1-p0)})
Al transformar exponencialmente ambos lados, obtenemos exp(log(a))=a,
exp(β1) = {p1/(1-p1)}/{p0/(1-p0)}
En otras palabras, β1 representa el logaritmo de la razón de probabilidades (cuando x cambia en 1), y exp(β1) representa la razón de probabilidades (cuando x cambia en 1). Esto es correcto, pero a veces puede ser difícil interpretar la razón de probabilidades. Si p1 y p0 son cercanos a cero,
{p1/(1-p1)}/{p0/(1-p0)}≒p1/p0
Por lo tanto, la razón de probabilidades (odds ratio) puede aproximarse mediante la razón de riesgos (razón de proporciones). Sin embargo, esta aproximación solo es válida si p1 y p0 son cercanos a cero. Una de las ventajas de la regresión logística es que p es mayor que 0 y menor que 1 cuando el lado derecho del modelo, β0+β1x, se encuentra entre -∞ y +∞. En otras palabras, independientemente del aspecto del lado derecho del modelo, p casi siempre conservará las propiedades de una probabilidad (entre 0 y 1). (Sin embargo, en un modelo de regresión logística, surgen problemas cuando la probabilidad es exactamente 0 o exactamente 1).
Dado que los odds ratios son difíciles de interpretar, el modelo correspondiente a los ratios de riesgo es fácil de interpretar para los analistas.
En otras palabras, el lado izquierdo es el logaritmo de la proporción en lugar del logit.
(Las fórmulas de transformación como logit (log(p/(1-p))) y logaritmo (log(p)) se denominan funciones de enlace).
En este caso, log(p1) - log(p0) = β1, por lo que
log(p1/p0) = β1
p1/p0 = exp(β1)
En otras palabras, en este modelo con una función de enlace logarítmica, β1 es la razón de riesgo logarítmica y exp(β1) es la razón de riesgo.
Sin embargo, este modelo presenta un problema. Si el lado derecho se convierte, por ejemplo, en β0 + β1x = 1.0, entonces log(p) = 1.0, o p = exp(1.0) ≒ 2.72, y aunque es una probabilidad, p a veces puede ser mayor que 1. Como resultado, si se realiza una estimación de máxima verosimilitud basada en la probabilidad de una distribución binomial, el cálculo puede fallar a mitad de camino.
Para evitar este problema, se podría considerar realizar una estimación de máxima verosimilitud basada en la distribución de Poisson, que funciona bien cuando la función de enlace es logarítmica. Sin embargo, al hacerlo, se calcularían errores estándar, pruebas e intervalos de confianza utilizando la distribución de Poisson (aunque en realidad se cree que siguen una distribución binomial). Por lo tanto, se propuso utilizar la varianza sándwich (varianza robusta/varianza empírica), que es robusta al especificar la distribución de probabilidad de la variable de respuesta. Cuando se especifica la varianza sándwich, incluso si se especifica una distribución de probabilidad incorrecta para la variable de respuesta, los errores estándar sándwich (y la matriz de covarianza sándwich) para los estimadores de parámetros serán razonables si el tamaño de la muestra es grande (son estimadores consistentes de los errores estándar reales).
Por cierto, probablemente haya dos razones por las que las razones de probabilidades son difíciles de interpretar. La primera es que son intuitivamente difíciles. La razón "a ÷ b" es fácil de entender, pero "(a ÷ (1-a)) ÷ (b ÷ (1-b))" es un poco difícil de entender intuitivamente. Además, mientras que las diferencias de riesgo y las razones de riesgo tienen una propiedad llamada colapsabilidad, las razones de probabilidades no la tienen.
Como se mencionó anteriormente, la regresión de Poisson modificada está ganando popularidad porque compensa la dificultad de interpretación de las razones de probabilidades. Por lo tanto, creo que, si estas no son particularmente difíciles de interpretar, sería adecuado seguir con el modelo de regresión logística tradicional. (Por ejemplo, si el propósito del análisis es la "predicción" o la "clasificación" y la interpretación de los coeficientes de regresión no es importante, creo que la regresión logística es más fácil de usar que la regresión de Poisson modificada).
Los métodos que mencioné el otro día (llamar a R o Python y programar desde cero con el lenguaje de script JMP) requieren programación, lo que elimina lo mejor de JMP (disfrutar del análisis estadístico interactivo sin programar). Sin embargo, incluso entre los usuarios de JMP, especialmente los que realizan análisis estadísticos rutinarios, la programación se realiza en JMP.
Con JMP Pro 19, podrá realizar una regresión de Poisson modificada con operaciones con el ratón. Actualmente, para hacerlo en JMP, necesitará programar un poco. Esto es mucho más complejo que usar JMP convencional.
A continuación se muestran unos pequeños pseudodatos como datos de muestra:
①Cómo llamar a R desde JMP y ejecutar una regresión de Poisson modificada
②Cómo ejecutar una regresión de Poisson modificada llamando a Python desde JMP
3) Cómo programar usted mismo una regresión de Poisson modificada desde cero utilizando el lenguaje de script JMP
Primero, me gustaría que le echaran un vistazo rápido para hacerse una idea aproximada de cuánto trabajo requiere. El código a continuación es solo una muestra creada por mí y no es, en ningún caso, un producto proporcionado por nuestra empresa (soy el único responsable, pero al ser una muestra, no ofrece garantía alguna).
Además, la siguiente explicación está diseñada para funcionar únicamente con los pequeños pseudodatos presentados aquí.
Aunque aquí utilizamos conjuntos de datos bastante pequeños a efectos de explicación, la varianza sándwich es en realidad un método válido para cantidades relativamente grandes de datos (es un método basado en la teoría asintótica).
Con respecto a 1 y 2, si los nombres de las columnas contienen caracteres japoneses, podrían no funcionar correctamente. Además, no estoy familiarizado con R ni Python, así que podría haber paquetes o métodos de programación más sofisticados. Respecto a la programación en Python que se describe a continuación, consulté ampliamente las sugerencias de Copilot.
Respecto al punto ③, los puntos siguientes no se tienen en cuenta en absoluto.
a) Manejo de valores faltantes: los valores faltantes en los datos generarán un error.
b) Manejo de errores por deficiencia de rango en la matriz de diseño: En modelos lineales generalizados, como el análisis de regresión lineal y la regresión de Poisson, si existe dependencia lineal entre las variables explicativas, las estimaciones de los parámetros no pueden determinarse de forma unívoca. No se realiza ningún procesamiento para abordar este problema.
c) Manejo de errores de separación: La regresión logística y la regresión de Poisson pueden presentar un problema llamado separación según los datos, lo que puede provocar que las estimaciones de los parámetros sean +∞ o -∞. Este problema no se aborda. Respecto al punto 3, priorizamos la legibilidad y hemos programado sin preocuparnos por la eficiencia del procesamiento (especialmente el consumo de memoria).
He adjuntado un archivo con el código JSL que combina los pasos 1 y 3. Para ejecutarlo, debes tener R instalado y poder llamarlo desde JMP.
También he adjuntado el archivo de código Python del paso 2. Para ejecutarlo, debe tener abierta una tabla de datos JMP, como se indica en las líneas 1 a 7 del archivo JSL anterior. Además, para ejecutarlo por primera vez, debe ejecutar las líneas 2 a 4.
■0) Primero, cree los datos de muestra que se utilizarán en ①, ② y ③ a continuación.
0-1) Seleccione Archivo (F) > Nuevo (N) > Script JSL (Ctrl + T).
0-2) Haga clic derecho en el editor de scripts y seleccione [Mostrar ventana de registro] en el menú emergente. Esto mostrará la ventana de registro en la parte inferior.
0-3) Copie y pegue el siguiente código en la parte superior del editor de scripts. Este código crea datos de muestra.
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) Haga clic derecho en el editor de scripts y seleccione [Ejecutar script] en el menú emergente.
0-5) Se creará una tabla de datos JMP con los valores de datos anteriores.
■①Cómo llamar a R desde JMP y ejecutar una regresión de Poisson modificada
1) Instale R y verifique que se pueda llamar desde JMP
Primero, para llamar a R desde JMP, R debe estar instalado en la máquina donde está instalado JMP (JMP, JMP Pro, JMP Student Edition).
Si R ya está instalado, o después de instalarlo, compruebe si se puede llamar desde JMP. A continuación, se explica cómo hacerlo mediante las operaciones de menú de la versión de Windows.
1-1) Copie y pegue el siguiente código en el editor de scripts:
R Init();
R Submit("
R.Version();
print('Hello, World!')
");
1-2) Haga clic derecho en el editor de scripts y seleccione [Ejecutar script] en el menú emergente.
1-3) Si el registro debajo del editor de scripts muestra la versión de R y una cadena como "¡Hola, mundo!", entonces ha llamado exitosamente a R desde JMP.
2) A continuación, envíe la tabla de datos de JMP a R, que está conectada a JMP mediante la función R Init() mencionada anteriormente. Puede enviar tablas de datos de JMP a R mediante la función R Send.
2-1) Con la tabla de datos abierta, ejecute el siguiente código (si ya está conectado a R, no necesita ejecutar R Init() al principio).
R Init();
dt = Current Data Table();
dt << New Column("ID", Nominal, Formula(Row()));
Wait(0);
R Send(dt);
La tabla de datos actual() obtiene una referencia a la tabla de datos JMP que actualmente es la tabla de datos actual.
"dt << New Column("ID", Nominal, Formula(Row()));" crea una columna de números secuenciales para cada fila. Si ya tiene una columna de números secuenciales (números de ID), esta línea no es necesaria.
El último R Send(dt) envía la tabla de datos a R.
3) Ejecute el código R.
Una vez que los datos se envían a R, el siguiente paso es estimar una regresión de Poisson modificada en R. A continuación, primero mostramos un ejemplo utilizando la función geeglm de la biblioteca geepack y luego la función rqlm del paquete rqlm.
La GEE (Estimación de Ecuaciones Generalizadas) se propuso originalmente como método de estimación para datos correlacionados, como datos de frecuencia, y se utilizó la varianza sándwich para calcular los errores estándar. Esta función de cálculo de la varianza sándwich puede utilizarse incluso para datos donde cada fila es independiente. Para ello, utilice una columna de ID única para cada fila (y utilice la independencia como estructura de correlación). Por las razones expuestas anteriormente, especifique la distribución de Poisson como la distribución de probabilidad de la variable de respuesta y el multilogaritmo como la función de enlace.
El siguiente código utiliza geeglm para realizar una regresión de Poisson modificada con X1 y X2 como variables explicativas (en términos de GEE, una función de enlace logarítmico, una distribución de Poisson y varianzas sándwich, donde cada fila se considera independiente). Las estimaciones de los parámetros y los intervalos de confianza se transforman exponencialmente.
Los resultados se muestran en el registro en la parte inferior de la ventana del script.
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)))
");
A continuación se muestra un ejemplo con el paquete rqlm. Este se programó con referencia a los siguientes materiales externos. Creo que es más fácil de especificar que el ejemplo de geepack mencionado anteriormente.
Noma, Hisashi (2025)
rqlm: Un paquete R para realizar análisis de regresión de Poisson modificada y regresión de mínimos cuadrados modificada
/* 以下の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
");
■①Cómo llamar a Python desde JMP y ejecutar una regresión de Poisson modificada
A partir de JMP 18, la instalación de JMP también instala Python, que se ejecuta solo dentro de JMP, por lo que no es necesario instalar Python por separado.
1-0) Lo siguiente se escribirá con código Python. Para crear código Python en JMP, seleccione Archivo (F) > Nuevo (N) > Código Python. También puede hacer clic derecho en el editor de scripts de Python y seleccionar "Mostrar registro en la ventana" para mostrar la ventana de registro en la parte inferior.
1-1) Instale los paquetes de Python necesarios. En este ejemplo, los paquetes requeridos son numpy, pandas y statsmodels. Copie y pegue el siguiente código en el editor de scripts de Python y ejecútelo.
from jmputils import jpip
jpip('install --upgrade', 'pip setuptools certifi')
jpip('install', 'numpy pandas statsmodels')
Este código solo necesita ejecutarse una vez y no es necesario ejecutarlo nuevamente.
1-2) Importar el paquete o función.
Copie y pegue el siguiente código en su editor de scripts de Python y ejecútelo.
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)
Leer una tabla de datos JMP en Python y convertirla en un marco de datos de pandas.
En este ejemplo, se lee la tabla de datos JMP actual (la tabla de datos JMP que está abierta actualmente) y los datos se convierten en un marco de datos de pandas.
El código siguiente debe modificarse dependiendo de sus datos.
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]
}
)
Para utilizar la función GEE de Python, parece necesario crear usted mismo variables indicadoras (variables ficticias).
En este ejemplo, X1 es una variable binaria. Creamos una variable ficticia que es 1 si X1=="B" y 0 en caso contrario.
data['X1B'] = (data['X1'] == 'B').astype(int)
print(data)
1-4)
Ejecute GEE. Especifique la distribución de Poisson como la distribución de probabilidad de la variable de respuesta y el logaritmo como la función de enlace.
Además, las estimaciones de los parámetros y sus intervalos de confianza se transformaron exponencialmente.
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()))
■③Cómo programar usted mismo una regresión de Poisson modificada desde cero utilizando el lenguaje de script JMP
La regresión de Poisson modificada primero obtiene estimaciones puntuales utilizando la estimación de máxima verosimilitud de la regresión de Poisson.
En este caso, se suele utilizar el método de Newton-Raphson. Este método es un cálculo iterativo, y para cada iteración, la función de log-verosimilitud debe diferenciarse una vez con respecto a los parámetros (pendiente, gradiente) y dos veces (matriz hessiana). La regresión de Poisson pertenece a una familia de modelos denominados "modelos lineales generalizados", que facilitan el cálculo de la pendiente y la matriz hessiana.
Además, si está programando desde cero usted mismo, necesitará crear sus propias variables indicadoras (variables ficticias) para las variables de escala nominal.
No entraré en detalles, pero si ignora el manejo de errores (cómo manejar valores faltantes, cómo manejar la deficiencia de rango, cómo manejar la separación), puede realizar una regresión de Poisson modificada con la siguiente cantidad de código:
Names Default To Here(1);
dt = Current Data Table();
/* 自分で指示変数（ダミー変数）を作成しないといけない */
x1 = Design(Column(dt, "X1")<
<GET as="" matrix=""></GET>
<DIV>Disculpen la extensión de la explicación. Para hacer esto en JMP Pro 18 (JMP Student Edition), necesitarán realizar la programación descrita anteriormente. El lenguaje de scripting JMP también permite crear aplicaciones con interfaz gráfica de usuario (GUI). Algunas funciones no disponibles en JMP están disponibles en JMP Marketplace (las extensiones que se ofrecen allí se proporcionan tal cual). Lamentablemente, la función de regresión de Poisson modificada no está disponible en JMP Marketplace. Si pueden esperar, les recomiendo esperar hasta el lanzamiento de JMP Pro 19.</DIV>
<DIV>
</DIV>
<DIV>La posibilidad de que un proceso determinado se realice en JMP se puede dividir en los siguientes niveles. En este caso, creo que corresponde al nivel d).</DIV>
<DIV>a) Funciones comúnmente disponibles en JMP y JMP Pro → Se pueden realizar en poco tiempo</DIV>
<DIV>
<DIV>b) Si utilizas las funciones de JMP o JMP Pro, podrás hacerlo rápidamente con algo de ingenio.</DIV>
</DIV>
<DIV>
<DIV>c) Características proporcionadas en JMP Marketplace → Tal como están, pero listas para usar</DIV>
<DIV>d) Una función simple que no se proporciona arriba → Se puede realizar mediante programación (especialmente en R o Python)</DIV>
<DIV>e) Funciones muy avanzadas no previstas anteriormente → Solo es posible con considerable experiencia</DIV>
<DIV>
</DIV>
<DIV>Finalmente, hay algo que me inquietó, así que quisiera añadir un comentario sobre la "corrección" que mencioné al principio. El análisis de regresión en sí, incluyendo la regresión lineal de mínimos cuadrados más simple, tiene diversos usos. Por ejemplo, el siguiente artículo lo divide en descripción, predicción y control/intervención.</DIV>
<DIV>Iwasaki, M. (2021)</DIV>
<DIV>Análisis de regresión múltiple desde la perspectiva de la inferencia causal estadística</DIV>
<DIV>Revista de la Sociedad Japonesa de Estadística, 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>La explicación anterior simplemente describe cómo realizar la regresión de Poisson modificada y no establece en modo alguno la "corrección" en el sentido de si es "correcto" realizar el análisis de regresión en la situación actual, si realizar la regresión de Poisson modificada es "correcto", ni siquiera si la situación actual es "correcta" en el sentido de que cumple con los supuestos principales del análisis de regresión. Para determinar si su análisis estadístico actual es "correcto" para su situación real (ya que yo mismo no tengo la capacidad de juzgarlo), recomiendo consultar con un experto en estadística aplicada en ese campo (nunca he consultado con uno, pero creo que podría sugerir alternativas al análisis en sí y también proporcionar instrucciones detalladas sobre cómo realizar la regresión de Poisson modificada).</DIV>
<DIV>
</DIV>
</DIV>
