Modern　neuroimaging　research　calls　for　statistical　methods　that　can　model　dynamic　relationships　between　a　functional　response　and　a　set　of　covariates.　Current　methods,　however,　remain　disparate　and　limited　in　their　ability　to　robustly　accommodate　real-world　data　and　integrate　smoothness　penalties.　In　this　work,　we　propose　an　M-estimation　framework　for　the　varying-coefficient　model　with　a　functional　response　that　encompasses　both　mean　and　quantile　regression.　To　accommodate　smoothness　regularization　and　circumvent　the　stringent　conditions　on　Fourier　coefficients　or　the　covariance　operator’s　eigenvalues　imposed　by　traditional　fixed-basis　representations,　we　assume　that　the　functional　coefficient　resides　in　a　reproducing　kernel　Hilbert　space.　We　show　that　our　proposed　estimator　is　minimax　rate　optimal　and　establish　convergence　properties　of　our　modified　alternating　direction　method　of　multipliers　algorithm.　We　further　propose　combining　a　weighted　M-estimator　and　a　copula　model　to　quantify　within-subject　spatial　dependence　to　improve　estimation　accuracy.　Simulation　studies　and　a　real-world　analysis　demonstrate　the　robustness　of　our　proposed　methods　to　outliers.　©　2024　ISI/BS.