Comment puis-je estimer le moment auquel 50% d'une variable binomiale aura transité?

8

J'ai les données suivantes, représentant l'état binaire de quatre sujets à quatre reprises, notez qu'il n'est possible que pour chaque sujet de passer de mais pas de :0110

testdata <- data.frame(id = c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4),
                       day = c(1,1,1,1,8,8,8,8,16,16,16,16,24,24,24,24,32,32,32,32),
                       obs = c(0,0,0,0,0,1,0,0,0,1,1,0,0,1,1,1,1,1,1,1))

Je peux le modéliser avec une régression logistique:

testmodel <- glm(formula(obs~day, family=binomial), data=testdata)

> summary(testmodel)


Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.018890   0.148077  -0.128 0.899907    
day          0.032030   0.007555   4.240 0.000493 ***

Premièrement, comment puis-je tenir compte des mesures répétées sur le même individu dans le modèle?

Deuxièmement, comment estimer, avec incertitude, le jour où la moitié des sujets auront fait la transition de ?01

David LeBauer
la source
1
Il semble qu'il y ait une forte dépendance dans ces données: à savoir, est-ce le cas que si obs = 1 pour le sujet le jour alors nécessairement obs = 1 pour le sujet le jour chaque fois que ? Si tel est le cas, vous n'avez vraiment que quatre valeurs de données - une pour chaque sujet - et l'une d'entre elles est censurée à droite. itisst
whuber
@quand vous avez raison sur la dépendance (au moins dans la présente analyse intra-annuelle); les données indiquent si un «débourrement» s'est produit ou non avant la date d'observation pour chacun des quatre arbres répliqués. Mais je ne sais pas ce que vous voulez dire sur les valeurs de données censurées à droite?
David LeBauer
1
Voici un résumé: le sujet 2 est passé dans l'intervalle [1,8]; c'est-à-dire 2 -> [1,8]. Aussi 3 -> [8,16], 4 -> [16,24] et 1 -> [24, infini]. Cette dernière signifie que le sujet 1 a été observé pendant 24 jours sans transition; c'est la valeur censurée. Vous pouvez définir cela comme un problème d'analyse de survie et l'analyser en conséquence. Soit dit en passant, cette dépendance signifie que les valeurs de p dans la régression logistique sont trompeusement faibles.
whuber
@whuber merci pour la perspicacité, mais cela signifie-t-il que mon approche est fondamentalement défectueuse étant donné que je ne suis pas intéressé par l'estimation des valeurs de p? De plus, aucune des données ne sera censurée à droite dans quelques semaines; Je suis en train de développer l'analyse avant que l'ensemble de données ne soit terminé. J'ai modifié les données du test afin qu'aucun des sujets ne soit correctement censuré.
David LeBauer
3
@DWin, @David Ce n'est pas une situation de mesure répétée. Le format des données ne fait que ressembler à ça. La mesure pour chaque sujet consiste en un seul intervalle pendant lequel une transition a été observée.
whuber

Réponses:

3

Comme il est devenu évident dans les commentaires à la question, les données ne comprennent que quatre observations du délai d'éclosion des bourgeons. (Ce serait une erreur de les analyser comme s'il s'agissait de 16 valeurs indépendantes.) Elles consistent en des intervalles de temps plutôt qu'en des temps exacts:

[1,8], [8,16], [16,24], [24,32]

Il existe plusieurs approches. Un appel, très général, est de prendre ces intervalles au mot: le vrai moment du débourrement pourrait être n'importe quoi dans chaque intervalle. Nous sommes donc amenés à représenter «l'incertitude» sous deux formes distinctes: l'incertitude d'échantillonnage (nous avons un échantillon vraisemblablement représentatif de l'espèce cette année) et l' incertitude d'observation (reflétée par les intervalles).

L'incertitude d'échantillonnage est gérée avec des techniques statistiques connues: on nous demande d'estimer la médiane et nous pouvons le faire de plusieurs façons, en fonction d'hypothèses statistiques, et nous pouvons fournir des intervalles de confiance pour l'estimation. Pour simplifier, supposons que le temps de débourrement ait une distribution symétrique. Parce qu'il est (vraisemblablement) non négatif, cela implique qu'il a une variance et suggère également que la moyenne de seulement quatre observations peut être distribuée approximativement normalement. De plus, la symétrie implique que nous pouvons utiliser la moyenne comme substitut de la médiane (ce qui est recherché dans la question d'origine). Cela nous donne accès à des méthodes standard, simples, d'estimation et d'intervalle de confiance.

