Johannes Fendel Yes, based on your description a multilevel model would work.

The number of Level-2 units (therapist) would be very small (n=8)--probably too small for a meaningful multilevel analysis. Cohen, Cohen, West, and Aiken (2003, Applied multiple regression/correlation analysis for the behavioral sciences) recommend using dummy coding for designs with few groups instead of multilevel modeling. That way, you could study whether there is a therapist effect (differences between therapists) in addition to a treatment effect using regression or AN(C)OVA rather than multilevel analysis. Either way, your statistical power for showing differences is probably limited given the rather small sample size.

Christian Geiser I don't see why 34 participants clustered within 8 therapists is problematic for a mixed effects model. Since the design is both unbalanced and contains both fixed (treatment) and random (participant, therapist) factors, a linear regression or ANCOVA is not going to be able to examine both in a single analysis - nor handle the complex clustered design. In my opinion, this clustered design (treatment nested within participant nested within therapist) calls for a more flexible model to estimate the slopes.

I agree with Blain. The low number of therapists will result in a low-accuracy estimate for the variance between therapists, but this is not of interest for the research question here.

Blaine Tomkins Jochen Wilhelm Check out the following article by McNeish and Stapleton (2016):

Daniel McNeish & Laura M. Stapleton (2016). Modeling Clustered Data with Very Few Clusters, Multivariate Behavioral Research, 51:4, 495-518, DOI: 10.1080/00273171.2016.1167008

Article Modeling Clustered Data with Very Few Clusters

Their simulations revealed that the fixed effects (dummy coding) approach performed better in terms of power than did multilevel modeling in designs with a small number of clusters like in the present example. I believe that power does matter here. McNeish and Stapleton also discuss other advantages of the fixed effects approach to representing the clusters.

