The majority of neuroimaging inference focuses on hypothesis testing rather than effect estimation. With concerns about replicability, there is growing interest in reporting standardized effect sizes from neuroimaging group-level analyses. Confidence sets are a recently developed approach to perform inference for effect sizes in neuroimaging but are restricted to univariate effect sizes and cross-sectional data. Thus, existing methods exclude increasingly common multigroup or nonlinear longitudinal associations of biological brain measurements with inter- and intra-individual variations in diagnosis, development, or symptoms. We broadly generalize the confidence set approach by developing a method for arbitrary effect sizes in longitudinal studies. Our method involves robust estimation of the effect size image and spatial and temporal covariance function based on generalized estimating equations. We obtain more efficient effect size estimates by concurrently estimating the exchangeable working covariance and using a nonparametric bootstrap to determine the joint distribution of effect size across voxels used to construct confidence sets. These confidence sets identify regions of the image where the lower or upper simultaneous confidence interval is above or below a given threshold with high probability. We evaluate the coverage and simultaneous confidence interval width of the proposed procedures using realistic simulations and perform longitudinal analyses of aging and diagnostic differences of cortical thickness in Alzheimer's disease and diagnostic differences of resting-state hippocampal activity in psychosis. This comprehensive approach along with the visualization functions integrated into the pbj R package offers a robust tool for analyzing repeated neuroimaging measurements.