This paper considers the problem of fitting regression models to historical fleet data with mixed effects, which arises in the context of statistical monitoring of data from a fleet (population) of similar units. A fleet-wide extension of the multivariable statistical process control approach is used to monitor for three different types of faults: a performance anomaly, a performance shift, and an anomalous unit. Our formulation requires the solution of a least-squares problem with very large numbers of both regressors (variables) and data measurements. For problems of interest, this least-squares problem cannot be solved using standard methods. We propose a method for solving the problem that is scalable to extremely large datasets, even ones that do not fit in to the memory of a single computer system. Our method can be parallelized, but also works serially on a single processor. This approach is demonstrated in a simulated example for monitoring a fleet of aircraft from historical cruise flight data.