L'incertitude d'observation peut être gérée avec les principes de l'arithmétique d'intervalle (souvent appelés «analyse des limites de probabilité» ): effectuer tous les calculs en utilisant toutes les configurations possibles de données cohérentes avec les observations. Voyons comment cela fonctionne dans un cas simple: l'estimation de la moyenne. Il est intuitivement clair que la moyenne ne peut être inférieure à = , obtenue en utilisant les plus petites valeurs dans chaque intervalle, et également que la moyenne ne peut être supérieure à = . Nous concluons:(1+8+16+24)/410.25(8+16+24+32)18

Mean=[10.25,18].

Cela représente un intervalle entier d'estimations: un résultat approprié d'un calcul avec des entrées d'intervalle!

Une limite de confiance supérieure à (unilatérale) de la moyenne de quatre valeurs est calculée à partir de leur moyenne et d'écart-type avec le Student t- distribution en tant que1αx=(x1,x2,x3,x4)ms

ucl(x,α)=x+tn1(α)s/n.

Contrairement au calcul de la moyenne, il n'est plus généralement le cas que l'intervalle des ucl soit limité par les ucl des valeurs limites. En effet, notez que l'ucl des limites d'intervalle inférieures, , est égal à , tandis que est encore plus petit. En maximisant et en minimisant l'ucl parmi toutes les combinaisons possibles de valeurs cohérentes avec les observations, nous constatons (par exemple) queucl((1,8,16,24),.025)28.0758ucl((8,11.676,16,24),.025)=25.8674

ucl(data,.025)=[25.8,39.3]

(c'est un intervalle de nombres représentant un ucl évalué par intervalle , pas un intervalle de confiance!) et, pour la limite de confiance inférieure,

lcl(data,.025)=[0,6.2].

(Ces valeurs ont été arrondies vers l'extérieur. Le est une valeur négative qui a été tronquée à en partant du principe que le temps médian des bourgeons ne peut pas être négatif.)00

En mots, on pourrait dire que

"Ces observations sont cohérentes avec des valeurs qui, si elles avaient été mesurées avec précision , pourraient entraîner une limite de confiance supérieure de 2,5% de la médiane pouvant atteindre 39,3 jours, mais pas plus. Elles sont cohérentes avec des valeurs (qui pourraient différer de la première) cela se traduirait par une limite de confiance inférieure de 2,5% aussi bas que 0. "

Ce que l'on doit en faire relève de la réflexion individuelle et dépend de l'application. Si l'on veut être raisonnablement sûr que le débourrement se produit avant 40 jours, alors ce résultat donne une certaine satisfaction (sous réserve des hypothèses sur la distribution du débourrement et l'indépendance des observations ). Si l'on veut estimer le débourrement au jour le plus proche, alors il est clair que davantage de données sont nécessaires. Dans d'autres circonstances, cette conclusion statistique en termes de limites de confiance à intervalles peut être frustrante. Par exemple, dans quelle mesure pouvons-nous être sûrs que le débourrement se produit dans 50% des spécimens avant 30 jours? C'est difficile à dire, car les réponses seront des intervalles.


Il existe d'autres façons de gérer ce problème. Je préfère particulièrement utiliser les méthodes du maximum de vraisemblance. (Pour les appliquer ici, nous aurions besoin d'en savoir plus sur la façon dont les seuils d'intervalle ont été établis. Il importe qu'ils aient été déterminés indépendamment des données ou non.) La présente question semble être une bonne occasion d'introduire des méthodes basées sur l'intervalle car elles ne semblent pas bien connues, même si dans certaines disciplines (évaluation des risques et analyse des algorithmes) elles ont été chaleureusement défendues par certains.

whuber
la source
Merci pour votre réponse. Les dates d'échantillonnage ont été choisies indépendamment des données (environ toutes les 1-2 semaines, quand j'ai eu la chance de me rendre sur place.
David LeBauer
J'en ai pensé autant, David, mais j'ai aussi pensé que votre capacité à faire des observations pouvait être liée aux conditions météorologiques et à d'autres facteurs qui pouvaient eux-mêmes influencer le moment du débourrement. Ainsi, bien que le processus de choix des dates d'échantillonnage ait pu être considéré comme indépendant du processus de débourrement, les deux pourraient encore avoir une forte dépendance
whuber
2
désolé, j'ai mal parlé. Mes dates d'échantillonnage ont été moins rigoureuses l'automne dernier; au printemps, toutes les dates étaient espacées de 10 jours, à l'exception des observations de première seconde avec dt = 13, mais aucun changement n'a été observé entre ces observations. À l'automne cependant, octobre-novembre était assez pluvieux; la sénescence des feuilles et les intervalles d'échantillonnage dépendaient des conditions météorologiques. (Je sais que la sénescence des feuilles dépend du temps de la biologie, cette information n'est pas dans les données).
David LeBauer
1

