We present a multivariate functional mixed efects 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 efects functions accounting for the dependence among bilateral measurements. The model is ftted by applying multivariate functional principal component analysis (mv-FPCA) and modelling the mv-FPCA scores using scalar linear mixed efects models. Simulation and bootstrap approaches are introduced to construct simultaneous confdence bands for the fxed efects 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 scientifc application, we observe a statistically signifcant efect of running speed on both joints. We observe strong within-subject correlations, refecting the highly idiosyncratic nature of running technique. Our approach is applicable to modelling multiple streams of smooth biomechanical data collected in complex experimental designs.