<p>Gaussian process regression is a common method for flexible yet fully probabilistic nonlinear regression. A frequent obstacle is its computational complexity, which scales poorly with the number of observations. The problem intensifies when Gaussian process models are applied simultaneously to multiple functions. We consider a multilevel Gaussian process regression model in which a common mean function and subject-specific deviations are jointly modeled as latent Gaussian processes. We derive exact, analytic, and computationally efficient expressions for the log-likelihood and the conditional posterior distributions when observations are sampled on either a completely or partially regular grid. Without using approximations, these expressions enable us to fit the model to large data sets that are currently computationally inaccessible with a standard implementation. We show through a simulation study that our analytic expressions are several orders of magnitude faster than a standard implementation, and we provide an implementation in the probabilistic programming language Stan.</p>