¿Cómo compruebo la tendencia en la suma de variables aleatorias?

Question

Tengo un conjunto de datos con mediciones de costes para ingresos hospitalarios de ~25.000 pacientes durante varios años. El coste parece tener una distribución logarítmica normal y al tomar el logaritmo se obtiene una distribución aproximadamente normal. Utilizando esta transformación, he analizado las diferencias en el coste medio a lo largo de los años utilizando un ANOVA de una vía.

¿Puedo hacer una prueba similar sobre los costes totales (es decir, la suma) a lo largo de los años?

Lo único que se me ocurre es probar una tendencia lineal mediante una regresión lineal de las sumas sobre el año. Sin embargo, esto me deja con sólo 8 en lugar de 25.000 puntos de datos, lo que intuitivamente me parece tirar mucha información.

Answer 1

Yo sugeriría modelar el coste total $Y$ en un año determinado como una variable aleatoria compuesta de Poisson. En concreto, podemos modelar el número de admisiones como una variable aleatoria de Poisson $N$ con media (desconocida) $\lambda$ y los costes de admisión como $X_1, \dots, X_N$ para que $$Y = \sum_{i=1}^N X_i$$ Aquí asumimos $X_1, X_2, \dots$ son iid (e independientes de $N$ ) con primer y segundo momento \begin{align*} \mu_1 &= E(X_i)\\ \mu_2 &= E(X_i^2) \end{align*} Si $\lambda$ es grande (lo cual es razonable suponer), la teoría asintótica muestra que $Y$ es aproximadamente gaussiano, con los siguientes parámetros: \begin{align*} E(Y) &= E(N)E(X) = \lambda\mu_1\\ \text{Var}(Y) &= E(N)E(X^2) = \lambda\mu_2 \end{align*} Por lo tanto, $(Y - E(Y))/\sqrt{\text{Var}(Y)}$ tiene aproximadamente una distribución gaussiana estándar. Podemos utilizar este hecho para derivar inferencias sobre $E(Y)$ . Escribamos $\theta = E(Y)$ ya que este es el parámetro de interés. Obsérvese que podemos expresar $\text{Var}(Y)$ en términos de $\theta$ como $$\text{Var}(Y) = \theta\frac{\mu_2}{\mu_1}$$ Aquí $\mu_1$ y $\mu_2$ puede ser aproximado por los momentos de la muestra: \begin{align*} \hat\mu_1 &= \frac1N\sum_{i=1}^n X_i \\ \hat\mu_2 &= \frac1N\sum_{i=1}^n X_i^2 \end{align*} Uniendo todo esto, obtenemos una cantidad pivotante aproximada $$\frac{Y-E(Y)}{\sqrt{\text{Var}(Y)}} = \frac{Y-\theta}{\sqrt{\theta\frac{\mu_2}{\mu_1}}} \approx \frac{Y-\theta}{\sqrt{\theta\frac{\hat\mu_2}{\hat\mu_1}}}$$ que tendrá aproximadamente una distribución gaussiana estándar. Si dejamos que $z_*$ sea un valor crítico apropiado de esta distribución (por ejemplo $z_* \approx 1.96$ para un 95% de confianza), entonces la ecuación $$\left|\frac{Y-\theta}{\sqrt{\theta\frac{\hat\mu_2}{\hat\mu_1}}}\right| = z_*$$ conduce a una ecuación cuadrática en $\theta$ $$\theta^2 - (2Y + z_*^2\frac{\hat\mu_2}{\hat\mu_1})\theta + Y^2 = 0$$ Las soluciones de esta ecuación son los dos puntos finales de un intervalo de confianza aproximado para $\theta$ . A continuación se presenta una implementación de este procedimiento en R, junto con simulaciones para estimar el nivel de cobertura del intervalo de confianza (basado en $\lambda=1000$ y utilizando una distribución lognormal estándar para $X_i$ ):

compound_poisson_CI = function(x, conf.level){
  y = sum(x)
  m1 = mean(x)
  m2 = mean(x^2)
  z = qnorm((1+conf.level)/2)
  b = 2*y + z^2*m2/m1
  (b + c(-1,1)*sqrt(b^2-4*y^2))/2
}

# simulate compound Poisson random variables Y based on lognormal samples
set.seed(0)
lam = 1000
mu = lam*exp(0.5)   # true mean of Y
num_sims = 1000000
cnt = 0
for(i in 1:num_sims){
  n = rpois(1, lam)
  x = rlnorm(n)
  CI = compound_poisson_CI(x, 0.95)
  if(mu>=CI[1] && mu<=CI[2]){
    cnt = cnt + 1
  }
}
print(cnt/num_sims)

