Despite interest in the joint modeling of multiple functional responses such as diffusion properties in neuroimaging, robust statistical methods appropriate for this task are lacking. To address this need, we propose a varying coefficient quantile regression model able to handle bivariate functional responses. Our work supports innovative insights into biomedical data by modeling the joint distribution of functional variables over their domains and across clinical covariates. We propose an estimation procedure based on the alternating direction method of multipliers and propagation separation algorithms to estimate varying coefficients using a B-spline basis and an |$L_2$| smoothness penalty that encourages interpretability. A simulation study and an application to a real-world neurodevelopmental data set demonstrates the performance of our model and the insights provided by modeling functional fractional anisotropy and mean diffusivity jointly and their association with gestational age and sex.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (