Let y(i) = (y(il),...,y(ini)) be a vector of n(i) responses from subject i, where i = l,...,m. The y(i) are assumed to follow the linear mixed effects model y(i) = X(i)alpha + Z(i)beta(i) + e(i), where alpha is a vector of p unknown fixed effects parameters; X-i is a known matrix of dimension (n(i) x p) linking alpha to y(i),beta(i) is a vector of q unobservable random effects; Z(i) is a known matrix of dimension (n(i) x q) linking beta(i) to Y-i, and e(i) is a vector of n(i) within individual random errors. We assume that the beta(i)'s and e(i)'s are independent and identically distributed with means beta and 0, respectively, and covariance matrices Sigma and V-i, respectively. This kind of model characterizes the common structure of repeated measures, growth curve, or longitudinal data. The maximun likelihood estimator for alpha and the restricted maximun likelihood estimators of the variance components under the normal model are usually determined by Newton-Raphson or EM algorithms. The computations of these algorithms are not simple. In this paper, we propose a computationally simple two-stage procedure, which unlike the likelihood based procedures, does not require the normality assumptions for the beta(i)'s and the e(i)'s. It is shown that the two-stage estimators for the parameters of the model are consistent estimators under the assumption that the fourth moments of the components of beta(i) and e(i) exist. A simulation study also shows that the two-stage estimators are not much inferior to the likelihood based estimators under the normal model.