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$?
 A: 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.
A: 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. 
Please pardon inelegant code


*

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


*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

A: 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).
A: 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.
A: Assuming that you will have more data of the same structure you will be able to use the actuarial (life table) method to estimate median survival.
