Grazie per la risposta.
Tuttavia, la regressione di Poisson modificata calcola gli intervalli di confidenza corretti utilizzando una varianza robusta, giusto?
Sì. Per essere più specifici sul termine "corretto", direi che è corretto quando la dimensione del campione è ampia. Ancora più specificamente, direi che i risultati stimati utilizzando un modello che presuppone una distribuzione di Poisson, anche se la variabile di risposta non segue effettivamente una distribuzione di Poisson, saranno corretti se la dimensione del campione è ampia. Usando il gergo statistico, direi che "indipendentemente dalla distribuzione di probabilità della variabile di risposta, la varianza sandwich è uno stimatore coerente per la matrice varianza-covarianza dello stimatore di un modello che presuppone un modello specifico (in questo caso, una distribuzione di Poisson)".
Ho capito bene che la varianza robusta e la varianza sandwich si riferiscono alla stessa cosa?
Sì. Almeno in questo contesto, intendo esattamente la stessa cosa con "varianza robusta" e "varianza sandwich". Un altro termine è "varianza empirica". Questa varianza sandwich può essere utilizzata per la stima M (in determinate condizioni). Anche la stima di massima verosimiglianza è inclusa nella stima M.
■Come hai detto, la regressione di Poisson modificata esegue prima la stima della massima verosimiglianza della regressione di Poisson, quindi calcola l'errore standard utilizzando la varianza robusta e infine calcola il test e l'intervallo di confidenza dall'errore standard inserito.
Nella stima di massima verosimiglianza della regressione di Poisson, si assume una distribuzione di Poisson e, pertanto, non si assume una distribuzione binomiale. Utilizzando la varianza robusta, anche se la distribuzione di probabilità reale della variabile di risposta non è una distribuzione di Poisson, ma il modello assume una distribuzione di Poisson, la stima è valida finché la dimensione del campione è ampia.
■In particolare in epidemiologia, ritengo che la regressione di Poisson modificata sia un metodo proposto con le seguenti motivazioni.
Se la probabilità che un evento si verifichi è p, un modello comunemente utilizzato è il modello logistico. Per semplicità, se è presente una sola variabile esplicativa, la formula del modello è la seguente:
Se il lato sinistro dell'equazione quando x=1 è log(p1/(1-p1)), e il lato sinistro quando x=0 è log(p0/(1-p0)), allora β1 è
β1 = log(p1/(1-p1)) - log(p0/(1-p0))
Per le proprietà dei logaritmi, questo è
β1 = log({p1/(1-p1)}/{p0/(1-p0)})
Trasformando esponenzialmente entrambi i lati, otteniamo exp(log(a))=a,
exp(β1) = {p1/(1-p1)}/{p0/(1-p0)}
In altre parole, β1 rappresenta il logaritmo del rapporto di probabilità (quando x varia di 1), ed exp(β1) rappresenta il rapporto di probabilità (quando x varia di 1). Questo va bene, ma a volte può essere difficile interpretare il rapporto di probabilità. Se p1 e p0 sono prossimi a zero,
{p1/(1-p1)}/{p0/(1-p0)}≒p1/p0
Pertanto, l'odds ratio può essere approssimato dal risk ratio (rapporto di proporzioni). Tuttavia, questa approssimazione è valida solo se p1 e p0 sono prossimi a zero. Uno dei vantaggi della regressione logistica è che p è maggiore di 0 e minore di 1 quando il lato destro del modello, β0+β1x, si muove tra -∞ e +∞. In altre parole, indipendentemente dall'aspetto del lato destro del modello, p manterrà quasi sempre le proprietà di una probabilità (tra 0 e 1). (Tuttavia, in un modello di regressione logistica, sorgono problemi quando la probabilità è esattamente 0 o esattamente 1.)
Poiché i rapporti di probabilità sono difficili da interpretare, il modello corrispondente ai rapporti di rischio è facile da interpretare per gli analisti.
In altre parole, il lato sinistro è il logaritmo della proporzione anziché il logit.
(Le formule di trasformazione come logit (log(p/(1-p))) e logaritmo (log(p)) sono chiamate funzioni di collegamento.)
In questo caso, log(p1) - log(p0) = β1, quindi
log(p1/p0) = β1
p1/p0 = exp(β1)
In altre parole, in questo modello con una funzione di collegamento logaritmica, β1 è il logaritmo del rapporto di rischio e exp(β1) è il rapporto di rischio.
Tuttavia, questo modello presenta un problema. Se il membro destro diventa, ad esempio, β0 + β1x = 1,0, allora log(p) = 1,0, ovvero p = exp(1,0) ≒ 2,72, e sebbene si tratti di una probabilità, p può talvolta superare 1. Di conseguenza, se si esegue una stima di massima verosimiglianza basata sulla verosimiglianza di una distribuzione binomiale, il calcolo potrebbe fallire a metà.
Per evitare questo problema, si potrebbe considerare di eseguire una stima di massima verosimiglianza basata sulla distribuzione di Poisson, che funziona bene quando la funzione di collegamento è logaritmica. Tuttavia, ciò comporterebbe il calcolo di errori standard, test e intervalli di confidenza utilizzando la distribuzione di Poisson (anche se in realtà si ritiene che seguano una distribuzione binomiale). Pertanto, è stato proposto di utilizzare la varianza sandwich (varianza robusta/varianza empirica), che è robusta quando si specifica la distribuzione di probabilità della variabile di risposta. Quando si specifica la varianza sandwich, anche se si specifica una distribuzione di probabilità errata per la variabile di risposta, gli errori standard sandwich (e la matrice di covarianza sandwich) per gli stimatori dei parametri saranno ragionevoli se la dimensione del campione è ampia (sono stimatori coerenti degli errori standard reali).
Per inciso, ci sono probabilmente due ragioni per cui gli odds ratio sono difficili da interpretare. La prima è che sono intuitivamente difficili. Il rapporto "a ÷ b" è facile da capire, ma "(a ÷ (1-a)) ÷ (b ÷ (1-b))" è un po' difficile da comprendere intuitivamente. Inoltre, mentre le differenze di rischio e i risk ratio hanno una proprietà chiamata collassabilità, gli odds ratio non hanno questa proprietà.
Come accennato in precedenza, la regressione di Poisson modificata sta guadagnando attenzione perché compensa lo svantaggio della difficoltà di interpretazione degli odds ratio, quindi personalmente ritengo che, se gli odds ratio non sono particolarmente difficili da interpretare, sarebbe opportuno attenersi al modello di regressione logistica tradizionale. (Ad esempio, se lo scopo dell'analisi è la "previsione" o la "classificazione" e l'interpretazione dei coefficienti di regressione non è importante, ritengo che la regressione logistica sia più facile da usare rispetto alla regressione di Poisson modificata.)
■I metodi che ho menzionato l'altro giorno (chiamare R o Python e programmare da zero usando il linguaggio di scripting JMP) richiedono programmazione, il che elimina il vantaggio migliore di JMP (godersi un'analisi statistica interattiva senza dover scrivere codice). Tuttavia, anche tra gli utenti JMP, soprattutto tra gli utenti più assidui che eseguono elaborazioni statistiche di routine, la programmazione viene eseguita in JMP.
Con JMP Pro 19, sarà possibile eseguire una regressione di Poisson modificata con operazioni tramite mouse. Attualmente, per farlo in JMP, è necessario un po' di programmazione. Si tratta di un percorso molto più complesso rispetto all'utilizzo di JMP standard.
Di seguito sono riportati alcuni piccoli pseudodati come dati campione:
①Come chiamare R da JMP ed eseguire una regressione di Poisson modificata
②Come eseguire una regressione di Poisson modificata chiamando Python da JMP
3) Come programmare da zero una regressione di Poisson modificata utilizzando il linguaggio di scripting JMP
Innanzitutto, vorrei che gli dessi una rapida occhiata per farti un'idea approssimativa di quanto lavoro sia stato necessario. Il codice qui sotto è solo un esempio creato da me personalmente e non è in alcun modo un prodotto fornito dalla nostra azienda come organizzazione (me ne assumo la piena responsabilità, ma trattandosi di un esempio, non offre alcuna garanzia).
Inoltre, la seguente spiegazione è concepita per funzionare solo con i piccoli pseudo-dati qui presentati.
Sebbene qui utilizziamo set di dati piuttosto piccoli a scopo esplicativo, la varianza sandwich è in realtà un metodo valido per quantità di dati relativamente grandi (è un metodo basato sulla teoria asintotica).
Per quanto riguarda ① e ②, se i nomi delle colonne contengono caratteri giapponesi, potrebbero non funzionare correttamente. Inoltre, non ho familiarità con R o Python, quindi potrebbero esserci pacchetti o metodi di programmazione più eleganti. Per quanto riguarda la programmazione Python qui sotto, ho fatto ampio riferimento ai suggerimenti di Copilot.
Per quanto riguarda ③, i seguenti punti non vengono presi in considerazione.
a) Gestione dei valori mancanti: i valori mancanti nei dati causeranno un errore.
b) Gestione degli errori per la carenza di rango nella matrice di progettazione: nei modelli lineari generalizzati come l'analisi di regressione lineare e la regressione di Poisson, se esiste una dipendenza lineare tra variabili esplicative, le stime dei parametri non possono essere determinate in modo univoco. Non viene eseguita alcuna elaborazione per risolvere questo problema.
c) Gestione degli errori per la separazione: la regressione logistica e la regressione di Poisson possono incontrare un problema chiamato separazione a seconda dei dati, che può causare stime dei parametri di +∞ o -∞. Questo problema non viene gestito. Per quanto riguarda il punto 3, diamo importanza alla leggibilità e abbiamo programmato senza alcuna preoccupazione per l'efficienza di elaborazione (in particolare il consumo di memoria).
Ho allegato un file con il codice JSL che combina i passaggi 1 e 3. Per eseguirlo, è necessario avere R installato ed essere in grado di richiamarlo da JMP.
Ho anche allegato il file di codice Python per il passaggio 2. Per eseguirlo, è necessario avere una tabella dati JMP aperta, come indicato nelle righe 1-7 del file JSL sopra. Inoltre, per eseguirlo per la prima volta, è necessario eseguire le righe 2-4.
■0) Per prima cosa, creare i dati campione da utilizzare nei punti ①, ② e ③ di seguito.
0-1) Selezionare File (F) > Nuovo (N) > Script JSL (Ctrl + T).
0-2) Fare clic con il pulsante destro del mouse all'interno dell'editor di script chiamato e selezionare [Mostra finestra di registro] dal menu a comparsa. Verrà visualizzata la finestra di registro in basso.
0-3) Copia e incolla il seguente codice nella parte superiore dell'editor di script. Questo codice crea dati di esempio.
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) Fare clic con il pulsante destro del mouse nell'editor di script e selezionare [Esegui script] dal menu a comparsa.
0-5) Verrà creata una tabella dati JMP con i valori dei dati sopra indicati.
■①Come chiamare R da JMP ed eseguire una regressione di Poisson modificata
1) Installa R e verifica che possa essere chiamato da JMP
Innanzitutto, per chiamare R da JMP, R deve essere installato sulla macchina in cui è installato JMP (JMP, JMP Pro, JMP Student Edition).
Se R è già installato, o dopo averlo installato, verificare se R può essere richiamato da JMP. Di seguito viene spiegato come farlo utilizzando le operazioni di menu per la versione Windows.
1-1) Copia e incolla il seguente codice nell'editor di script:
R Init();
R Submit("
R.Version();
print('Hello, World!')
");
1-2) Fare clic con il tasto destro del mouse nell'editor di script e selezionare [Esegui script] dal menu a comparsa.
1-3) Se il registro sotto l'editor di script visualizza la versione di R e una stringa come "Hello, World!", significa che hai chiamato correttamente R da JMP.
2) Successivamente, invia la tabella dati JMP a R, che è connessa a JMP tramite la funzione R Init() descritta sopra. Puoi inviare una tabella dati JMP a R utilizzando la funzione R Send.
2-1) Con la tabella dati aperta, esegui il seguente codice (se sei già connesso a R, non è necessario eseguire R Init() all'inizio).
R Init();
dt = Current Data Table();
dt << New Column("ID", Nominal, Formula(Row()));
Wait(0);
R Send(dt);
Current Data Table() ottiene un riferimento alla tabella dati JMP che è attualmente la tabella dati corrente.
"dt << New Column("ID", Nominal, Formula(Row()));" crea una colonna di numeri sequenziali per ogni riga. Se si dispone già di una colonna di numeri sequenziali (numeri ID), questa riga non è necessaria.
L'ultimo R Send(dt) invia la tabella dati a R.
3) Eseguire il codice R.
Una volta inviati i dati a R, il passo successivo è stimare una regressione di Poisson modificata in R. Di seguito, mostriamo prima un esempio che utilizza la funzione geeglm della libreria geepack e poi la funzione rqlm del pacchetto rqlm.
La GEE (Generalized Equation Estimation) è stata originariamente proposta come metodo di stima per i casi in cui esiste una correlazione, come nel caso dei dati di frequenza, e la varianza sandwich è stata utilizzata per calcolare gli errori standard. Questa funzione di calcolo della varianza sandwich può essere utilizzata anche per dati in cui ogni riga è indipendente. A tale scopo, utilizzare una colonna ID univoca per ogni riga (e utilizzare l'indipendenza come struttura di correlazione). Per i motivi sopra indicati, specificare la distribuzione di Poisson come distribuzione di probabilità della variabile di risposta e il multilogaritmo come funzione di collegamento.
Il codice seguente utilizza geeglm per eseguire una regressione di Poisson modificata con X1 e X2 come variabili esplicative (in termini GEE, una funzione di collegamento logaritmica, una distribuzione di Poisson e varianze sandwich, con ogni riga considerata indipendente). Le stime dei parametri e gli intervalli di confidenza vengono trasformati esponenzialmente.
I risultati vengono visualizzati nel registro nella parte inferiore della finestra dello 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)))
");
Di seguito è riportato un esempio che utilizza il pacchetto rqlm. Questo è stato programmato facendo riferimento ai seguenti materiali esterni. Penso che sia più facile da specificare rispetto all'esempio geepack di cui sopra.
Noma, Hisashi (2025)
rqlm: un pacchetto R per eseguire analisi di regressione di Poisson modificata e di regressione dei minimi quadrati modificata
/* 以下の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
");
■①Come chiamare Python da JMP ed eseguire una regressione di Poisson modificata
A partire da JMP 18, l'installazione di JMP installa anche Python, che viene eseguito solo all'interno di JMP, quindi non è necessario installare Python separatamente.
1-0) Quanto segue verrà scritto utilizzando codice Python. Per creare codice Python in JMP, selezionare File (F) > Nuovo (N) > Codice Python. È anche possibile fare clic con il pulsante destro del mouse nell'editor di script Python e selezionare "Mostra finestra di registro" per visualizzare la finestra di registro in basso.
1-1) Installa alcuni pacchetti Python necessari. In questo esempio, i pacchetti richiesti sono numpy, pandas e statsmodels. Copia e incolla il seguente codice nell'editor di script Python ed eseguilo.
from jmputils import jpip
jpip('install --upgrade', 'pip setuptools certifi')
jpip('install', 'numpy pandas statsmodels')
Questo codice deve essere eseguito solo una volta e non è necessario eseguirlo di nuovo.
1-2) Importare il pacchetto o la funzione.
Copia e incolla il seguente codice nel tuo editor di script Python ed eseguilo.
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)
Leggere una tabella dati JMP in Python e convertirla in un data frame pandas.
In questo esempio, la tabella dati JMP corrente (la tabella dati JMP attualmente aperta) viene letta e i dati vengono convertiti in un data frame pandas.
Il codice seguente deve essere modificato in base ai dati.
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]
}
)
Per utilizzare la funzione GEE di Python, sembra necessario creare autonomamente delle variabili indicatrici (variabili fittizie).
In questo esempio, X1 è una variabile binaria. Creiamo una variabile fittizia che vale 1 se X1=="B" e 0 altrimenti.
data['X1B'] = (data['X1'] == 'B').astype(int)
print(data)
1-4)
Eseguire GEE. Specificare la distribuzione di Poisson come distribuzione di probabilità della variabile di risposta e il logaritmo come funzione di collegamento.
Inoltre, le stime dei parametri e i relativi intervalli di confidenza sono stati trasformati in modo esponenziale.
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()))
■③Come programmare autonomamente una regressione di Poisson modificata partendo da zero utilizzando il linguaggio di scripting JMP
La regressione di Poisson modificata trova innanzitutto stime puntuali utilizzando la stima di massima verosimiglianza della regressione di Poisson.
In questo caso, viene spesso utilizzato il metodo di Newton-Raphson. Il metodo di Newton-Raphson è un calcolo iterativo e, per ogni iterazione, la funzione di verosimiglianza logaritmica deve essere differenziata una volta rispetto ai parametri (pendenza, gradiente) e due volte rispetto alla matrice Hessiana. La regressione di Poisson appartiene a una famiglia di modelli chiamati "modelli lineari generalizzati", che semplificano il calcolo della pendenza e della matrice Hessiana.
Inoltre, se si programma da zero, sarà necessario creare le proprie variabili indicatrici (variabili fittizie) per le variabili di scala nominale.
Non entrerò nei dettagli, ma se si ignora la gestione degli errori (come gestire i valori mancanti, come gestire la carenza di rango, come gestire la separazione), è possibile eseguire una regressione di Poisson modificata con la seguente quantità di codice:
Names Default To Here(1);
dt = Current Data Table();
/* 自分で指示変数（ダミー変数）を作成しないといけない */
x1 = Design(Column(dt, "X1")<
<GET as="" matrix=""></GET>
<DIV>Mi scuso per la lunga spiegazione. Per fare ciò in JMP Pro 18 (JMP Student Edition), è necessario eseguire la programmazione indicata sopra. Il linguaggio di scripting JMP consente anche di creare applicazioni con interfaccia grafica. Alcune funzionalità non disponibili in JMP sono disponibili su JMP Marketplace (le estensioni offerte sono fornite così come sono). Purtroppo, la funzionalità di regressione di Poisson modificata non è disponibile su JMP Marketplace. Se potete aspettare, vi consiglio di attendere il rilascio di JMP Pro 19.</DIV>
<DIV>
</DIV>
<DIV>La possibilità o meno di eseguire un determinato processo in JMP può essere suddivisa nei seguenti livelli. In questo caso, credo che corrisponda al livello d).</DIV>
<DIV>a) Funzioni comunemente disponibili in JMP e JMP Pro → Possono essere eseguite in pochissimo tempo</DIV>
<DIV>
<DIV>b) Se si utilizzano le funzioni di JMP o JMP Pro → È possibile farlo rapidamente con un po' di ingegno</DIV>
</DIV>
<DIV>
<DIV>c) Funzionalità fornite su JMP Marketplace → Così come sono, ma pronte per l'uso</DIV>
<DIV>d) Una funzione semplice che non è fornita sopra → Può essere realizzata tramite programmazione (specialmente in R o Python)</DIV>
<DIV>e) Funzioni altamente avanzate non fornite sopra → Possibili solo con notevole competenza</DIV>
<DIV>
</DIV>
<DIV>Infine, c'è un ultimo punto che mi ha infastidito, quindi vorrei aggiungere un commento sulla "correttezza" che ho menzionato all'inizio. L'analisi di regressione in sé, inclusa la più semplice regressione lineare con i minimi quadrati, ha una varietà di utilizzi. Ad esempio, il seguente articolo la suddivide in descrizione, previsione e controllo/intervento.</DIV>
<DIV>Iwasaki, M. (2021)</DIV>
<DIV>Analisi di regressione multipla dal punto di vista dell'inferenza causale statistica</DIV>
<DIV>Rivista della Società di Statistica Giapponese, 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 spiegazione di cui sopra descrive semplicemente come eseguire la regressione di Poisson modificata e non stabilisce in alcun modo la "correttezza" nel senso che sia "corretto" eseguire l'analisi di regressione nella situazione attuale, se eseguire la regressione di Poisson modificata sia "corretto" o anche se la situazione attuale sia "corretta" nel senso che soddisfa le principali ipotesi formulate dall'analisi di regressione. Se la vostra attuale analisi statistica sia "corretta" per la vostra situazione attuale (dato che io stesso non sono in grado di giudicare) è qualcosa che consiglierei di consultare uno statistico applicato in quel campo (non ne ho mai consultato uno, ma credo che possa essere in grado di suggerire alternative all'analisi stessa e anche di fornire istruzioni dettagliate su come eseguire la regressione di Poisson modificata).</DIV>
<DIV>
</DIV>
</DIV>
