Können Gewichte und Versatz zu ähnlichen Ergebnissen bei der Poisson-Regression führen?

In "A Practioner's Guide to Generalized Linear Models" in Absatz 1.83 heißt es:

"Im speziellen Fall eines multiplikativen Poisson-GLM kann gezeigt werden, dass die Modellierungsanspruchszahlen mit einem Versatzterm gleich dem Logarithmus der Exposition zu identischen Ergebnissen führen wie die Modellierungsanspruchshäufigkeiten mit vorherigen Gewichten, die gleich der Exposition jeder Beobachtung eingestellt sind. ""

Ich kann keine weiteren Referenzen zu diesen Ergebnissen finden, daher habe ich einige empirische Tests durchgeführt, bei denen ich keinen Beweis dafür finden konnte, dass die Aussage korrekt ist. Kann jemand einen Einblick geben, warum diese Ergebnisse richtig / falsch sein können.

Zu Ihrer Information, ich habe den folgenden R-Code verwendet, um die Hypothese zu testen, bei der ich für die beiden genannten Fälle keine ähnlichen Ergebnisse erzielen konnte:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y'=y, 'X'=X, 'offset' = offset)
formula = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula, family = "poisson", df, weights = offset) 
#Third model using poisson model on y/offset as reference
dfNew = df
dfNew$y = dfNew$y/offset
fit3 = glm(formula, family = "poisson", dfNew)

#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, fit3$coefficients, c(intercept,coefs))

Die Koeffizientenschätzungen, die sich aus der Ausführung dieses Codes ergeben, sind nachstehend aufgeführt:

 >  
    (Intercept)       X.1       X.2       X.3        X.4       X.5       X.6
[1,]    1.998277 0.2923091 0.4586666 0.1802960 0.11688860 0.7997154 0.4786655
[2,]    1.588620 0.2708272 0.4540180 0.1901753 0.07284985 0.7928951 0.5100480
[3,]    1.983903 0.2942196 0.4593369 0.1782187 0.11846876 0.8018315 0.4807802
[4,]    2.000000 0.2909240 0.4576965 0.1807591 0.11658183 0.8005451 0.4780123
              X.7       X.8       X.9      X.10
[1,]  0.005772078 0.9154808 0.9078758 0.3512824
[2,] -0.003705015 0.9117014 0.9063845 0.4155601
[3,]  0.007595660 0.9181014 0.9076908 0.3505173
[4,]  0.005881960 0.9150350 0.9084375 0.3511749
>

und wir können beobachten, dass die Koeffizienten nicht identisch sind.

generalized-linear-model modeling poisson-regression weights offset BDP1
quelle

Sie sollten nicht wirklich rm(list=ls() )in den R-Code aufnehmen, den Sie hier posten! Das könnte jemanden überraschen, der es betreibt, und ihn wütend auf dich machen. Ich habe es entfernt. Ich habe auch bearbeitet, um die Ergebnisse der Ausführung des Codes aufzunehmen. Wenn Sie dies ursprünglich getan hätten, hätten Sie möglicherweise eine schnellere Antwort erhalten, da nur wenige Leser den Code selbst ausführen.

kjetil b halvorsen

@kjetilbhalvorsen Danke für das Feedback, wird in Zukunft tun.

BDP1

@ BDP1 Derzeit testet Ihr R-Code nicht die Behauptung, auf die Sie sich beziehen. Ich schlage vor, dass Sie den Gewichtsbegriff für fit3 hinzufügen und der Frage direkt einen Nachtrag hinzufügen.

Aivanov

Antworten:

(Mit Ihrem R-Code können Sie "poisson" durch "quasipoisson" ersetzen, um alle generierten Warnungen zu vermeiden. Ansonsten ändert sich nichts anderes. Siehe (*) unten). Ihre Referenz verwendet den Begriff "multiplikatives glm", der meiner Meinung nach nur ein glm mit Protokollverknüpfung bedeutet, da ein Protokollverknüpfung als multiplikatives Modell betrachtet werden kann. Ihr eigenes Beispiel zeigt, dass die Behauptung falsch ist, zumindest so, wie wir sie interpretiert haben (da die geschätzten Parameter nicht gleich sind). Sie könnten den Autoren schreiben und sie fragen, was sie bedeuten. Im Folgenden werde ich argumentieren, warum die Behauptung falsch ist.

