正如对该问题的评论所表明的那样，数据仅包含四个对萌芽时间的观察。（将它们当作 16 个独立值来分析是错误的。）它们由时间间隔而不是确切时间组成：

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

有几种方法可以采取。一个吸引人的、高度通用的方法是按照他们的话来接受这些间隔：萌芽的真实时间可以是每个间隔内的任何时间。因此，我们被引导以两种不同的形式来表示“不确定性”：采样不确定性（我们今年可能有该物种的代表性样本）和观测不确定性（由间隔反映）。

使用熟悉的统计技术处理抽样不确定性：我们被要求估计中位数，我们可以根据统计假设以多种方式进行估计，并且我们可以提供估计的置信区间。为简单起见，假设萌芽时间呈对称分布。因为它（可能）是非负的，这意味着它具有方差，并且还表明即使只有四个观察值的平均值也可能近似正态分布。此外，对称性意味着我们可以使用均值作为中位数的替代值（在原始问题中寻求）。这使我们可以使用标准、简单的估计和置信区间方法。

观察不确定性可以用区间算术原理（通常称为“概率界限分析”）处理：使用与观察一致的所有可能的数据配置执行所有计算。 让我们看看这在一个简单的情况下是如何工作的：估计平均值。直观上可以看出均值不小于 =，通过使用每个区间的最小值来实现，并且均值不大于 =。我们得出结论：$(1+8+16+24)/4$$10.25$$(8+16+24+32)$$18$

$$\text{Mean} = [10.25, 18].$$

这代表了整个估计区间：使用区间输入计算的适当结果！

的平均值的上限（单边）置信限是根据它们的平均值和样本标准差计算的，学生 t-分布为$1-\alpha$$\mathbf{x} = (x_1, x_2, x_3, x_4)$$m$$s$

$$\text{ucl}(\mathbf{x}, \alpha) = x + t_{n-1}(\alpha) s / \sqrt{n}.$$

与平均值的计算不同，通常情况下，ucl 的区间不再受限制值的 ucl 的限制。实际上，请注意区间下限的 ucl，等于，而更小。通过在与观察结果一致的所有可能的值组合中最大化和最小化 ucl，我们发现（例如）$\text{ucl}((1,8,16,24), .025)$$28.0758$$\text{ucl}((8, 11.676, 16, 24), .025) = 25.8674$

$$\text{ucl}(\text{data},.025) = [25.8, 39.3]$$

（这是代表区间值 ucl 的数字区间，而不是置信区间！）并且，对于置信下限，

$$\text{lcl}(\text{data},.025) = [0, 6.2].$$

（这些值已经向外四舍五入是负值，在中位芽时间不能为负的前提下$0$$0$

换句话说，我们可以说

“这些观察结果与值一致，如果它们被精确测量，可能会导致中位数的 2.5% 置信上限高达 39.3 天，但不会更高。它们与值一致（可能与第一个不同）这将导致 2.5% 的置信限低至 0。”

对此的理解是个人思考的问题，取决于应用程序。如果想合理地确定萌芽发生在 40 天之前，那么这个结果会让人满意（以萌芽分布和观察独立性的假设为条件）。如果要估计最近一天的芽爆发，那么显然需要更多的数据。在其他情况下，这个关于区间值置信限的统计结论可能令人沮丧。 例如，我们对 50% 的标本在 30 天前发生芽爆的信心有多大？很难说，因为答案将是间隔。

---

还有其他方法可以处理这个问题。我特别喜欢使用最大似然法。（要在这里应用它们，我们需要更多地了解区间截点是如何建立的。它们是否独立于数据确定很重要。）目前的问题似乎是引入基于区间的方法的好机会，因为它们似乎并不为人所知，尽管在某些学科（风险评估和算法分析）中它们受到了一些人的热烈提倡。

这是一种不使用逻辑回归的简单方法，但确实尝试使用上述建议。摘要统计数据的计算可能天真地假设日期是正态分布的。

请原谅不优雅的代码

1. 编写一个函数来估计每个个体的萌芽日期：使用一年中的日期，介于每个个体的最后一次观察 0 和第一次观察 1 之间。

   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. 计算汇总统计

   #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
   

我们知道受试者的转换时间（从状态 0 到状态 1）介于两个边界之间：。一个近似值是假设可能以均匀的概率在此范围内取值。重新采样值，我们可以得到的近似分布：$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))

结果（重复）：

    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 

因此，该中位数的 95% 置信区间的近似值为 16 (5 – 28)。

编辑：当观察数量很少（包括 n=4 本身）时，请参阅 whuber 对这种方法的限制的评论。

您可以使用符合逻辑回归的离散时间风险模型（使用人员周期数据集）。请参阅应用纵向数据分析 - 软件 和书籍第 10-12 章。

艾莉森还讨论了

你的数据集很小。