Robust estimation of multivariate reference bands

I have a subjectivly healthy population of approx 1200 individuals with three measurements on the continous scale, we can call them y1, y2 and y3. All of these are strongly related to age in a non-linear fashion and I would like to create age-dependent uni- and multivariate 95% reference bands. The variables y1, y2 and y3 are not normally distributed and, at least one, are almost impossible to tranform to normality (the problem really is that the residuals becomes non-normal or heteroskedastic no matter what I try). I'm guessing the reason for this are the presence of outliers. I can't really exclude the outliers from the material as I have no evidence of data entry error or disease. Robust methodology seems like the path to take.

Would the following strategy be "statistically" sound:

1. Perform robust regression (mmregress in stata with appropriate fp terms) on untransformed data (three regressions with y1, y2 and y3 as dependent variable and age as independent variable). Find outliers by plotting robust mahalanobi's distance of residuals (mcd in stata) against expected chi2-distribution and visually identify ouliers (mark them, but keep them).
2. Transform y1,y2,y3 to approx multivariate normal distribution (mboxcox in stata).
3. Perform robust regression with transformed y1,y2,y3 with appropriate fp terms of age. Calculate residuals for the three regressions (observed minus predicted value).
4. Check assumption: are the residuals univariatly normal distributed for all ages? Are they multivariatly normally distributed? (if not so for the full populations maybe it is sufficient if the assumptions hold for the outlier-free sample)
5. Check assumption: are the variances and covariances constant (lowess of squared residuals and residual products) with regard to age? (if not so for the full populations maybe it is sufficient of the assumptions hold for the outlier-free sample)
6. Use mcd to obtain the robust variance-covariance matrix of the residuals.
7. Take sqrt of diagonal terms of the variance-covariance-matrix and use as a "robust" SD for the univariate reference intervals (predicted +/- 1.96*SD) and then transform back to original scale.
8. Calculate the 95% multivariate reference region by using (y-py)'C-1(y-py)=7.815 where y and py are 3 element column vectors with observed and predicted values respectivly, C-1 the inverted "robust" variance covariance matrix and 7.815 from the chi2 distribution for a 95% reference regions and 3 df. Transform back to original scale.

Main hesitations: Is it ok to relax the assumptions by just checking for them in the outlier free subsample? Is it ok to use the robustly estimated variance-covariance matrix for calculation of univariate SD and reference intervals? Any thoughts? Objections? Suggestions?

• Aren't you asking for bands rather than intervals, since (if I understand correctly) even given a fixed $y_i$, you want a different interval for every age? – Kodiologist Jun 5 '16 at 15:42
• You are right. Univariatly I want separate equations for upper and lower 95% age-related reference centiles for y1, y2 and y3 (for instance: [lower limit of y1]=retransform([predicted transform(y1) for age]-1.96*sd)). And on top of that a age-related 95% multivariate reference region. To succeed with this I need transformations that makes the residuals normally distributed and the variance and covariance of the residuals independent of age). Then I can use the same SDs and the same variance-covariance matrix for all ages and only regress the predicted mean. – user301448 Jun 5 '16 at 17:53