R: Relationship between hourly mortality and lifespan given that the latter is approximated by an exponential distribution? I have 33 records of lifespans (in hours) for worker bees
bee[,2]
[1]  7.1  2.3  9.6 25.8 14.6 12.8 20.9 30.0 71.1 36.9  3.9 18.2 24.3 55.8 84.7
[16] 13.6 10.8 18.2 19.0 34.8 27.1 54.2 34.8 65.8 33.3 26.5  9.5 44.2 35.7  2.3
[31]  4.0 25.9 41.3

Assume that the lifespans are approximated by an exponential distribution, how can I estimate the hourly mortality rate of the bees using exponential approximation and maximum likelihood?
I was thinking about using the following code:
p <- seq(from = 0.01, to = 1, by = 0.001)
x <- bee[,2]
like <- vector()
for (i in 1:length(p)){like[i] <- sum(dexp(x, rate = p[i]))}
plot(like ~ p)
(p[which(like == max(like))]) # which gives 0.111

but I am awfully confused about the rate parameter in the probability density function of  exponential distribution (I thought the rate corresponding to max likelihood should be 1/mean(x) here?) and I can't tell how hourly mortality fits in anywhere using an maximum likelihood approach. 
 A: The hazard function will be $\frac{f(t)}{1-F(t)}$.
Assuming lifetime to be exponential, the formula is given at that link; it's $\lambda$, the rate parameter of the exponential; you estimate that rate parameter $\lambda$ from the data exactly the same way you would for any exponential distribution. 
The hourly mortality rate is the integral of the hazard rate across an hour, but
since the hazard is constant and times are measured in hours, it's trivial.
--
Discussion on the suitability of the model:
If lifespan were truly exponential, hazard would be constant. But that doesn't look exponential to me.
While it's possible to get a nonparametric estimate of the hazard function, with so little data, I'd probably start with a Weibull model, since it includes your suggested exponential life model as a special case and can accommodate both increasing and decreasing hazard functions.
survreg(Surv(beelife)~1)
Call:
survreg(formula = Surv(beelife) ~ 1)

Coefficients:
(Intercept) 
   3.416515 

Scale= 0.7265582 

Loglik(model)= -140.5   Loglik(intercept only)= -140.5
n= 33


 plot(ecdf(beelife))
 lines((0:90),pweibull(0:90,1/.726558,exp(3.4165)),col=2)


The fit looks pretty reasonable.
Further, if we look at the standard error of the log(scale) coefficient, it suggests that exponential (which corresponds to log-scale of 0 if I understand everything right) seems inconsistent with the data:
             Value Std. Error     z         p
(Intercept)  3.417      0.133 25.66 3.50e-145
Log(scale)  -0.319      0.137 -2.32  2.01e-02

Scale= 0.727 

Weibull distribution
Loglik(model)= -140.5   Loglik(intercept only)= -140.5
Number of Newton-Raphson Iterations: 6 
n= 33 

If I did it right, the fitted hazard function $h(t) = f(t)/(1-F(t))$ should look like this:

Another model I might consider, which also includes the exponential as a special case would be the gamma. It also has increasing and decreasing hazards. However, I doubt it would do any better in this case.
A: The question asks about the relationship between hourly mortality and exponential distribution of lifespans.
The universe divides time into short ticks.  For each bee a die is rolled.  The die has many sides, most of them blank--but one is black.  Death comes for the bee whenever the black side shows.
We are about to simulate this with the computer, but being lazy I don't want to slice time as finely as Death does.  So instead of rolling a die every second, say, maybe I will roll a die every hour: but I will make it 3600 times more likely to show its black side to make up for the paucity of rolls.  Because the chance of death in one second or even one hour is so small, this is a good approximation.
To exploit the computer's ability to do a lot of calculations with arrays at once, I will simulate an entire hive one bee at a time.  Let the die roll again and again, once per hour, until the bee dies.  Then continue with a newly born (reincarnated?) bee, repeating until many bees are born and die.  The results will be the same as if all the bees had coexisted.
Here is R code.  Its input is the array of bee lifetimes, bees, used to estimate the exponential rate.
set.seed(17)                          # (Allows reproducible results)
n <- 10^5                             # Number of time slices
units <- 1                            # Number of hours per time slice
rate.hat <- 1 / mean(bees)            # ML estimate of the death rate
deaths <- runif(n) < rate.hat * units # Roll Death's die
times <- which(deaths != 0)           # Note when bees are killed
lifetimes <- (diff(c(0,times))-1) * units # The time differences are the bee lifetimes

