I have a dataset with cost measurements for hospital admissions of ~25,000 patients over multiple years. Cost seems to be log-normal distributed and taking the logarithm results in an approximately normal distribution. Using this transformation, I have analysed differences in the average cost across years using one-way ANOVA.

Can I do a similar test on the total costs (i.e. sum) across years?

The only thing I could come up with is testing for a linear trend using a linear regression of the sums on year. However, this leaves me with only 8 instead of 25,000 data points, which intuitively strikes me as throwing away a lot of information.

  • $\begingroup$ Do you have exact date information of the incurred costs, or only the ? $\endgroup$ – kjetil b halvorsen Jul 30 '17 at 9:53
  • $\begingroup$ For each patient (and thus each cost observation) I have the date of admission and the length of stay. So far I have simply aggregated them by year of admission $\endgroup$ – prock Jul 30 '17 at 11:04
  • $\begingroup$ If you have differences across years, the shape of the marginal distribution isn't really relevant $\endgroup$ – Glen_b Jul 30 '17 at 12:45
  • $\begingroup$ Maybe then order the observations in time order and useCUSUM or some other method for dete ting structural change. Search this site. $\endgroup$ – kjetil b halvorsen Jul 30 '17 at 12:51
  • $\begingroup$ @Glen_b: I am not quite sure I am following. Could you do me the favour and elaborate what this would mean in this example? $\endgroup$ – prock Jul 30 '17 at 13:45

I would suggest modeling the total cost $Y$ in a given year as a compound Poisson random variable. Specifically, we can model the number of admissions as a Poisson random variable $N$ with (unknown) mean $\lambda$, and the costs of admissions as $X_1, \dots, X_N$, so that $$Y = \sum_{i=1}^N X_i$$ Here we assume $X_1, X_2, \dots$ are iid (and independent of $N$) with first and second moments \begin{align*} \mu_1 &= E(X_i)\\ \mu_2 &= E(X_i^2) \end{align*} If $\lambda$ is large (which is reasonable to assume), asymptotic theory shows that $Y$ is approximately Gaussian, with the following parameters: \begin{align*} E(Y) &= E(N)E(X) = \lambda\mu_1\\ \text{Var}(Y) &= E(N)E(X^2) = \lambda\mu_2 \end{align*} Therefore, $(Y - E(Y))/\sqrt{\text{Var}(Y)}$ has approximately a standard Gaussian distribution. We can use this fact to derive inferences about $E(Y)$. Let us write $\theta = E(Y)$ since this is the parameter of interest. Notice that we can express $\text{Var}(Y)$ in terms of $\theta$ as $$\text{Var}(Y) = \theta\frac{\mu_2}{\mu_1}$$ Here $\mu_1$ and $\mu_2$ may be closely approximated by the sample moments: \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*} Putting this together, we obtain an approximate pivotal quantity $$\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}}}$$ which will have approximately a standard Gaussian distribution. If we let $z_*$ be an appropriate critical value from this distribution (e.g., $z_* \approx 1.96$ for 95% confidence), then the equation $$\left|\frac{Y-\theta}{\sqrt{\theta\frac{\hat\mu_2}{\hat\mu_1}}}\right| = z_*$$ leads to a quadratic equation in $\theta$ $$\theta^2 - (2Y + z_*^2\frac{\hat\mu_2}{\hat\mu_1})\theta + Y^2 = 0$$ The solutions to this equation are the two endpoints of an approximate confidence interval for $\theta$. Here's an implementation of this procedure in R, together with simulations to estimate the coverage level of the confidence interval (based on $\lambda=1000$, and using a standard lognormal distribution for $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
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

[1] 0.947843

This shows that the approximation is working well, as the estimated coverage level 0.948 is close to the nominal level of 95%. For larger values of $\lambda$ (as in your scenario), we would expect it to be even closer, although this may also depend on the skewness of the distribution of the costs $X_i$.

