Estimation des ratios de risque ajustés dans les données binaires à l'aide de la régression de Poisson

9

Je m'intéresse à l'estimation d'un rapport de risque ajusté, analogue à la façon dont on estime un rapport de cotes ajusté en utilisant la régression logistique. Certaines publications (par exemple, ceci ) indiquent que l'utilisation de la régression de Poisson avec les erreurs standard de Huber-White est une façon basée sur un modèle de le faire

Je n'ai pas trouvé de littérature sur la façon dont l'ajustement pour les covariables continues affecte cela. La simulation simple suivante montre que ce problème n'est pas aussi simple:

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 

Dans ce cas, le rapport de risque réel est 2, qui est récupéré de manière fiable lorsque l'effet covariable est faible. Mais, lorsque l'effet covariable est important, cela est faussé. Je suppose que cela se produit parce que l'effet de covariable peut repousser la limite supérieure (1) et cela contamine l'estimation.

J'ai regardé mais je n'ai trouvé aucune littérature sur l'ajustement pour les covariables continues dans l'estimation du rapport de risque ajusté. Je connais les messages suivants sur ce site:

mais ils ne répondent pas à ma question. Y a-t-il des articles à ce sujet? Y a-t-il des mises en garde connues à appliquer?

kjetil b halvorsen
la source
1
Peut vous intéresser: aje.oxfordjournals.org/content/162/3/199.full
StatsStudent
De plus, cette Q&A stats.stackexchange.com/questions/18595/… peut aider.
mdewey

Réponses:

1

Je ne sais pas si vous avez toujours besoin d'une réponse à cette question, mais j'ai un problème similaire dans lequel j'aimerais utiliser la régression de Poisson. En exécutant votre code, j'ai constaté que si je configurais le modèle comme

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

plutôt que comme votre modèle de régression de Poisson, le même résultat se produit: le OR estimé est ~ 1,5 à l'approche de ce 1. Donc, je ne suis pas sûr que votre exemple fournisse des informations sur un problème possible avec l'utilisation de la régression de Poisson pour les résultats binaires.

David F
la source
1
Le problème lié à l'ajustement d'un modèle logit, bien qu'il n'entraîne pas de risques prévus supérieurs à 1, est que le rapport de cotes est un estimateur biaisé du rapport de risque et que le biais augmente considérablement à mesure que le résultat devient plus répandu. Vous pouvez spécifier binomial(link=log)de s'adapter à un modèle de risque relatif, mais il converge rarement en raison d'une sur-prédiction des résultats.
AdamO
1

Je trouve que l'utilisation du maximum de vraisemblance directe avec la fonction de probabilité appropriée améliore considérablement l'estimation du risque relatif. Vous pouvez spécifier directement la fonction de risque tronquée comme taux prévu pour le processus.

entrez la description de l'image ici

Habituellement, nous utilisons la Hesse pour créer des IC pour l'estimation. Je n'ai pas exploré la possibilité de l'utiliser comme matrice "B" (viande) dans l'erreur Huber White et d'utiliser les risques ajustés pour obtenir la matrice "A" (pain) ... mais je soupçonne que cela pourrait fonctionner! De manière plus pratique, vous pouvez utiliser un bootstrap pour obtenir des erreurs de modèle qui sont robustes à une relation moyenne-variance mal spécifiée.

## 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))

Donne:

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

Le coefficient moyen vous donne ce que vous voulez.

AdamO
la source