Dear Christian Geiser , thank you for the article. It is interesting, but I must admit that I am not able to fully interpret the results presented therein (particularily the table about power; firstly I think that investigating the power through simulation with such tiny effects in a quite complex model is not much more than digging in the noise, and secondly I have not really understood the important results tables (9+10). I also don't understand why the results in table 11 do not nealy match the parameters in table 3 (as I understood were these parameters chose to approximately match findings from the real data). Maybe some of my confusion will dissolve after reading the article more often and more carefully, but this takes too much time for now.

The results in table 11 look very similar for all methods. I did not recognize gross differences that I had expected from the main message of the paper. On the one hand the authos say that FEM provides more power for study desigs as the one analyzed, on the other hand the FEM was the only one that failed to reject the hypothesis of zero treatment effect.

Christian Geiser power is certainly important, not just here but always. However, statistical power only reflects the probability of rejecting the null. Just because model A rejects H0 with fewer participants and/or fewer assumptions than model B does not imply that model A more accurately reflects the processes that generated the data.

Perhaps the reason a model containing solely fixed factors rejects the null more easily (i.e. has more power) in a particular situation is because it does not take into account the additional variability introduced by random factors. Sure, omitting this additional variability will give us more power, but at the expense of missing the real variability that likely exists across different therapists, patients, etc.,.

Not only is power important, but also error of the estimates. Looking at Table 11 - the mixed models actually have slightly lower standard errors than the FEM. For example:

ELL: ML (2.37), REML (2.38), FEM (2.53)

Treat: ML(2.37), REML (2.54), FEM (2.82)

ELL x Treat : ML (3.19), REML (3.25), FEM (3.38)

What I don't see in this article is a direct statistical comparison among the alternative models. As Jochen Wilhelm said above, there doesn't seem to be much of a difference in the estimates. My hunch is if you directly compared these models FEM would not show a statistically better fit. Thus, I'm not convinced the authors have supported their claim that FEM is the best approach in this situation.

Jochen Wilhelm Blaine Tomkins I am tagging Daniel McNeish in case he is active on Researchgate and wants to chime in.

My thought process here is the following: Yes, no doubt you can do multilevel modeling in this case. However, why would you want to go through the complexities of running a multilevel model with n = 34 at Level 1 and n = 8 at Level 2 when the sole purpose of the analysis is to account for clustering (non-independence of observations) and when you can accomplish this with a much simpler method (single-level regression with dummy coding of Level-2 units), for which, on top of everything else, there is evidence in the literature (as per McNeish & Stapleton's 2006 study) that this approach may show better performance than more complex methods when the number of clusters is very small as in the present example?

I am not surprised that most estimates in the simulation and application are very close across approaches. On the contrary, I would expect this and I'd be rather surprised if they were very different. These are different (but statistically near equivalent) ways of accounting for the same process. They should produce similar estimates (although with potential differences in the degree of parameter estimate and standard error bias). For that reason, I also don't think that a comparison of fit statistics would be relevant or useful here or that "better fit" of one approach relative to another would tell us much (in other words, I'd expect most approaches to trivially result in more or less the same fit). I'm assuming that that's the reason why fit statistics were not considered in the article. However, I may be missing something and perhaps Daniel McNeish could speak to this issue.

Christian Geiser With respect to performance it depends on whether the goal is fitting or prediction. To quote Richard McElreath, "fitting is easy; prediction is hard. Future data will not be exactly like past data, and so any model that is unaware of this fact tends to make worse predictions than it could"

Presumably, the study above is intended to generalize beyond this group of therapists and patients. As such, our model should account for this random variability.

I can't phrase it any better than McElreath, so permit me to crib his words once more:

"Multilevel regression deserves to be the default form of regression. Papers that do not use multilevel models should have to justify not using a multilevel model approach. Certainly some data and contexts do not need the multilevel treatment. But most contemporary studies in the social and natural sciences, whether experimental or not, would benefit from it. Perhaps the most important reason is that even well-controlled treatments interact with unmeasured aspects of the individuals, groups, or populations studied. This leads to variation in treatment effects, in which individuals or groups vary in how they respond to the same circumstance. Multilevel models attempt to quantify the extent of this variation, as well as identify which units in the data responded in which ways" (Ibid.)

Blaine Tomkins , is the comparison McElreath between multilevel versus ignoring the clustering, or as Christian Geiser (and if he joins in likely Daniel McNeish ) between random versus fixed effect with small numbers of clusters? Is that from McElreath's book?

Daniel Wright Yes, the above is from McElreath's "Statistical Rethinking". I'd say his view is moreso concerned with the first point: ignoring random variation present in clustered designs altogether. However, I think his argument holds regardless of the number or levels of clusters. There will always be some amount of random variability in a clustered design and a mixed model will capture more information than an OLS regression.

Whether that additional information is useful or meaningful is, of course, situation dependent.

Blaine Tomkins , Daniel Wright , Christian Geiser, Jochen Wilhelm, Thank you all!

In R, using lmer with "individual" as random intercept gives identical p-values as lm with "individual" as additional fixed factor. This is independent of the actual distribution of inter-individual differences. This works for as few as two clusters!

Yes, lmer() occasionally runs into singularity issues, in which cases the resulting p-value is lower than that of the corresponding fixed-effect model. And yes, if the inter-individual variance is not interesting, the fixed-effects model may be the simpler solution to fit - but it is more difficult to get "population predictions".

Sounds like a multi-level analysis would be perfect, especially considering the unbalanced design.

Thank you all for your thoughts! As we saw, the data structure is a bit more complicated. I propose two options for a multilevel analysis and I would be grateful for your thoughts on that:

Option 1:

· Level 3: Therapists (n = 7)

· Level 2: Individuals (n =34)

· Level 1: Condition (Treatment A and Treatment B; number of observation = 68)

As described above, Level 3 (=therapists) is unbalanced: First, the therapists did not perform the same number of treatments overall (i.e., Therapist 1 performed 25 treatments, Therapist 2 performed 10 treatments, and so on…). Second, each therapist did not perform Treatment A or Treatment B with equal frequency. Moreover, one therapist never performed both treatments for one individual. In other words, the individuals (Level 2) are nested in two different therapists (Level 3). I tried the following R code with the nlme-package:

Model_1

I don't get your experimental design. "the individuals (Level 2) are nested in two different therapists (Level 3)." - that means that patients are assigned to (exactly) one therapist, so for instance patient #3 for therapist #5 is a different patient than patient #3 for any other therapist ?

Maybe you could attach the R data.frame, if you like with altered or random values for the response variable.

Jochen Wilhelm

Thank you! Each individual is assigned to two different therapists: one therapist for Treatment A and one Therapist for Treamtent B - so for instance, individual #3 for therapist #1 for Treatment A and therapist #4 for Treatment B. Please find a .csv attached!

I think the key problem is that you try to fit a random slope for "condition". This does not work in your design.

What works is:

lme(outcome ~ condition, random = ~1|therapist/individual)

I don't know how the original data look like, but for the (fake?) data you provided the outcome is very bimodal with values either close to 0 or close to 1. It might be more appropriate to use a logit transformation, or a beta-regression model (usually, a quasi-bionmial model would do as well, but I don't think that they are available as mixed-effect models).