The question

I have a dataset which I think requires a multivariate multilevel analysis. I am unsure both of the appropriate model and of how to fit it with R. I have come up with a tentative model, but my understanding of the math is so superficial that I cannot tell whether my analysis is "right" or whether it includes blatant errors. I would appreciate any insight on the model design or the model specification in R.

The study design

The question is whether the architectural design of a clinic will influence the outcome of a pathology for permanent residents in this clinic. We have collected data on 13 symptoms for 8 patients per clinic in 21 clinics.

There is a patient-level IV medication and a clinic-level IV architecture. All variables are continuous-ish.

The 13 symptoms are correlated +.20 on average, which I think indicates a multivariate multilevel analysis is appropriate.

The data

To run the multivariate analysis with nlme I have standardized my DVs, stacked these 13 DVs in a single column, and added a categorical dummy variable to flag which row corresponds to which symptom.

It looks like this:

 Clinic Patient Symptom    Score    Medication  Architecture
 1            1   EP1      0.12         1               3.2  
 1            1   EP2      0.11         1               3.2  
 1            1   EP3      0.13         1               3.2  
 1            2   EP1      0.56         4               3.2  
 1            2   EP2      0.67         4               3.2  
 1            2   EP3      0.23         4               3.2  
 2            3   EP1      0.22         3               5.1  
 2            3   EP2      0.25         3               5.1  
 2            3   EP3      0.14         3               5.1  
 2            4   EP1      0.78         6               5.1  
 2            4   EP2      0.89         6               5.1  
 2            4   EP3      0.11         6               5.1  

The model design

  • To run the analysis as multivariate, I use both symptom and symptom:architecture as IVs and I remove the intercept in both the fixed and random parts of the model. I do not include the main effect of architecture as an IV.
  • The effect of medication should be the same within all clinics, so there is no random effect for this variable.
  • I do not want to constrain equality between the effect of architecture on the different symptoms.
  • Due to the multivariate nature of the analysis, I expect the residuals to be correlated, with different correlations between the 13 different symptoms; therefore I specify the covariance structure of residuals as corSymm (non-zero but unstructured, if I get this correctly).
  • I also expect heteroscedasticity between the different symptoms (there should be more variance on certain symptoms), so I add the option weights as (~ 1|symptoms).

The end result

This is the model I come up with:

model1 = lme(fixed = Score ~ symptom + medication:symptom + architecture:symptom + medication:architecture:symptom - 1,  
+ random = ~ symptom - 1 | patient/clinic,  
+ correlation = corSymm,  
+ weights=varIdent(form= ~ 1|symptoms)  
+ method = "ML")  

In order to test the effect of the architectural variables, I would then compare this model to the following constrained model, dropping all the terms related to architecture:

model2 = lme(fixed = Score ~ symptom + medication:symptom - 1,
+ random = ~ symptom - 1 | patient/clinic,
+ correlation = corSymm,
+ weights=varIdent(form= ~ 1|symptoms)
+ method = "ML")

I would then run this comparison with the command anova(model1, model2) and compare the log-likelihood of the two models.

Overall, do these model design and r specification look correct to you? Thank you so much for your help!


I am curious why you have mentioned it as a multivariate analysis, when 'Score' appears to be the only outcome in your model.

To answer your question, there are no right or wrong models. The model fit can be compared using log-likelihood. Personally, I feel that adding one term at a time, and checking model fit is a better strategy than comparing model fit of a 'full model' with terms that you desire to test, with a comparatively basic model. The effect of variable of interest (such as symptoms) can be checked one at a time.

Lastly, collinear variables will still be problematic in your above formulation.


Woah, I wasn't actually expecting an answer to that one anymore. Thank you for taking the time to reflect on this.

Of course, four years later the project is long buried, but here are some things I've learned along the way in case someone else comes across these issues:

  • This was indeed a multivariate analysis since Score represented the score on various symptoms, as indexed by the column "symptom" in the example. The patients completed several different measures for several different symptoms, which were then collapsed in a single "score" column to allow for simultaneous testing of the different symptoms. In essence, the multilevel framework is used to yield a multivariate analysis based on three measurement levels: type of measure within patient within clinic. This solution is suggested by Goldstein (multilevel statistical models, chapter 4). I never came across a published example of someone doing it that way; even if correct, it is counter-intuitive enough that it might be difficult to publish.
  • I'm still not sure why some sources recommend removing the intercept in the multivariate multilevel case. I think it is due to the way some of them use dummy coding to recode categorical variables. With more current software, where categorical variables can be imported as is, I would probably leave the intercept in the model.
  • I definitely agree with your strategy of progressively building more complex models, at least for initial model specification and sanity checking. In psychology though you'll often want to run the actual test of the hypothesis by comparing the most complex model with the most-complex-without-the-IV-model, so as to return to a familiar Type-III sum-of-squares framework. This complex model would have been a bit dangerous in the present case, however, given that some of the control variables were heavily unbalanced in terms of sample size - which would probably have messed up the test of the main effects of the predictors of interest.
  • If I had to redo this analysis now, I would definitely make the model less complex by removing the unnecessary gizmos that are only there to improve model fit: heteroscedasticity and correlated residuals. While I think they are theoretically useful in that case, they tend to impede model convergence with such a small sample size, and with multiple measures (which was the case here with 13 symptoms) they make computation time prohibitive even on very powerful computers.

Thanks for your feedback!


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.