We address statistical issues regarding modeling a collection of longitudinal response trajectories characterized by the presence of subject‐specific phase‐dependent effects and variation. To accommodate these two time‐varying individual characteristics, we employ a geometric stochastic differential equation for modeling based on a Brownian motion process and develop a two‐step paradigm for statistical analysis. This paradigm reverses the order of statistical inference in random effects model. We first extract individual information about phase‐dependent treatment effects and volatility parameters for all subjects. Then, we derive the association relationship between the parameters characterizing the individual longitudinal trajectories and the corresponding covariates by means of multiple regression analysis. The stochastic differential equation model and the two‐step paradigm together provide significant advantages both in modeling flexibility and in computational efficiency. The modeling flexibility is due to the easy adaptation of temporal change points for subject‐specific phase transition in treatment effects, whereas the computational efficiency benefits in part from the independent increment property of Brownian motion that avoids high‐dimensional integration. We demonstrate our modeling approach and statistical analysis on a real data set of longitudinal measurements of disease activity scores from a rheumatoid arthritis study. Copyright © 2011 John Wiley & Sons, Ltd.