Here is the first part of the deaths array, showing the results of the beginning die rolls (with "*" denoting the black face):
paste(ifelse(head(deaths, 100),"*", "."), collapse="")


...........*
  ......................................................*
  ..............*
  ..................

The first few times at which the black face appeared on the die:
head(times)


12  67  82 120 134 147

Thus the first bee was felled at the 12th tick (so I count its lifetime as just 11 ticks), the second one at the 67th tick (for a lifetime of 54 ticks), and so on: these are the positions of the stars in the previous output.
The figures plot these deaths (each is a vertical slash on the left, showing each lifetime as a horizontal gap between slashes) and the distribution of the lifetimes (on the right).  The latter has an exponential distribution function superimposed.  It's an excellent fit.

The code to plot the histogram is
hist(lifetimes, freq=FALSE, ylim=c(0, rate.hat), breaks=25, xlab="Hours")
curve(dexp(x, rate=rate.hat), col="Red", add=TRUE)

Why is the histogram exponential?  Because on average the number of bees alive after a certain age is reduced by the same fraction during the next time slice.  That means the histogram has to decrease exponentially.
By the way, the reason why 1/mean(bees) is the maximum likelihood estimator is that the exponential probability distribution function is of the form $\kappa \exp(-\kappa t)$ for lifetimes $t$.  To understand this, recall that a PDF is a density: it gives us probability per unit time.  Since $t$ is time (in hours), $\kappa$ must be a probability (of death) per unit time: it's the rate.  Its logarithm equals $\log(\kappa) - \kappa t$.  Therefore the log likelihood for a dataset of lifetimes $(t_1, t_2, \ldots, t_n)$ is
$$\log(L(\kappa)) = n \log(\kappa) - \sum_{i=1}^n \kappa\, t_i.$$
Calculus shows us (by taking the derivative with respect to $\kappa$ and setting that to zero) that $$\hat\kappa = \frac{n}{\sum_{i=1}^n t_i}$$ is the only critical point (and obviously is where the likelihood is maximized).  This is the reciprocal of the mean lifetime.  That is how hourly mortality fits in with the maximum likelihood machinery.
A: Here's an approach
Get the histogram of the data
bl <- c(7.1, 2.3, 9.6, 25.8, 14.6, 12.8, 20.9, 30.0, 71.1, 36.9, 3.9, 18.2, 24.3, 55.8, 84.7,
          13.6, 10.8, 18.2, 19.0, 34.8, 27.1, 54.2, 34.8, 65.8, 33.3, 26.5, 9.5, 44.2, 35.7, 2.3,
          4.0, 25.9, 41.3)

b1h <- hist(bl, freq = FALSE, xlim = c(0, quantile(bl, 0.999)), plot=TRUE)
df <- data.frame(density=b1h$density, mids=b1h$mids)

Plot the histogram as points
plot(b1h$mids,b1h$density,cex=1,xlab='Life (hr)',ylab='Probability')

Fit an exponential distribution to the histogram
fit <- nls(density ~ lambda*exp(-lambda*mids), start=list(lambda=.1), data=df)
lines(df$mids,predict(fit),col='red')

Alternatively, fit a exponential distribution to the data
fit1 <- fitdistr(bl, "exponential") 
curve(dexp(x, rate = fit1$estimate), col = "green", add = TRUE)

The mean life is 1/lambda, and the mortality rate is lambda
cat('from nls (red) the mean life is',1/coef(fit),'and mortality rate is',coef(fit),'per year')
cat('from fitdistr (green) the mean life is',1/coef(fit1),'and mortality rate is',coef(fit1),'per year')

The results are
from nls (red) the mean life is 35.47136 and mortality rate is 0.02819176 per year
from fitdistr (green) the mean life is 27.84848 and mortality rate is 0.0359086 per year


Function fitdisr provides the maximum likelihood estimate for the parameters.  The nls function provides a least square estimate.  Note, for an exponential distribution, the maximum likelihood estimate is given by n/sum(xi), ref derivation