[1] 0.947843

Esto demuestra que la aproximación funciona bien, ya que el nivel de cobertura estimado 0,948 se acerca al nivel nominal del 95%. Para valores mayores de $\lambda$ (como en su escenario), esperaríamos que estuviera aún más cerca, aunque esto también puede depender de la asimetría de la distribución de los costes $X_i$ .

Aplicando este método a cada año, se podría construir una secuencia de intervalos de confianza, dando el coste total esperado para cada año, y estos se podrían juntar en un gráfico. Si el movimiento de los totales de un año a otro es grande en relación con la anchura de los intervalos de confianza, esto sería una prueba de los cambios subyacentes a lo largo del tiempo que no se deben simplemente a las variaciones del azar. Creo que mostrar los intervalos de confianza de esta manera hace que los resultados sean más fáciles de interpretar, en comparación con mostrar simplemente los totales de los años junto con un valor P de una prueba. Sin embargo, si está especialmente interesado en probar una hipótesis nula de ausencia de cambios en el coste total esperado a lo largo de los años, creo que podría hacerse construyendo una prueba de razón de verosimilitud basada en esta misma configuración del modelo; si lo desea, hágamelo saber y podría ampliarlo.

---

Editar : Así es como podemos construir una prueba de razón de verosimilitud generalizada utilizando esta configuración del modelo. Durante un período de $k$ años, tenemos una secuencia de costes observados $Y_1, \dots, Y_k$ . Podemos tratarlas como variables aleatorias independientes, cada una con una distribución normal con parámetros $E(Y_i)=\mu_i$ y $\text{Var}(Y_i) = \mu_ic_i$ , donde $\mu_i$ se desconoce, pero donde cada $c_i$ se trata como una constante conocida, calculada como la relación de los momentos de la muestra $c_i = \hat\mu_2/\hat\mu_1$ utilizando la muestra del año $i$ . Escribir $Y = (Y_1,\dots,Y_k)$ y $\mu = (\mu_1,\dots,\mu_k)$ la función de probabilidad se puede escribir $$f_Y(y|\mu) = \prod_{i=1}^k \frac1{\sqrt{2\pi\mu_ic_i}}\text{exp}\left(-\frac1{2\mu_ic_i}(y_i-\mu_i)^2\right)$$ Ahora nos interesa comparar dos conjuntos de modelos: 1) los modelos $\Theta_0$ donde $\mu$ se limita a ser una constante $\mu_i = b_0$ en todos los años, y 2) los modelos $\Theta_1$ donde $\mu$ sigue una tendencia lineal $\mu_i = b_0 + b_1i$ . En la prueba de razón de verosimilitud generalizada, determinamos la máxima verosimilitud para cada uno de estos dos conjuntos y luego consideramos la razón: $$\Lambda = \frac{\underset{\mu\in\Theta_0}{\max} f_Y(y|\mu)}{\underset{\mu\in\Theta_1}{\max} f_Y(y|\mu)}$$ La estadística $-2\log\Lambda$ se llama desviación y bajo las condiciones del teorema de Wilks, tiene aproximadamente un $\chi^2(1)$ distribución. En general, sería una $\chi^2(n)$ distribución, donde $n$ es el número de parámetros libres adicionales en $\Theta_1$ en comparación con $\Theta_0$ . A partir de ahí, podemos construir una prueba. Aquí es conveniente trabajar directamente con $-2\log f_Y(y|\mu)$ en lugar de la probabilidad $f_Y(y|\mu)$ por lo que calculamos $$-2\log f_Y(y|\mu) = \sum_{i=1}^k\left(\log(2\pi c_i) + \log(\mu_i) + \frac{(y_i-\mu_i)^2}{\mu_ic_i}\right)$$ Podemos utilizar la optimización numérica para encontrar las máximas probabilidades, utilizando un ajuste de mínimos cuadrados ordinarios para darnos nuestra conjetura inicial, que por lo general resulta estar ya muy cerca. Aquí se explica cómo implementarlo en R:

library(trustOptim)
library(tidyverse)

log_like = function(mu, y, C){
  # Inputs:
  # mu = vector of n values giving parameter values (expected total cost each year)
  # y = vector of n values giving observed data (actual total cost each year)
  # C = vector of n values giving ratio of variance to mean (approximated as ratio of 1st and 2nd moments)
  # Output: -2*log-likelihood, up to constant difference
  sum(log(mu) + (y - mu)^2/(mu*C))
}