Applying this method to each year, you could construct a sequence of confidence intervals, giving the expected total cost for each year, and these could be put together into a plot. If the movement of the totals from year to year is large relative to the width of the confidence intervals, then this would be evidence of underlying changes over time that are not simply due to chance variations. I think that displaying confidence intervals like this makes the results easier to interpret, compared to just showing the year totals together with a P-value from a test. However, if you are particularly interested in testing a null hypothesis of no change in expected total cost across years, I believe this could be done by constructing a likelihood ratio test based on this same model setup; if you want, let me know and I could expand on that.

Edit: Here is how we can construct a generalized likelihood ratio test using this model setup. Over a period of $k$ years, we have a sequence of observed costs $Y_1, \dots, Y_k$. We can treat these as independent random variables, each with a normal distribution having parameters $E(Y_i)=\mu_i$ and $\text{Var}(Y_i) = \mu_ic_i$, where $\mu_i$ is unknown but where each $c_i$ is treated as a known constant, computed as the ratio of the sample moments $c_i = \hat\mu_2/\hat\mu_1$ using the sample from year $i$. Writing $Y = (Y_1,\dots,Y_k)$ and $\mu = (\mu_1,\dots,\mu_k)$, the likelihood function can be written $$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)$$ Now we are interested in comparing two sets of models: 1) the models $\Theta_0$ where $\mu$ is constrained to be a constant $\mu_i = b_0$ across all years, and 2) the models $\Theta_1$ where $\mu$ follows a linear trend $\mu_i = b_0 + b_1i$. In the generalized likelihood ratio test, we determine the maximum likelihood for each of these two sets and then consider the ratio: $$\Lambda = \frac{\underset{\mu\in\Theta_0}{\max} f_Y(y|\mu)}{\underset{\mu\in\Theta_1}{\max} f_Y(y|\mu)}$$ The statistic $-2\log\Lambda$ is called the deviance, and under the conditions of Wilks' theorem, it has approximately a $\chi^2(1)$ distribution. In general it would be a $\chi^2(n)$ distribution, where $n$ is the number of additional free parameters in $\Theta_1$ compared to $\Theta_0$. Using this, we can construct a test. Here it is convenient to work directly with $-2\log f_Y(y|\mu)$ instead of the likelihood $f_Y(y|\mu)$, so we calculate $$-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)$$ We can use numerical optimization to find the maximum likelihoods, using an ordinary least squares fit to give us our initial guess, which usually turns out to already be very close. Here is how to implement it in R:


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)

From there, assuming you have your data loaded up into vectors y and C, with Y being the observed values of $Y_1,\dots,Y_k$, and C being the $c_i$, you can get the P-value as follows:

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)

To make sure that this is implemented correctly, we can run simulations under a case where the null hypothesis is true and make sure that the resulting P-values have close to a uniform distribution:

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() +
    x = 'P-value',
    y = 'Empirical CDF'

Here the size of the test is indistinguishable from its nominal value (such as $\alpha=.01$ or $\alpha=.05$), at least within the margin of error that we have based on 10000 simulations. So the test is working properly in a case where the null hypothesis is true. A similar method could be used to explore the power of the test under various conditions. As a quick sanity check, let's just take one simulation and throw in an artificial obvious trend to make sure the test rejects:

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

Also, if instead of testing against a linear trend $\mu_i = b_0 + b_1i$ you would prefer to test against the full model (similar to an ANOVA) where all the means $\mu_i$ are free parameters, you can do that by replacing the line

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


model_alt = model_mle(y ~ i, C, data = data.frame(y=y, i=as.factor(1:k))
  • $\begingroup$ I like the approach and thanks for the extensive explanation, it makes intuitive sense although I could not have formalised it that way. Regarding the p-value, I due agree that a plot including CIs would be a better way to report variation. However, changes in total cost are a side-note rather than the main part of the study and I am unsure whether there would be space for such a graph. Given the target audience of clinicians, a p-value in the main text with the plot in the supplement might be the way to go. How would I go about constructing a likelihood ratio test from your above model? $\endgroup$ – prock Aug 1 '17 at 10:54
  • $\begingroup$ I've added a section walking through how to implement the likelihood ratio test. $\endgroup$ – Brent Kerby Aug 2 '17 at 17:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.