Como posso estimar o tempo em que 50% de uma variável binomial terá transitado?

Eu tenho os seguintes dados, representando o estado binário de quatro sujeitos quatro vezes, observe que só é possível para cada sujeito fazer a transição mas não 1 → 0 :$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))

Eu posso modelá-lo com uma regressão logística:

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

Primeiro, como posso explicar medidas repetidas no mesmo indivíduo dentro do modelo?

Segundo, como posso estimar, com incerteza, o dia em que 1/2 dos sujeitos fará a transição de ?$0\to 1$

Como ficou evidente nos comentários à pergunta, os dados consistem em apenas quatro observações de tempo para o início da explosão. (Seria um erro analisá-los como se fossem 16 valores independentes.) Eles consistem em intervalos de tempo e não em tempo exato:

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

Existem várias abordagens que se pode adotar. Um recurso atraente e altamente geral é levar em consideração esses intervalos: o tempo real de rompimento de botões pode ser qualquer coisa dentro de cada intervalo. Somos, assim, levados a representar a "incerteza" de duas formas distintas: incerteza de amostragem (temos uma amostra presumivelmente representativa das espécies este ano) e incerteza observacional (refletida pelos intervalos).

A incerteza de amostragem é tratada com técnicas estatísticas familiares: somos solicitados a estimar a mediana e podemos fazê-lo de várias maneiras, dependendo das suposições estatísticas, e podemos fornecer intervalos de confiança para a estimativa. Por uma questão de simplicidade, vamos supor que o tempo para brotar burst tenha uma distribuição simétrica. Por ser (presumivelmente) não-negativo, isso implica que há uma variação e também sugere que a média de apenas quatro observações pode ser aproximadamente normalmente distribuída. Além disso, a simetria implica que podemos usar a média como substituta da mediana (o que é procurado na pergunta original). Isso nos dá acesso a métodos padrão, simples, estimativas e intervalo de confiança.

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

Isso representa um intervalo inteiro de estimativas: um resultado apropriado de um cálculo com entradas de intervalo!

$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$ainda é menor. Ao maximizar e minimizar o ucl entre todas as combinações possíveis de valores consistentes com as observações, descobrimos (por exemplo) que

(é um intervalo de números que representa um ucl avaliado com intervalo , não um intervalo de confiança!) e, para o limite de confiança mais baixo,

$0$$0$

Em palavras, podemos dizer que

"Essas observações são consistentes com valores que, se tivessem sido medidos com precisão , poderiam resultar em um limite superior de confiança de 2,5% da mediana tão alta quanto 39,3 dias, mas não mais alto. Eles são consistentes com valores (que podem diferir do primeiro) isso resultaria em um limite de confiança de 2,5% inferior a 0 ".

O que se deve fazer disso é uma questão de contemplação individual e depende da aplicação. Se alguém quiser ter certeza razoável de que o rebentamento de brotos ocorre antes de 40 dias, esse resultado dará alguma satisfação ( condicional às suposições sobre a distribuição do rebento de brotos e a independência das observações ). Se alguém quiser estimar o número de brotos até o dia mais próximo, são claramente necessários mais dados. Em outras circunstâncias, essa conclusão estatística em termos de limites de confiança com valor de intervalo pode ser frustrante. Por exemplo, quão confiantes podemos estar que a brotação ocorre em 50% das amostras antes de 30 dias? É difícil dizer, porque as respostas serão intervalos.

Existem outras maneiras de lidar com esse problema. Sou especialmente favorável ao uso de métodos de máxima verossimilhança. (Para aplicá-los aqui, precisaríamos saber mais sobre como os pontos de corte de intervalo foram estabelecidos. Importa se eles foram determinados independentemente dos dados ou não.) A presente pergunta parece ser uma boa oportunidade para introduzir métodos baseados em intervalos, porque elas não parecem ser bem conhecidas, embora em certas disciplinas (avaliação de riscos e análise de algoritmos) tenham sido calorosamente defendidas por algumas pessoas.

Aqui está uma abordagem simples que não usa regressão logística, mas tenta usar as sugestões acima. O cálculo das estatísticas de resumo pressupõe, talvez ingenuamente, que a data seja normalmente distribuída.

Perdoe código deselegante

1. escreva uma função para estimar o dia de brotação de cada indivíduo: use o dia do ano a meio caminho entre a última observação de 0 e a primeira observação de 1 para cada indivíduo.

   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. Calcular estatísticas resumidas

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

Resultado (repetido):

    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 

Assim, uma aproximação com intervalo de confiança de 95% dessa mediana é 16 (5 - 28).

EDIT: Veja o comentário do whuber sobre a limitação deste método quando o número de observações é pequeno (incluindo o próprio n = 4).

Você pode usar um modelo discreto de risco de tempo adequado à regressão logística (usando um conjunto de dados de período por pessoa). Consulte Análise longitudinal de dados aplicada - software e capítulos 10-12 do livro .

Allison também discute

Seu conjunto de dados é minúsculo.

Supondo que você tenha mais dados da mesma estrutura, poderá usar o método atuarial (tabela de vida) para estimar a sobrevida média.