R: Relación entre la hora de la mortalidad y la esperanza de vida dado que esto es aproximada por una distribución exponencial?

Question

Tengo 33 los registros de las esperanzas de vida (en horas) para las abejas obreras

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

Suponga que la duración de la vida son aproximar por una distribución exponencial, ¿cómo puedo calcular el salario por hora tasa de mortalidad de las abejas utilizando exponencial de aproximación y de máxima verosimilitud?

Yo estaba pensando en usar el siguiente código:

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

pero estoy muy confundido acerca de la tasa parámetro en la función de densidad de probabilidad de la distribución exponencial (pensé que la tasa correspondiente a max probabilidades debe ser 1/media(x) de aquí?) y no sé cómo cada hora de mortalidad encaja en cualquier lugar utilizando un enfoque de máxima verosimilitud.

Answer 1

La función de riesgo se $\frac{f(t)}{1-F(t)}$.

Suponiendo que toda la vida a ser exponencial, la fórmula es dado en el enlace; es $\lambda$, la tasa parámetro de la exponencial; se estima que la tasa parámetro $\lambda$ a partir de los datos exactamente de la misma manera que lo haría para cualquier distribución exponencial.

La hora de la tasa de mortalidad es la integral del riesgo de la tasa a través de una hora, pero debido a que el peligro es constante y los tiempos se miden en horas, es trivial.

--

La discusión sobre la idoneidad del modelo:

Si la esperanza de vida eran realmente exponencial, el peligro sería constante. Pero que no se ve exponencial para mí.

Aunque es posible obtener una estimación no paramétrica de la función de riesgo, con tan pocos datos, probablemente me inicio con un modelo Weibull, ya que incluye su sugirió exponencial modelo de vida como un caso especial y puede acomodar tanto el aumento y la disminución de funciones de riesgo.

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)

El ajuste parece bastante razonable.

Además, si nos fijamos en el error estándar de la log(escala) coeficiente, se sugiere que la exponencial (que corresponde a la sesión escala de 0 si entiendo todo a la derecha) parece inconsistente con los datos:

             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 

Si lo hizo bien, la equipada con función de riesgo $h(t) = f(t)/(1-F(t))$ debería tener este aspecto:

Otro modelo que yo podría considerar, que también incluye la exponencial como un caso especial sería el de gamma. También ha ido en aumento y la disminución de riesgos. Sin embargo, creo que no podría hacer nada mejor en este caso.

Answer 2

La pregunta indaga acerca de la relación entre la hora de la mortalidad y la distribución exponencial de la duración de la vida.

El universo se divide el tiempo en corto garrapatas. Para cada abeja un dado es lanzado. El dado tiene muchos lados, la mayoría de ellas en blanco ... pero una es de color negro. La muerte viene de la abeja siempre el lado negro de la muestra.

Estamos a punto de simular esta con el equipo, pero la pereza no quiero rebanada de tiempo tan finamente como hace la Muerte. Así que en lugar de rodar una mueren cada segundo, digo, tal vez me va a rodar una cada hora mueren: pero voy a hacer 3600 veces más propensos a mostrar su lado negro para compensar la escasez de rollos. Debido a que la posibilidad de la muerte en una segunda o incluso una hora es tan pequeño, esta es una buena aproximación.

Para explotar la capacidad de la computadora para hacer un montón de cálculos con matrices a la vez, voy a simular toda una colmena una abeja en un momento. Vamos a la tirada de dados una y otra vez, una vez por hora, hasta que la abeja muere. Luego continuar con un recién nacido (reencarnar?) abeja, repitiendo hasta muchas abejas nacen y mueren. El resultado será el mismo como si todas las abejas han coexistido.

Aquí es R código. Su entrada es la matriz de abeja vidas, bees, utilizado para la estimación de la velocidad exponencial.

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

Aquí está la primera parte de la deaths de la matriz, mostrando los resultados de comienzo morir rollos (con "*" denota la cara negra):

paste(ifelse(head(deaths, 100),"*", "."), collapse="")

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

El primer par de veces en el que el color negro de la cara apareció en el morir:

head(times)

12 67 82 120 134 147

Así, la primera abeja fue derribado en el 12 de garrapata (así que me cuente su vida como sólo 11 de las garrapatas), el segundo en la 67ª sesión de la garrapata (para toda la vida de 54 garrapatas), y así sucesivamente: estas son las posiciones de las estrellas en la salida anterior.

Las cifras de la parcela estas muertes (cada uno es una barra vertical a la izquierda, mostrando cada uno de toda la vida como un espacio horizontal entre barras) y la distribución de los tiempos de vida (a la derecha). El último tiene una distribución exponencial la función superpuesta. Es un excelente ajuste.

El código para graficar el histograma es

hist(lifetimes, freq=FALSE, ylim=c(0, rate.hat), breaks=25, xlab="Hours")
curve(dexp(x, rate=rate.hat), col="Red", add=TRUE)

¿Por qué es el histograma exponencial? Debido a que en promedio el número de abejas vivas después de cierta edad se reduce en la misma fracción durante el siguiente intervalo de tiempo. Eso significa que el histograma tiene a disminuir exponencialmente.

Por cierto, la razón por la 1/mean(bees) es el estimador de máxima verosimilitud es que la exponencial de la función de distribución de probabilidad es de la forma $\kappa \exp(-\kappa t)$ toda la vida $t$. Para entender esto, recordemos que un PDF es un densidad: nos da la probabilidad por unidad de tiempo. Desde $t$ es el tiempo (en horas), $\kappa$ debe ser una probabilidad (de la muerte) por unidad de tiempo: es la tasa de interés. Su logaritmo es igual a $\log(\kappa) - \kappa t$. Por lo tanto, el registro de probabilidad de un conjunto de datos de tiempos de vida de $(t_1, t_2, \ldots, t_n)$ es

$$\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}$$ es el único punto crítico (y que, obviamente, es donde la probabilidad es maximizada). Este es el recíproco de la media de toda la vida. Que es como hora de mortalidad encaja con el de máxima verosimilitud de la maquinaria.

Answer 3

He aquí un enfoque

Obtener el histograma de los datos

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)

Graficar el histograma como puntos de

plot(b1h$mids,b1h$density,cex=1,xlab='Life (hr)',ylab='Probability')

Ajuste de una distribución exponencial para el histograma

fit <- nls(density ~ lambda*exp(-lambda*mids), start=list(lambda=.1), data=df)
lines(df$mids,predict(fit),col='red')

Alternativamente, se ajuste a una distribución exponencial a los datos

fit1 <- fitdistr(bl, "exponential") 
curve(dexp(x, rate = fit1$estimate), col = "green", add = TRUE)

La vida media es de 1/lambda, y la tasa de mortalidad es de 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')

Los resultados son

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

Función fitdisr proporciona la estimación por máxima verosimilitud de los parámetros. El nls función proporciona una estimación de mínimos cuadrados. Nota, para una distribución exponencial, la estimación de máxima verosimilitud está dada por n/suma(xi), ref derivación