Wie kann ich den Zeitpunkt abschätzen, zu dem 50% einer Binomialvariablen übergegangen sind?

Ich habe die folgenden Daten, die den Binärzustand von vier Subjekten zu vier Zeiten darstellen. Beachten Sie, dass es nur möglich ist, dass jedes Subjekt , nicht aber 1 → 0 übergeht :$0\to 1$$1\to 0$

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

Ich kann es mit einer logistischen Regression modellieren:

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 ***

Erstens, wie kann ich wiederholte Messungen an derselben Person innerhalb des Modells berücksichtigen?

Zweitens, wie kann ich mit Unsicherheit den Tag abschätzen, an dem die Hälfte der Probanden den Übergang von vollzogen hat?$0\to 1$

Wie aus den Kommentaren zur Frage hervorgeht, bestehen die Daten nur aus vier Beobachtungen der Zeit bis zum Knospenausbruch. (Es wäre ein Fehler, sie so zu analysieren, als wären sie 16 unabhängige Werte.) Sie bestehen eher aus Zeitintervallen als aus genauen Zeiten:

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

Es gibt verschiedene Ansätze. Eine ansprechende, sehr allgemeine ist es, diese Intervalle beim Wort zu nehmen: Die wahre Zeit des Knospenausbruchs kann alles innerhalb jedes Intervalls sein. Wir werden daher dazu gebracht, "Unsicherheit" in zwei getrennten Formen darzustellen: Stichprobenunsicherheit (wir haben dieses Jahr eine vermutlich repräsentative Stichprobe der Arten) und Beobachtungsunsicherheit (die sich in den Intervallen widerspiegelt).

Die Stichprobenunsicherheit wird mit bekannten statistischen Techniken behandelt: Wir werden gebeten, den Median zu schätzen, und wir können dies abhängig von statistischen Annahmen auf verschiedene Arten tun, und wir können Konfidenzintervalle für die Schätzung bereitstellen. Nehmen wir zur Vereinfachung an, dass die Zeit bis zum Knospenausbruch symmetrisch verteilt ist. Da es (vermutlich) nicht negativ ist, impliziert dies eine Varianz und legt nahe, dass der Mittelwert von nur vier Beobachtungen ungefähr normalverteilt sein kann. Darüber hinaus impliziert Symmetrie, dass wir den Mittelwert als Ersatz für den Median verwenden können (der in der ursprünglichen Frage gesucht wird). Dies gibt uns Zugang zu Standard-, einfachen Schätz- und Konfidenzintervallmethoden.

$(1+8+16+24)/4$$10.25$$(8+16+24+32)$$18$

Dies stellt ein ganzes Intervall von Schätzungen dar: ein geeignetes Ergebnis einer Berechnung mit Intervalleingaben!

$1-\alpha$$\mathbf{x} = (x_1, x_2, x_3, x_4)$$m$$s$

$\text{ucl}((1,8,16,24), .025)$$28.0758$$\text{ucl}((8, 11.676, 16, 24), .025) = 25.8674$ist noch kleiner. Durch Maximieren und Minimieren des ucl unter allen möglichen Wertekombinationen, die mit den Beobachtungen übereinstimmen, finden wir (zum Beispiel), dass

( Dies ist ein Zahlenintervall, das ein Intervall mit ucl-Wert darstellt, kein Konfidenzintervall!) und für die untere Konfidenzgrenze

$0$$0$

In Worten könnten wir das sagen

"Diese Beobachtungen stimmen mit Werten überein, die, wenn sie genau gemessen worden wären , zu einer oberen 2,5% -Konfidenzgrenze des Medians von bis zu 39,3 Tagen führen könnten, jedoch nicht höher. Sie stimmen mit Werten überein (die von den ersten abweichen könnten). das würde zu einer unteren Konfidenzgrenze von 2,5% von nur 0 führen. "

Was man daraus machen soll, ist eine Frage der individuellen Betrachtung und hängt von der Anwendung ab. Wenn man einigermaßen sicher sein möchte, dass ein Knospenausbruch vor 40 Tagen auftritt, gibt dieses Ergebnis eine gewisse Befriedigung ( abhängig von den Annahmen über die Verteilung des Knospenausbruchs und der Unabhängigkeit der Beobachtungen ). Wenn man den Knospenausbruch auf den nächsten Tag abschätzen möchte, werden deutlich mehr Daten benötigt. Unter anderen Umständen kann diese statistische Schlussfolgerung in Bezug auf Intervallwert-Konfidenzgrenzen frustrierend sein. Zum Beispiel, wie sicher können wir das Austrieb tritt in 50% der Proben vor 30 Tagen sein? Es ist schwer zu sagen, weil die Antworten Intervalle sein werden.

Es gibt andere Möglichkeiten, um dieses Problem zu lösen. Ich bevorzuge besonders die Verwendung von Maximum-Likelihood-Methoden. (Um sie hier anzuwenden, müssten wir mehr darüber wissen, wie die Intervallschnittpunkte festgelegt wurden. Es ist wichtig, ob sie unabhängig von den Daten bestimmt wurden oder nicht.) Die vorliegende Frage scheint eine gute Gelegenheit zu sein, intervallbasierte Methoden einzuführen, weil Sie scheinen nicht bekannt zu sein, obwohl sie in bestimmten Disziplinen (Risikobewertung und Analyse von Algorithmen) von einigen Menschen nachdrücklich befürwortet wurden.

Hier ist ein einfacher Ansatz, der keine logistische Regression verwendet, sondern versucht, die obigen Vorschläge zu verwenden. Bei der Berechnung der zusammenfassenden Statistiken wird möglicherweise naiv davon ausgegangen, dass das Datum normal verteilt ist.

Bitte entschuldigen Sie den uneleganten Code

1. Schreiben Sie eine Funktion, um den Tag des Knospenbruchs für jede Person zu schätzen: Verwenden Sie den Tag des Jahres auf halbem Weg zwischen der letzten Beobachtung von 0 und der ersten Beobachtung von 1 für jede Person.

   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. Berechnen Sie die zusammenfassende Statistik

   #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
   

$t_1$id=1$24<t_1<32$$t_1$$t_i$$\text{median}(t_i)$

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

Ergebnis (wiederholt):

    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 

Somit beträgt eine Annäherung mit einem 95% -Konfidenzintervall dieses Medians 16 (5 - 28).

EDIT: Siehe Whubers Kommentar zur Einschränkung dieser Methode, wenn die Anzahl der Beobachtungen gering ist (einschließlich n = 4 selbst).

Sie können ein diskretes Zeitrisikomodell verwenden, das mit der logistischen Regression übereinstimmt (unter Verwendung eines Personendatensatzes). Siehe Angewandte Longitudinal Datenanalyse - Software und Buchkapitel 10-12.

Allison diskutiert auch

Ihr Datensatz ist jedoch winzig.

Unter der Annahme, dass Sie mehr Daten derselben Struktur haben, können Sie die versicherungsmathematische Methode (Lebenstabelle) verwenden , um das mediane Überleben abzuschätzen.