We present a multivariate functional mixed effects model for kinematic data from a large number of recreational runners. The runners’ sagittal plane hip and knee angles are modelled jointly as a bivariate function with random effects functions accounting for the dependence among bilateral measurements. The model is fitted by applying multivariate functional principal component analysis (mv-FPCA) and modelling the mv-FPCA scores using scalar linear mixed effects models. Simulation and bootstrap approaches are introduced to construct simultaneous confidence bands for the fixed effects functions, and covariance functions are reconstructed to summarise the variability structure in the data and thoroughly investigate the suitability of the proposed model. In our scientific application, we observe a statistically significant effect of running speed on both joints. We observe strong within-subject correlations, reflecting the highly idiosyncratic nature of running technique. Our approach is applicable to modelling multiple streams of smooth biomechanical data collected in complex experimental designs