[R-sig-ME] lme4 for a 4 level hierarchy : revisited

[Kim Pearce]

Hello everyone,

I am rephrasing the question I sent to the list yesterday for clarity and wondered if anyone has any further views?  

My example yesterday was hypothetical and I was questioning whether the syntax I provide below (including the term (1|region/hospital/ward)) is appropriate for a 4 level hierarchical model.  This term is an extension of that detailed for a 3 level hierarchy in Table 2 (page 7) of the paper "Fitting Linear Mixed-Effects Models Using lme4" by Douglas Bates, Martin Mächler, Ben Bolker and Steve Walker. In this document, the authors stipulate that a term of the form (1|g1/g2) would mean that a model would be generated with the intercept varying among (i) g1 and (ii) g2 within g1  (i.e. intercepts would be generated (i) for each level of g1 and (ii) for each level produced by the g1 x g2 combination).

I discuss a mixed effects logistic regression below with a 4 level hierarchy but the term (1|region/hospital/ward) would be appropriate too if my response had been of a continuous nature and we were instead using lmer() rather than glmer().

In my hypothetical example, we have nurses which are nested within wards and the wards are nested within hospitals and hospitals are nested within regions of the UK. Each nurse has a binary response recorded (Y). Nurses are at level 1, wards at level 2, hospitals at level 3 and regions at level 4.  Say each ward has one category of an ordered 3 category intervention applied where a higher ordered category of intervention was hypothesised to affect a nurse's binary response.

Hypothetically, I wish to generate a model with a fixed effect for "Intervention" and intercepts for (i) each region, (ii) each region x hospital combination and (iii) each region x hospital x ward combination.  Additionally, an intercept would be included in the fixed-effects portion of the model by default.  Would the following syntax be appropriate?

Note intervention appears in the model as a factor below.

m1<-glmer(Y~interventionf+(1|region/hospital/ward),data=data,family=binomial,control=glmerControl(optimizer="bobyqa"))

Various online sources have hinted that we could approach an example such as this by considering this as a 2 level hierarchy with nurses (level 1) nested within wards (level 2) and having variables which identify hospitals and regions merely added as predictors in the model - but this seems suspect to me as we would not be considering the full nested nature of the design.

Many thanks for your views.
Kind regards,
Kim

Dr Kim Pearce PhD, CStat, Fellow HEA
Newcastle University