# How can I estimate the time at which 50% of a binomial variable will have transitioned?

I have the following data, representing the binary state of four subjects at four times, note that it is only possible for each subject to transition $0\to 1$ but not $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))


I can model it with a logistic regression:

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


First, how can I account for repeated measures on the same individual within the model?

Second, how can I estimate, with uncertainty, the day on which 1/2 of the subjects will have made the transition from $0\to 1$?

-
It appears there is a strong dependence in these data: namely, is it the case that if obs=1 for subject $i$ on day $t$ then necessarily obs=1 for subject $i$ on day $s$ whenever $s \ge t$? If this is so, then you really have only four data values--one for each subject--and one of them is censored on the right. –  whuber Apr 20 '11 at 18:13
@whuber you are correct about the dependence (at least in the present within-year analysis); the data represent whether or not 'bud burst' has occurred prior to the observation date for each of four replicate trees. But I am not sure what you mean about on of the data values being censored on the right? –  David Apr 20 '11 at 18:19
Here's a summary: subject 2 transitioned in the interval [1,8]; that is, 2-->[1,8]. Also 3-->[8,16], 4-->[16,24], and 1-->[24, infinity]. The latter means subject 1 was observed for 24 days without transitioning; it's the censored value. You can frame this as a survival analysis problem and analyze it accordingly. Incidentally, this dependence means the p-values in the logistic regression are misleadingly low. –  whuber Apr 20 '11 at 18:24
@whuber thank you for the insight, but does this mean that my approach if fundamentally flawed given that I am not interested in estimating p-values? Also, none of the data will be right-censored in a few weeks; I am just developing the analysis before the dataset is complete. I have altered the test data so that none of the subjects are right censored. –  David Apr 20 '11 at 18:37
@DWin, @David This is not a repeated measurement situation. The data format only makes it look like that. The measurement for each subject consists of a single interval during which a transition was observed. –  whuber Apr 20 '11 at 21:32

As became evident in comments to the question, the data consist of only four observations of time to bud burst. (It would be a mistake to analyze them as if they were 16 independent values.) They consist of intervals of times rather than exact times:

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


There are several approaches one might take. An appealing, highly general one is to take these intervals at their word: the true time of bud burst could be anything within each interval. We are thus led to represent "uncertainty" in two separate forms: sampling uncertainty (we have a presumably representative sample of the species this year) and observational uncertainty (reflected by the intervals).

Sampling uncertainty is handled with familiar statistical techniques: we are asked to estimate the median and we can do so in any number of ways, depending on statistical assumptions, and we can provide confidence intervals for the estimate. For simplicity, let's suppose time to bud burst has a symmetrical distribution. Because it is (presumably) non-negative, this implies it has a variance and also suggests the mean of even just four observations may be approximately normally distributed. Moreover, symmetry implies we can use the mean as a surrogate for the median (which is sought in the original question). This gives us access to standard, simple, estimates and confidence interval methods.

Observation uncertainty can be handled with principles of interval arithmetic (often called "probability bounds analysis"): perform all calculations using all possible configurations of data consistent with the observations. Let's see how this works in a simple case: estimating the mean. It is intuitively clear that the mean can be no smaller than $(1+8+16+24)/4$ = $10.25$, achieved by using the smallest values in each interval, and also that the mean can be no greater than $(8+16+24+32)$ = $18$. We conclude:

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

This represents an entire interval of estimates: an appropriate result of a computation with interval inputs!

A $1-\alpha$ upper (one-sided) confidence limit of the mean of four values $\mathbf{x} = (x_1, x_2, x_3, x_4)$ is computed from their mean $m$ and sample standard deviation $s$ with the Student t-distribution as

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

Unlike the calculation of the mean, it is no longer generally the case that the interval of ucl's is bounded by the ucl's of the limiting values. Indeed, note that the ucl of the lower interval limits, $\text{ucl}((1,8,16,24), .025)$, equals $28.0758$, whereas $\text{ucl}((8, 11.676, 16, 24), .025) = 25.8674$ is smaller yet. By maximizing and minimizing the ucl among all possible combinations of values consistent with the observations, we find (for example) that

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

