Schätzung der angepassten Risikoverhältnisse in binären Daten unter Verwendung der Poisson-Regression

9

Ich bin daran interessiert, ein angepasstes Risikoverhältnis zu schätzen, analog dazu, wie man ein angepasstes Quotenverhältnis unter Verwendung der logistischen Regression schätzt. Einige Literaturstellen (z. B. diese ) weisen darauf hin, dass die Verwendung der Poisson-Regression mit Huber-White-Standardfehlern eine modellbasierte Methode ist, um dies zu erreichen

Ich habe keine Literatur darüber gefunden, wie sich die Anpassung an kontinuierliche Kovariaten darauf auswirkt. Die folgende einfache Simulation zeigt, dass dieses Problem nicht so einfach ist:

arr <- function(BLR,RR,p,n,nr,ce)
{
   B = rep(0,nr)
   for(i in 1:nr){
   b <- runif(n)<p 
   x <- rnorm(n)
   pr <- exp( log(BLR) + log(RR)*b + ce*x)
   y <- runif(n)<pr
   model <- glm(y ~ b + x, family=poisson)
   B[i] <- coef(model)[2]
   }
   return( mean( exp(B), na.rm=TRUE )  )
}

set.seed(1234)
arr(.3, 2, .5, 200, 100, 0)
[1] 1.992103
arr(.3, 2, .5, 200, 100, .1)
[1] 1.980366
arr(.3, 2, .5, 200, 100, 1)
[1] 1.566326

In diesem Fall beträgt das wahre Risikoverhältnis 2, das zuverlässig wiederhergestellt wird, wenn der kovariate Effekt gering ist. Wenn der kovariate Effekt jedoch groß ist, wird dies verzerrt. Ich gehe davon aus, dass dies entsteht, weil der kovariate Effekt gegen die Obergrenze (1) drücken kann und dies die Schätzung kontaminiert.

Ich habe nachgesehen, aber keine Literatur zur Anpassung an kontinuierliche Kovariaten bei der Schätzung des angepassten Risikoverhältnisses gefunden. Mir sind folgende Beiträge auf dieser Seite bekannt:

aber sie beantworten meine Frage nicht. Gibt es irgendwelche Papiere dazu? Gibt es bekannte Vorsichtsmaßnahmen, die ausgeübt werden sollten?

r references biostatistics poisson-regression relative-risk kjetil b halvorsen
quelle
1
Kann
StatsStudent
Auch diese Fragen und Antworten unter stats.stackexchange.com/questions/18595/… können hilfreich sein.
Mdewey

Antworten:

1

Ich weiß nicht, ob Sie noch eine Antwort auf diese Frage benötigen, aber ich habe ein ähnliches Problem, bei dem ich die Poisson-Regression verwenden möchte. Beim Ausführen Ihres Codes habe ich festgestellt, dass, wenn ich das Modell als eingerichtet habe

model <- glm(y ~ b + x, family=binomial(logit)

Anstatt wie bei Ihrem Poisson-Regressionsmodell tritt dasselbe Ergebnis auf: Der geschätzte OR beträgt ~ 1,5, wenn sich ce 1 nähert. Ich bin mir also nicht sicher, ob Ihr Beispiel Informationen zu einem möglichen Problem bei der Verwendung der Poisson-Regression für binäre Ergebnisse enthält.

David F.
quelle
1
Das Problem bei der Anpassung eines Logit-Modells, obwohl es nicht zu vorhergesagten Risiken von mehr als 1 führt, besteht darin, dass das Odds Ratio ein voreingenommener Schätzer des Risikoverhältnisses ist und dass die Voreingenommenheit dramatisch zunimmt, wenn das Ergebnis vorherrscht. Sie können angeben binomial(link=log), ob ein relatives Risikomodell tatsächlich angepasst werden soll, es konvergiert jedoch selten, da das Ergebnis zu stark vorhergesagt wird.
AdamO
1

Ich finde, dass die Verwendung der direkten maximalen Wahrscheinlichkeit mit der richtigen Wahrscheinlichkeitsfunktion die Schätzung des relativen Risikos erheblich verbessert. Sie können die abgeschnittene Risikofunktion direkt als die vorhergesagte Rate für den Prozess angeben.

Geben Sie hier die Bildbeschreibung ein

Normalerweise verwenden wir das Hessische, um CIs für die Schätzung zu erstellen. Ich habe nicht die Möglichkeit untersucht, dies als "B" -Matrix (Fleisch) im Huber White-Fehler zu verwenden und die angepassten Risiken zu verwenden, um die "A" -Matrix (Brot) zu erhalten ... aber ich vermute, dass es funktionieren könnte! Es ist praktikabler, einen Bootstrap zu verwenden, um Modellfehler zu erhalten, die für eine falsch spezifizierte Mittelwert-Varianz-Beziehung robust sind.

## the negative log likelihood for truncated risk function
negLogLik <- function(best, X, y) { 
  pest <- pmin(1, exp(X %*% best))
  -sum(dpois(x = y, lambda = pest, log=TRUE))
}

set.seed(100)

sim <- replicate(100, {
  n <- 200
  X <- cbind(1, 'b'=rbinom(n, 1, 0.5), 'x'=rnorm(n))
  btrue <- c(log(0.3), log(2), 1)
  ptrue <- pmin(1, exp(X %*% matrix(btrue)))
  y <- rbinom(n, 1, ptrue) ## or just take y=ptrue for immediate results
  nlm(f = logLik, p = c(log(mean(y)),0,0), X=X, y=y)$estimate
})

rowMeans(exp(sim))

Gibt:

> rowMeans(exp(sim))
[1] 0.3002813 2.0680780 3.0888280

Der mittlere Koeffizient gibt Ihnen, was Sie wollen.

AdamO
quelle