Longitudinal imaging studies are essential to understanding the neural development of neuropsychiatric disorders, substance use disorders, and normal brain. Using appropriate image processing and statistical tools to analyze the imaging, behavioral, and clinical data is critical for optimally exploring and interpreting the findings from those imaging studies. However, the existing imaging processing and statistical methods for analyzing imaging longitudinal measures are primarily developed for cross-sectional neuroimaging studies. The simple use of these cross-sectional tools to longitudinal imaging studies will significantly decrease the statistical power of longitudinal studies in detecting subtle changes of imaging measures and the causal role of time-dependent covariate in disease process. The main objective of this paper is to develop longitudinal statistics toolbox, called LSTGEE, for the analysis of neuroimaging data from longitudinal studies. We develop generalized estimating equations for jointly modeling imaging measures with behavioral and clinical variables from longitudinal studies. We develop a test procedure based on a score test statistic and a resampling method to test linear hypotheses of unknown parameters, such as associations between brain structure and function and covariates of interest, such as IQ, age, gene, diagnostic groups, and severity of disease. We demonstrate the application of our statistical methods to the detection of the changes of the fractional anisotropy across time in a longitudinal neonate study. Particularly, our results demonstrate that the use of longitudinal statistics can dramatically increase the statistical power in detecting the changes of neuroimaging measures. The proposed approach can be applied to longitudinal data with multiple outcomes and accommodate incomplete and unbalanced data, i.e., subjects with different number of measurements.