# Temporal analysis of variation in random effects

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.

• I am note understand well what are "changes in hospital-level variation in the outcome during the period." What kind of data do you have? Year, hopital id, patient id, died within 30 days? Is that it? And you want to check whether during the 6 years of data collecting, the mortality within 30 days changed in some of the hospitals? – Elvis Jan 18 '12 at 13:40
• @PSellaz: How about an update? Did you sink or swim? – Aaron left Stack Overflow Sep 13 '12 at 19:44
• @Aaron sorry, I only just saw your comment. In the end, I just ran separate random intercepts models for each year of data, using glmer. I did bootstrapping to compute confidence intervals for the random effects and these confidence overlapped to a very high degree from year to year, with no apparent trend, from which I concluded there was no change over the period. The fixed effects were also very stable from year to year. – P Sellaz Mar 6 '13 at 12:33
• Do you still have the same dataset, or has it been updated since you first posted this (over 8 years ago !!!) ? – Robert Long Jan 17 at 11:05
• @RobertLong yes the dataset now covers 2003 - 2017 inclusive – P Sellaz Jan 17 at 12:11

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.

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.

• Thank you. What about using a 3 level hierarchical approach, with patients within years within hospitals. Would this make sense ? – P Sellaz Dec 18 '11 at 23:08
• Yes, maybe. How many hospitals are patients are there ? – Robert Long Dec 19 '11 at 7:47
• Around 100 hospitals, on average around 500 patients per hospital (but of course it's not balanced) – P Sellaz Dec 19 '11 at 19:44
• I did try the 3-level approach - patients within years within hospitals, but I couldn't get it to work. I think this is a multiple membership model - patients "belong" to a year, but each year "belongs" to every hospital. I didn't know how to model this... – P Sellaz Mar 6 '13 at 12:32

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.