Temporal analysis of variation in random effects

Question

I am looking at patient data where the main outcome of interest is mortality within 30 days following hospitalisation with an emergency condition. I am working on data from 2003-2017, with approximately 100,000 observations per year. I had previously understood that this was to be a "simple" multilevel logistic analysis, with which I have some experience, but I have just learned that the main interest lies in the changes in hospital-level variation in the outcome during the period. So there is also a time series element to it, but the data are not repeat measurements at the patient level (in fact patients with multiple admissions are to be excluded).

The clinical lead said he thought that the data should be sliced into a data set for each year, and the same model run on each data set and then look at how hospital-level variation has changed. This strikes me as not ideal. On the other hand I don't know what is ideal. Any suggestions or links to similar studies would be most appreciated.

Answer 1

New answer: 2020 !

the main interest lies in the changes in hospital-level variation

You have 15 years of data, with over 100 hospitals and around 100,000 observations per year, so an average of around 1000 observations per hospital per year.

I think there is only one approach that will answer the research question, and that is to divide the data into subsets. How many subsets will depend on what frequency you wish to the changes in variation. Yearly would be one option, quarterly another or even monthly.

You would then fit a model on each subset of the data, with random intercepts for hospitals. Since your outcome is binary, this would be a generalised linear mixed model, with binomial family and logit link. You would use the same model formula for each subset. Then you simply extract the hospital level variation for each time period, ie. the variance or standard deviation of the hospital intercept and present that in whatever way is appropriate - graphically would be my choice. If you think that there is likely to be a seasonal component, then you probably want to use quarterly subsets.

Another approach (mentioned in another answer) is to model the time:hospital interaction as random. Here you would fit a model on the whole dataset but additionally with random intercepts for for the time variable interacted with hospital. Again you can choose whatever period for the time variable makes the most sense. You could also fit a model with just the hospital as random and use a likelihood ratio test to determine which model fits best. However, this will not answer the question of how the hospital level variation changes over time. The same applies to using a correlation structure, such as AR1, because this also will not say anything about changes in hospital-level variation.

Answer 2

Some more information is needed to figure out the best solution here, so I'm simply answering a number of scenarios with example R code.

Modeling the outcome

If the outcome is binary, use family = binomial(). If it is count data, use family = poisson() since it's a fixed time interval. You could also consider aggregating binary data to counts and then use Poisson.

I'll just assume binomial from here and on.

Modeling a fixed hospital:time interaction

If you have few hospitals (say, less than five), there will be very little information to infer the random effects hyperparameters. In that case, they may just be modeled as fixed:

fit_full = glm(status ~ time * hospital, data = df, family = binomial())

Then the resulting inference on the time:hospital interaction could be of interest.

Modeling random hospital:time interaction

One of the only practical implications of modeling a term as random is that it applies shrinkage. That is, data points far from the model's prediction are regarded as partially random fluctuations with the true value being closer to the mean. Read more here.

Using the same model as above, but allowing for random slopes for each hospital:

fit_full = lme4::glmer(status ~ time + (1 + time|hospital), family = binomial())

To test the random slope, you can do a Likelihood Ratio Test (LRT) by comparing to a (nested) model that does not contain this term:

fit_null = lme4::glmer(status ~ time + (1|hospital), family = binomial())
summary(anova(fit_full, fit_null))

Personally, I have a preference for Bayesian inference and you could use the brms package which is much like glmer. As a quick fix, you can also compute a BIC-based Bayes Factor (cf. Wagenmakers et al. (2007):

exp((BIC(fit_full) - BIC(fit_null))/2)

Modeling time series

I know of no other packages that can do the above and model some autocorrelation than brms (and perhaps nlmer::lme). brms may take hours to fit, though. For AR(1), it would be something like:

fit = brms::brm(status ~ time + (1 + time|hospital), data = df, family = bernoulli(), autocor = cor_ar(~1, p = 1))

If your dates only come in integers (2013, 2014, 2015, 2016, 2017), then there is may be too little information to estimate autoregressive coefficient(s) and you may consider leaving it out. You do have a lot of data, so this may now be necessary. Your time variable would need a finer resolution for an autoregressive model to be identifiable.

Answer 3

You've really been thrown in the deep end ! It doesn't seem like a time series problem, but does seem though it could naturally be modelled as a multilevel regression. As a first step (after the usual data exploration etc. of course) I would probably fit a generalised linear mixed effects model. To include a time compnent, you could then add a time variate (1=2003, 2=2004 etc). There are probably better ways to build time into it - I imagine others will have a better odea on that.