Sei der Poisson-Parameter und die Gewichte. Sei der lineare Prädiktor ohne Offset und dann der lineare Prädiktor mit dem Offset. Die Poisson-Wahrscheinlichkeitsfunktion ist Dann wird die Log-Likelihood-Funktion für das Modell mit Offset während die Log-Likelihood-Funktion für das Modell mit Gewichten $\lambda_i$ $\omega_i$ $\eta_i$ $\eta_i+\log(\omega_i)$

f (y_{i}) = e^{- λ_{i}} λ_{i}^{y_{i}} / y_{i}!

$f(y_i) = e^{-\lambda_i} \lambda_i^{y_i}/y_i !$

ℓ = - \sum_{i} ω_{i} e^{η_{i}} + \sum_{i} y_{i} η_{i} + \sum_{i} y_{i} \log ω_{i} - \sum_{i} \log y_{i}!

$\ell = -\sum_i \omega_i e^{\eta_i} + \sum_i y_i \eta_i +\sum_i y_i\log \omega_i - \sum_i \log y_i!$

ℓ^{w} = - \sum_{i} ω_{i} e^{η_{i}} + \sum_{i} y_{i} ω_{i} η_{i} - \sum_{i} ω_{i} \log y_{i}!

$\ell^w = -\sum_i \omega_i e^{\eta_i}+\sum_i y_i \omega_i \eta_i -\sum_i \omega_i \log y_i!$ und das ist eindeutig nicht dasselbe. Was diese Autoren meinten, ist mir also nicht klar.

(*) Anmerkung aus der Hilfe der glmFunktion von R :

Nicht-NULL-Gewichte können verwendet werden, um anzuzeigen, dass unterschiedliche Beobachtungen unterschiedliche Dispersionen aufweisen (wobei die Werte in Gewichten umgekehrt proportional zu den Dispersionen sind). oder äquivalent, wenn die Elemente von 'Gewichten' positive ganze Zahlen w_i sind, ist jede Antwort y_i der Mittelwert von w_i Beobachtungen des Einheitsgewichts. Für ein binomiales GLM werden vorherige Gewichte verwendet, um die Anzahl der Versuche anzugeben, wenn die Antwort der Anteil der Erfolge ist: Sie würden selten für ein Poisson-GLM verwendet.

Ein Blick auf die Bedeutung der Gewichtungsargumente erklärt dies und gibt der Funktion der Poisson-Familie, die einen konstanten Skalierungsparameter annimmt, wenig Bedeutung, während die Gewichtungsargumente modifizieren . Dies gibt der Funktion der Quasiposson-Familie mehr Bedeutung. Siehe die Antwort auf die Eingabe "Gewicht" in den Funktionen glm und lm in R. Die dort gegebene Antwort hilft auch zu erkennen, warum die Wahrscheinlichkeit im gewichteten Fall die oben angegebene Form annimmt. $\phi=1$ $\phi$

Die hier gegebene Antwort könnte relevant sein: Wie ist eine Poisson-Ratenregression gleich einer Poisson-Regression mit entsprechendem Offset-Term? und ist sehr interessant.

kjetil b halvorsen
quelle

Danke für die Antwort. Das Durchlaufen der Wahrscheinlichkeitsbeiträge verdeutlicht viel. Bei der Suche nach weiteren Antworten auf Ihre Antwort habe ich die folgende Frage gefunden: stats.stackexchange.com/questions/66791/… In der gezeigt wird, dass durch Teilen der ursprünglichen Antwort durch den Versatz und Festlegen der Versatzvariablen als Gewichte Die gleichen Ergebnisse können wie beim Basismodell erzielt werden, bei dem der Logarithmus (Offset) mit einem konstanten Koeffizienten von 1 in den linearen Prädiktor eintritt.

