In large multilevel studies effects of interest are often evaluated for a number of more or less related outcomes. For instance, the present work was motivated by the multiplicity issues that arose in the analysis of a cluster-randomized, crossover intervention study evaluating the health benefits of a school meal programme. We propose a novel and versatile framework for simultaneous inference on parameters estimated from linear mixed models that were fitted separately for several outcomes from the same study, but did not necessarily contain the same fixed or random effects. By combining asymptotic representations of parameter estimates from separate model fits we could derive the joint asymptotic normal distribution for all parameter estimates of interest for all outcomes considered. This result enabled the construction of simultaneous confidence intervals and calculation of adjusted p-values. For sample sizes of practical relevance we studied simultaneous coverage through simulation, which showed that the approach achieved acceptable coverage probabilities even for small sample sizes (10 clusters) and for 2–16 outcomes. The approach also compared favourably with a joint modelling approach. We also analysed data with 17 outcomes from the motivating study, resulting in adjusted p-values that were appreciably less conservative than Bonferroni adjustment.