Voici une approche simple qui n'utilise pas de régression logistique, mais tente d'utiliser les suggestions ci-dessus. Le calcul des statistiques récapitulatives suppose, peut-être naïvement, que la date est normalement distribuée.

Veuillez pardonner le code inélégant

  1. écrire une fonction pour estimer le jour du débourrement pour chaque individu: utiliser le jour de l'année à mi-chemin entre la dernière observation de 0 et la première observation de 1 pour chaque individu.

    budburst.day <- function(i){
       data.subset <- subset(testdata, subset =
                             id == i, 
                             na.rm = TRUE)
       y1 <- data.subset$day[max(which(data.subset$obs==0))]
       y2 <- data.subset$day[min(which(data.subset$obs==1))]
       y <- mean(c(y1, y2), na.rm = TRUE)
       if(is.na(y) | y<0 | y > 180) y <- NA
       return(y)
    }
    
  2. Calculer des statistiques sommaires

    #calculate mean
    mean(unlist(lapply(1:4, budburst.day)))
    [1] 16.125  
    
    #calculate SE = sd/sqrt(n)
    sd(unlist(lapply(1:4, budburst.day)))/2
    [1] 5.06777
    
David LeBauer
la source
0

Nous savons que le temps de transition (de l'état 0 à l'état 1) du sujet était entre deux limites: . Une approximation consiste à supposer que peut avoir pris des valeurs dans cette plage avec une probabilité uniforme . En rééchantillonnant les valeurs de , nous pouvons obtenir une distribution approximative de :t1id=124<t1<32t1timedian(ti)

t = replicate(10000, median(sample(c(runif(1, 24, 32),  # id=1
                                     runif(1,  1,  8),  # id=2
                                     runif(1,  8, 16),  # id=3
                                     runif(1, 16, 24)), # id=4
                                   replace=TRUE)))
c(quantile(t, c(.025, .25, .5, .75, .975)), mean=mean(t), sd=sd(t))

Résultat (répété):

    2.5%       25%       50%       75%     97.5%      mean        sd 
4.602999 11.428310 16.005289 20.549056 28.378774 16.085808  6.243129 
4.517058 11.717245 16.084075 20.898324 28.031452 16.201022  6.219094 

Ainsi, une approximation avec un intervalle de confiance à 95% de cette médiane est de 16 (5 - 28).

EDIT: Voir le commentaire de whuber sur la limitation de cette méthode lorsque le nombre d'observations est petit (y compris n = 4 lui-même).

GaBorgulya
la source
@GaBorgulya Je pense que vous avez une faute de frappe; médiane (IC à 95%) = 16 (5,28)
David LeBauer
Vous feriez mieux avec un ajustement ML d'une forme de distribution raisonnable aux données d'intervalle suivi d'une estimation de la médiane de la distribution.
whuber
@whuber La "distribution raisonnable" est la question clé elle-même.
GaBorgulya
1
Je suis d'accord. Il me semble qu'il doit y avoir des approches non paramétriques, telles que les lissages du noyau, qui fonctionnent avec des données à intervalles.
whuber
4
Avec quatre points de données indépendants, il est impossible d'obtenir un IC à 95% pour la médiane: quelle que soit la vraie médiane, les quatre points seront tous supérieurs à celui-ci avec une probabilité = 6,25% et ils seront tous moins que cela avec la même probabilité. Par conséquent, la plage des quatre observations ne doit pas couvrir la médiane 12,5% du temps. Cela me rend méfiant de votre approche. 1/24
whuber
0

Vous pouvez utiliser un modèle d'aléa temporel discret avec régression logistique (en utilisant un ensemble de données de période-personne). Voir Analyse des données longitudinales appliquées - logiciel et chapitres 10-12 du livre .

Allison discute également

Cependant, votre ensemble de données est minuscule.

B_Miner
la source
1
Merci pour votre réponse; bien que l'exemple de jeu de données soit minuscule, le vrai jeu de données a 100 sujets mesurés sur 6 dates
David LeBauer
-1

En supposant que vous disposerez de plus de données de la même structure, vous pourrez utiliser la méthode actuarielle (table de mortalité) pour estimer la survie médiane.

GaBorgulya
la source
1
Bonne idée! - Mais pourriez-vous peut-être expliquer comment obtenir des IC pour la médiane d'une table de mortalité?
whuber