BDP1

Ich habe auch versucht abzuleiten, dass die Wahrscheinlichkeitsbeiträge dieses neuen Modells den Beiträgen des neu erwähnten Modells entsprechen, konnte dies jedoch angesichts vonBegriff, der für letztere erscheint.

(y_{i} / w_{i})!

$(y_i/w_i)!$

BDP1

Zumindest im Experiment im bereitgestellten R-Skript scheint diese Behauptung wahr zu sein. Dies kann leicht durch Hinzufügen eines Gewichts = Offset-Arguments in der Zeile fit3 = glm (...)

BDP1

Ich verstehe nicht, was Sie hier sagen. Wenn Sie der Meinung sind, dass Ihr Experiment dies bestätigt, warum bearbeiten Sie diese Änderung nicht, um das OP zu erläutern, und erklären Sie, warum Sie der Meinung sind, dass dies die Behauptung bestätigt. Ich habe versucht zu erklären, warum es nicht stimmt, was ist falsch an meiner Argumentation?

kjetil b halvorsen

Das Zitat aus dem Practioner's Guide ist tatsächlich korrekt, es wurde von OP einfach nicht korrekt implementiert. Dies wurde von Alan Chalk in einer anderen Antwort hier hervorgehoben, von mir stats.stackexchange.com/questions/430001 und auch von Cokes stats.stackexchange.com/questions/264071 (obwohl die richtige Schlussfolgerung in der letzten Zeile von begraben ist diese Antwort).

Gordon Smyth

Es tut mir leid, dass ich oben nicht einfach einen Kommentar hinzugefügt habe, aber ich habe nicht genug Wiederholungen.

Die ursprüngliche Behauptung - aber ein wenig modifiziert - ist tatsächlich richtig.

Die folgenden zwei Modelle geben in R mit einem Poisson glm mit Log-Link genau die gleiche Antwort:

Modell y, verwenden Sie einen Versatz als Protokoll (Versatz)
Modell y / Versatz, verwenden Sie Gewichte gleich Versatz

Wenn Sie Ihren Originalcode leicht anpassen, werden identische Antworten angezeigt:

n=1000
m=10

# Generate random data
X = matrix(data = rnorm(n*m)+1, ncol = m, nrow = n)

intercept = 2
coefs = runif(m)
offset = runif(n)
## DGP: exp of Intercept + linear combination X variables + log(offset)
mu = exp(intercept + X%*%coefs + log(offset))
y = rpois(n=n, lambda=mu)

df = data.frame('y' = y,
                'y_over_offset' = y/offset,
                'X' = X,
                'offset' = offset)

formula_offset = paste("y ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))
formula_weights = paste("y_over_offset ~",paste(colnames(df)[grepl("X", colnames(df))], collapse = "+"))

#First model using log(offset) as offset
fit1  = glm(formula_offset, family = "poisson", df, offset = log(offset))
#Second model using offset as weights for individual observations
fit2 = glm(formula_weights, family = "poisson", df, weights = offset) 


#Combine coefficients with the true coefficients
rbind(fit1$coefficients, fit2$coefficients, c(intercept,coefs))

Hoffentlich sollte das identische Antworten geben.

Es ist möglich zu zeigen, dass die beiden Modelle statistisch äquivalent sind (irgendwo gibt es ein CAS-Papier, das dies zeigt - ich werde einen Link posten, wenn ich Zeit habe).

Wenn Sie eine bestrafte Regression durchführen, kann die Art und Weise, wie verschiedene Pakete wie glmnet und H2o die Abweichung für die beiden verschiedenen Arten der Definition eines Modells messen, zu unterschiedlichen Ergebnissen führen.

Alan Chalk
quelle

Bei einer kurzen Frage können Sie sehen, dass Sie mit fit2 keinen AIC haben und die Diagramme zwischen den beiden Anpassungen leicht unterschiedlich sind (fit1). Handlung (fit2) - wissen Sie auch, warum dies passiert?]

Charl Francois Marais