Utilisation de l'offset dans le modèle binomial pour tenir compte du nombre accru de patients

18

Deux questions connexes de moi. J'ai un cadre de données qui contient le nombre de patients dans une colonne (plage de 10 à 17 patients) et 0 et 1 indiquant si un incident s'est produit ce jour-là. J'utilise un modèle binomial pour régresser la probabilité d'incident sur le nombre de patients. Cependant, je voudrais tenir compte du fait que lorsqu'il y a plus de patients, il y aura inévitablement plus d'incidents parce que le temps total passé par les patients dans la salle est plus élevé ce jour-là.

J'utilise donc un modèle binomial décalé comme celui-ci (code R):

glm(Incident~Numbers, offset=Numbers, family=binomial, data=threatdata)

Mes questions sont:

  1. Est-il correct d'avoir exactement les mêmes variables de prédiction et de décalage? Je veux compenser l'augmentation tonique de la probabilité d'incident et voir s'il reste quelque chose, essentiellement. Cela a du sens pour moi, mais je suis un peu prudent au cas où je me trompe.

  2. Le décalage est-il spécifié correctement? Je sais que dans les modèles de poisson, il se lirait

    offset=log(Numbers)
    

Je ne sais pas s'il y a un équivalent ici et je n'arrive pas à trouver de compensations binomiales avec Google (le problème majeur étant que je reçois toujours des binômes négatifs, ce qui bien sûr n'est pas bon).

Chris Beeley
la source
2
N'est-ce pas ce que vous cherchez à ajuster, précisément ce que vous voulez mesurer - c'est-à-dire comment la probabilité «d'incident» augmente avec le nombre de patients?
B_Miner
1
Je dois faire écho au point de B_Miner. Je pense que vous êtes confus quand / pourquoi le décalage est utilisé dans cette situation. Votre modèle, sans décalage, vous donnera des valeurs ajustées de probabilité d'incident en fonction du nombre de patients. Si vous êtes intéressé par une forme fonctionnelle différente, envisagez des transformations (comme le journal ou l'exponentiation de #) en fonction de ce qui est scientifiquement intéressant.
AdamO
Pouvez-vous clarifier quelque chose sur les incidents? Un incident est-il lié à un patient ou quelque chose au sujet de la salle dans son ensemble? S'ils sont liés à des patients, est-il possible qu'il y ait> 1 incident S'il n'y a pas de patients, est-il impossible d'avoir un incident?
atiretoo
1
Apparemment, ma réponse «ne contient pas suffisamment de détails». J'ai fourni un développement théorique, du code exécutable et des réponses à vos deux questions, alors pouvez-vous peut-être clarifier ce qui est nécessaire de plus?
conjugateprior
1
Désolé, Conjugué Prior, votre réponse est excellente. La chose "pas assez de détails" était la balise ajoutée à la prime (c'est-à-dire qu'elle était là avant que vous ne postiez). Je vais accepter la fin de la prime au cas où quelqu'un produirait une réponse encore meilleure, mais cela est peu probable et la vôtre est très utile, merci.
Chris Beeley

Réponses:

17

Si vous êtes intéressé par la probabilité d'un incident avec N jours de patients en service, vous souhaitez un modèle comme:

mod1 <- glm(incident ~ 1, offset=patients.on.ward, family=binomial)

le décalage représente les essais, incidentest soit 0 soit 1, et la probabilité d'un incident est constante (pas d'hétérogénéité dans la tendance à générer des incidents) et les patients n'interagissent pas pour provoquer des incidents (pas de contagion). Alternativement, si la probabilité d'un incident est faible, ce qui est pour vous (ou si vous avez limité le nombre d'incidents sans nous le mentionner), vous préférerez peut-être la formulation de Poisson

log.patients.on.ward <- log(patients.on.ward)
mod2 <- glm(incident ~ 1, offset=log.patients.on.ward, family=poisson)

où les mêmes hypothèses s'appliquent. Le décalage est enregistré car le nombre de patients en salle a un effet proportionnel / multiplicateur.

En développant le deuxième modèle, vous pensez peut-être qu'il y a plus d' incidents que ce qui serait prévu autrement simplement en raison de l'augmentation du nombre de patients. Autrement dit, peut-être que les patients interagissent ou sont hétérogènes. Alors tu essaies

mod3 <- glm(incident ~ 1 + log.patients.on.ward, family=poisson)

Si le coefficient sur log.patients.on.wardest significativement différent de 1, où il a été fixé mod2, alors quelque chose peut en effet être faux avec vos hypothèses d'absence d'hétérogénéité et de contagion. Et bien que vous ne puissiez bien sûr pas distinguer ces deux (ni l'une ni l'autre des autres variables manquantes), vous avez maintenant une estimation de l'augmentation du nombre de patients en salle augmente le taux / la probabilité d'un incident au-delà de ce que vous auriez attendre du hasard. Dans l'espace des paramètres c'est 1-coef(mod3)[2]avec un intervalle dérivable de confint.

Alternativement, vous pouvez simplement travailler avec la quantité de journal et son coefficient directement. Si vous voulez simplement prédire la probabilité d'incident en utilisant le nombre de patients en salle, alors ce modèle serait un moyen simple de le faire.

Questions

  1. Est-il correct d'avoir des variables dépendantes dans votre décalage? Cela me semble être une très mauvaise idée, mais je ne vois pas que vous le deviez.

  2. Le décalage dans les modèles de régression de Poisson pour exposureest en effet log(exposure). Il est peut-être déroutant que l'utilisation des offsetmodèles de régression binomiale de R soit essentiellement un moyen d'indiquer le nombre d'essais. Il peut toujours être remplacé par une variable dépendante définie comme cbind(incidents, patients.on.ward-incidents)et sans décalage. Considérez-le comme ceci: dans le modèle Poisson, il entre sur le côté droit derrière la fonction de liaison logarithmique, et dans le modèle binomial, il entre sur le côté gauche devant la fonction de liaison logitique.

conjugateprior
la source
18

Offsets dans les régressions de Poisson

Commençons par examiner pourquoi nous utilisons un décalage dans une régression de Poisson. Nous voulons souvent en raison de cela pour contrôler l'exposition. Soit le taux de base par unité d'exposition et t le temps d'exposition dans les mêmes unités. Le nombre prévu d'événements sera λ × t .λtλ×t

Dans un modèle GLM, nous modélisons la valeur attendue à l'aide d'une fonction de lien , c'est-à-direg

g(λti)=log(λti)=β0+β1x1,i+

tiixii

Nous pouvons simplifier simplifier l'expression ci-dessus

log(λ)=log(ti)+β0+β1x1,i+

log(ti)

Régression binomiale

Dans une régression binomiale, qui utilise généralement un lien logit, c'est:

g(pi)=logit(pi)=log(pi1pi)=β0+β1x1,i+

pi

pii

pi,jjiij=1Ni(1pi,j)Nii

pi=1j=1Ni(1pi,j).

pi=1(q)Ni,
q=1pp

pig(pi)log((q)N1)

Par conséquent, nous ne pouvons pas utiliser de décalage dans ce cas.

p

Rider_X
la source
2
+1, bienvenue sur le site, @Rider_X. J'espère que nous pouvons espérer davantage de réponses de ce type à l'avenir.
gung - Rétablir Monica
1
@gung - Merci! Je n'avais pas beaucoup entendu parler de ce que je pensais être une réponse utile, donc je ne suis pas revenu beaucoup. Je vais devoir changer cela. Cordialement.
Rider_X
2
+1 J'apprécie vraiment les réponses qui expliquent la théorie et le raisonnement, plutôt que (ou en plus de) montrer quel code et quelles commandes utiliser.
whuber
9

Cette réponse se décompose en deux parties, la première une réponse directe à la question et la seconde un commentaire sur le modèle que vous proposez.

La première partie se rapporte à l'utilisation de Numberscomme décalage en plus de l'avoir sur le côté droit de l'équation. L'effet de cette opération sera simplement de soustraire 1 du coefficient estimé de Numbers, inversant ainsi l'effet du décalage, et ne changera pas autrement les résultats. L'exemple suivant, avec quelques lignes de sortie non pertinentes supprimées, le démontre:

library(MASS)
Numbers <- rpois(100,12)
p <- 1 / (1 + exp(0.25*Numbers))
y <- rbinom(100, Numbers, p)
Incident <- pmin(y, 1) 

> summary(glm(Incident~Numbers, family="binomial"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3121  -1.0246  -0.8731   1.2512   1.7465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  0.99299    0.80624   1.232   0.2181  
Numbers     -0.11364    0.06585  -1.726   0.0844 . <= COEFFICIENT WITH NO OFFSET TERM
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 135.37  on 99  degrees of freedom
Residual deviance: 132.24  on 98  degrees of freedom
AIC: 136.24

> summary(glm(Incident~Numbers, offset=Numbers, family="binomial"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.3121  -1.0246  -0.8731   1.2512   1.7465  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.99299    0.80624   1.232    0.218    
Numbers     -1.11364    0.06585 -16.911   <2e-16 *** <= COEFFICIENT WITH OFFSET TERM
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 342.48  on 99  degrees of freedom
Residual deviance: 132.24  on 98  degrees of freedom
AIC: 136.24

Notez que tout est le même sauf le coefficient de nombres et la déviance nulle (et la statistique t, car il teste toujours contre 0 au lieu de -1.)

La deuxième partie concerne le modèle que vous construisez. Étant donné que les incidents ne sont pas enregistrés comme le nombre d'incidents dans une journée, mais eu des incidents dans une journée, la probabilité d'observer un 1 le jourt est 1-(1-pt)Nt, où Nt est le nombre de patients le jour t et pt est la probabilité par patient d'un incident le jour t. La fonction de lien habituelle, le logit, paramètrerait ceci commeJournal(1-(1-pt)Nt)/NtJournal(1-pt). Cela indique que la relation entre la probabilité d’observer un 1 le jourt et Ntpeut ne pas être bien modélisé par une fonction linéaire sur l'échelle logit. (Cela peut être le cas de toute façon, car on pourrait s'attendre à un "seuil" approximatif en dessous duquel la qualité des soins aux patients est OK mais au-dessus duquel la qualité des soins aux patients chute rapidement.) Inverser la définition des probabilités afin de déplacer laNt dans le dénominateur au lieu du numérateur vous laisse toujours cette exponentielle maladroite à l'intérieur du journal.

On pourrait également soupçonner que la probabilité par patient varie d'un patient à l'autre, ce qui conduirait à un modèle hiérarchique plus complexe, mais je n'entrerai pas dans les détails ici.

Dans tous les cas, étant donné cela et la plage limitée du nombre de patients que vous observez, plutôt que d'utiliser un modèle linéaire sur l'échelle logit, il pourrait être préférable de ne pas être paramétrique sur la relation et de grouper le nombre de patients en trois ou quatre groupes, par exemple, 10-11, 12-13, 14-15 et 16-17, construisent des variables muettes pour ces groupes, puis exécutent la régression logistique avec les variables muettes sur le côté droit. Cela permettra de mieux saisir les relations non linéaires telles que "le système est surchargé autour de 16 patients et les incidents commencent à augmenter de manière significative." Si vous aviez un éventail de patients beaucoup plus large, je suggérerais un modèle additif généralisé, par exemple, «gam» du package «mgcv».

jbowman
la source
0

Il semble plus simple de spécifier un lien de journal et de conserver le décalage comme pour un modèle de Poisson.

un arrêt
la source
2
Je suis sûr que vous avez raison, mais pour moi, comment est-ce un poisson? On dirait que l'OP a un ensemble de données avec un résultat binaire. Serait-ce glm (Incident ~ Numbers, offset = log (Numbers), family = poisson, data = menacedonnées) ??
B_Miner