(that's an interval of numbers representing an interval-valued ucl, not a confidence interval!) and, for the lower confidence limit,

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

(These values have been rounded outwards. The $0$ is a negative value that was truncated to $0$ on the premise that the median bud time cannot be negative.)

In words, we might say that

"These observations are consistent with values that, had they been precisely measured, could result in an upper 2.5% confidence limit of the median as high as 39.3 days, but no higher. They are consistent with values (which might differ from the first) that would result in a lower 2.5% confidence limit as low as 0."

What one is to make of this is a matter for individual contemplation and depends on the application. If one wants to be reasonably sure that bud burst occurs before 40 days, then this result gives some satisfaction (conditional on the assumptions about bud burst distribution and independence of the observations). If one wants to estimate bud burst to the nearest day, then clearly more data are needed. In other circumstances, this statistical conclusion in terms of interval-valued confidence limits may be frustrating. E.g., how confident can we be that bud burst occurs in 50% of specimens before 30 days? It's hard to say, because the answers will be intervals.

There are other ways to handle this problem. I especially favor using maximum likelihood methods. (To apply them here, we would need to know more about how the interval cutpoints were established. It matters whether they were determined independently of the data or not.) The present question appears to be a good opportunity to introduce interval-based methods because they do not seem to be well known, even though in certain disciplines (risk assessment and analysis of algorithms) they have been warmly advocated by some people.

-
thank you for your answer. The sampling dates were chosen independently of the data (approximately every 1-2 weeks, when I had a chance to get out there. –  David Mar 15 '12 at 0:19
I figured as much, David, but it also occurred to me that your ability to make observations could be related to weather conditions and other factors that themselves might influence the time of bud burst. So, although the process of choosing the sampling dates may have been considered to be independent of the process of bud burst, the two could still have a strong statistical dependence. –  whuber Mar 15 '12 at 15:08
sorry, I mis-spoke. My sampling dates were less rigorous last fall; in the spring all of the dates were 10 days apart, excluding the first-second observations with dt = 13, but there were no changes between these observations. In the fall however, Oct-Nov was quite rainy; both leaf senescence and the sampling intervals were dependent on weather. (I know that leaf senescence is dependent on weather from biology, this information is not in the data). –  David Mar 15 '12 at 18:47

We know that the $t_1$ transition time (from state 0 to state 1) of subject id=1 was between two boundaries: $24<t_1<32$. An approximation is to assume that $t_1$ may have taken values within this range with uniform probability. Resampling the $t_i$ values we can get an approximate distribution of $\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))


Result (repeated):

    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


Thus an approximation with 95% confidence interval of this median is 16 (5 – 28).

EDIT: See whuber's comment on the limitation of this method when the number of observations is small (including n=4 itself).

-
@GaBorgulya I think you have a typo; median(95%CI) = 16(5,28) –  David Apr 20 '11 at 19:45
@David Corrected, thanks. –  GaBorgulya Apr 20 '11 at 19:58
You would do better with an ML fit of a reasonable distributional form to the interval data followed by an estimate of the median of the distribution. –  whuber Apr 20 '11 at 21:33
@whuber The "reasonable distribution" is the key question itself. –  GaBorgulya Apr 20 '11 at 21:44
With four independent data points it is impossible to obtain a 95% CI for the median: no matter what the true median is, the four points will all be greater than it with probability $1/2^4$ = 6.25% and they will all be less than it with the same probability. Therefore the range of the four observations must fail to cover the median 12.5% of the time. This makes me suspicious of your approach. –  whuber Apr 20 '11 at 21:54

Here is a simple approach that does not use logistic regression, but does attempt to use the suggestions above. Calculation of summary stats assumes, perhaps naively, that the date is normally distributed.

1. write a function to estimate the day of budbreak for each individual: use the day of year half way between the last observation of 0 and the first observation of 1 for each individual.

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 summary statistics

#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

-

You could use discrete time hazard model fit with logistic regression (using a person period data set). See Applied Longitudinal Data Analysis - software and Book Chapters 10-12.

Allison also discusses

You data set is tiny though.

-
thank you for your answer; although the example dataset is tiny, the real dataset has 100 subjects measured on 6 dates –  David Jun 27 '11 at 16:33