log_like_gradient = function(mu, y, C){
  # Output: gradient of -2*log-likelihood
  1/mu + (1-(y/mu)^2)/C
}

model_mle = function(formula, C, ...){
  L = lm(formula, ...) # Apply ordinary least squares to get initial guess
  A = model.matrix(L)  # ... and to extract the model matrix
  n = length(y)
  f = function(p){ log_like(A %*% p, y, C) }
  f_gr = function(p){ t(A) %*% log_like_gradient(A %*% p, y, C) }
  trust.optim(coef(L), f, f_gr, method = 'SR1', control=list(report.level = 0))
}

deviance_test_pvalue = function(model_null, model_alt){
  deviance = model_null$fval - model_alt$fval
  pchisq(deviance, length(model_alt$gradient) - length(model_null$gradient), lower.tail = F)
}

A partir de ahí, asumiendo que tienes tus datos cargados en vectores y y C con Y siendo los valores observados de $Y_1,\dots,Y_k$ y C siendo el $c_i$ se puede obtener el valor P de la siguiente manera:

model_null = model_mle(y ~ 1, C, data = data.frame(y=y))
model_alt = model_mle(y ~ i, C, data = data.frame(y=y, i=1:k))
P = deviance_test_pvalue(model_null, model_alt)

Para asegurarnos de que esto se implementa correctamente, podemos realizar simulaciones bajo un caso en el que la hipótesis nula es verdadera y asegurarnos de que los valores P resultantes tienen una distribución cercana a la uniforme:

set.seed(0)
lam = 1000
mu = lam*exp(0.5)   # true mean of Y
k = 10  # 10 years of simulated data per simulation
num_sims = 10000
P = numeric(num_sims)
for(i in 1:num_sims){
  y = numeric(k)
  C = numeric(k)
  for(j in 1:k){
    n = rpois(1, lam)
    x = rlnorm(n)
    y[j] = sum(x)
    C[j] = mean(x^2)/mean(x)
  }
  model_null = model_mle(y ~ 1, C, data = data.frame(y=y))
  model_alt = model_mle(y ~ i, C, data = data.frame(y=y, i=1:k))
  P[i] = deviance_test_pvalue(model_null, model_alt)
}
P = sort(P)
data = tibble(P = P, ecdf = 1:length(P)/length(P))
min_alpha = .01
z = qnorm(.975)
data %>%
  filter(P >= min_alpha) %>%
  ggplot() +
  geom_ribbon(aes(x = P, ymin = ecdf-z*sqrt(ecdf*(1-ecdf)/num_sims), ymax = ecdf+z*sqrt(ecdf*(1-ecdf)/num_sims)), fill = 'gray') +
  geom_step(aes(x = P, y = ecdf)) +
  annotate('segment', x = min_alpha, y = min_alpha, xend = 1, yend = 1, color = 'blue') +
  scale_x_log10(limits = c(min_alpha, 1)) +
  scale_y_log10() +
  labs(
    x = 'P-value',
    y = 'Empirical CDF'
  )

Aquí el tamaño de la prueba es indistinto de su valor nominal (como $\alpha=.01$ o $\alpha=.05$ ), al menos dentro del margen de error que tenemos basado en 10000 simulaciones. Por tanto, la prueba funciona correctamente en el caso de que la hipótesis nula sea cierta. Se podría utilizar un método similar para explorar la potencia de la prueba en diversas condiciones. Como comprobación rápida de la cordura, tomemos una simulación e introduzcamos una tendencia artificial obvia para asegurarnos de que la prueba se rechaza:

y = y+100*(1:k)
model_null = model_mle(y ~ 1, C, data = data.frame(y=y))
model_alt = model_mle(y ~ i, C, data = data.frame(y=y, i=1:k))
deviance_test_pvalue(model_null, model_alt)

[1] 1.844717e-26

Además, si en lugar de comprobar una tendencia lineal $\mu_i = b_0 + b_1i$ preferiría hacer la prueba con el modelo completo (similar a un ANOVA) donde todas las medias $\mu_i$ son parámetros libres, puede hacerlo sustituyendo la línea

model_alt = model_mle(y ~ i, C, data = data.frame(y=y, i=1:k))

con

model_alt = model_mle(y ~ i, C, data = data.frame(y=y, i=as.factor(1:k))