Rfit: Rank-based Estimation for Linear Models

In the nineteen seventies, Jurečková and Jaeckel proposed rank estimation for linear models. Since that time, several authors have developed inference and diagnostic methods for these estimators. These rank-based estimators and their associated inference are highly efficient and are robust to outliers in response space. The methods include estimation of standard errors, tests of general linear hypotheses, confidence intervals, diagnostic procedures including stu-dentized residuals, and measures of influential cases. We have developed an R package, Rfit, for computing of these robust procedures. In this paper we highlight the main features of the package. The package uses standard linear model syntax and includes many of the main inference and diagnostic functions.


Introduction
Rank-based estimators were developed as a robust, nonparametric alternative to traditional likelihood or least squares estimators.Rank-based regression was first introduced by Jurečková (1971) and Jaeckel (1972).McKean and Hettmansperger (1978) developed a Newton step algorithm that led to feasible computation of these rank-based estimates.Since then a complete rank-based inference for linear models has been developed that is based on rank-based estimation analogous to the way that traditional analysis is based on least squares (LS) estimation; see Chapters 3-5 of the monograph by Hettmansperger and McKean (2011) and Chapter 9 of Hollander and Wolfe (1999).Furthermore, robust diagnostic procedures have been developed with which to ascertain quality of fit and to locate outliers in the data; see McKean and Sheather (2009) for a recent discussion.Kloke et al. (2009) extended this rank-based inference to mixed models.Thus rank-based analysis is a complete analysis analogous to the traditional LS analysis for general linear models.This rank-based analysis generalizes Wilcoxon procedures for simple location models and, further, it inherits the same high efficiency that these simple nonparametric procedures possess.In addition, weighted versions can have high breakdown (up to 50%) in factor space (Chang et al., 1999).In this paper, we discuss the Rfit package that we have developed for rank-based (R) estimation and inference for linear models.We illustrate its use on examples from simple regression to k-way factorial designs.
The geometry of rank-based estimation is similar to that of LS.In rank-based regression, however, we replace Euclidean distance with another measure of distance which we refer to as Jaeckel's dispersion function; see Hettmansperger and McKean (2011) for details.For a brief overview see McKean (2004).
Jaeckel's dispersion function depends on the choice of a score function.
As discussed in Hettmansperger and McKean (2011), the rank-based fit and associated analysis can be optimized by a prudent choice of scores.If the form of the error distribution is known and the associated scores are used, then the the analysis is fully efficient.In Rfit we have included a library of score functions.The default option is to use Wilcoxon (linear) scores, however it is straightforward to create user-defined score functions.We discuss score functions further in a later section.
Others have developed software for rank-based estimation.Kapenga et al. (1995) developed a Fortran package and Crimin et al. (2008) developed a web interface (cgi) with Perl for this Fortran program.Terpstra and McKean (2005) developed a set of R functions to compute weighted Wilcoxon (WW) estimates including the high breakdown point rank-based (HBR) estimate proposed by Chang et al. (1999).See McKean et al. (2009) for a recent review.Rfit differs from the WW estimates in that its estimation algorithms are available for general scores and it uses a standard linear models interface.
The package Rfit allows the user to implement rank-based estimation and inference described in Chapters 3-5 of Hettmansperger and McKean (2011) and Chapter 9 of Hollander and Wolfe (1999).There are other robust packages in R. For example, the R function rlm of the R package MASS (Venables and Ripley, 2002) computes M estimates for linear models based on the ψ functions of Huber, Hampel, and Tukey (bisquare).The CRAN Task View on robust statistical methods offers robust procedures for linear and nonlinear models including methods based on M, M-S, and MM estimators.These procedures, though, do not obtain rank-based estimates and associated inference as do the procedures in Rfit.
Rfit uses standard linear model syntax so that those familiar with traditional parametric analysis can easily begin running robust analyses.In this paper, discussion of the assumptions are kept to a minimum and we refer the interested reader to the literature.All data sets used in demonstrating Rfit are included in the package.
The rest of the paper is organized as follows.

Rank-regression
As with least squares, the goal of rank-based regression is to estimate the vector of coefficients, β, of a general linear model of the form: where y i is the response variable, x i is the vector of explanatory variables, α is the intercept parameter, and e i is the error term.We assume that the errors are iid with probability density function (pdf) f (t).
For convenience, we write (1) in matrix notation as follows where y = [y 1 , . . ., y n ] T is the n × 1 vector of responses, X = [x 1 , . . ., x n ] T is the n × p design matrix, and e = [e 1 , . . ., e n ] T is the n × 1 vector of error terms.The only assumption on the distribution of the errors is that it is continuous; in that sense the model is general.Recall that the least squares estimator is the minimizor of Euclidean distance between y and ŷLS = X βLS .To obtain the R estimator, a different measure of distance is used that is based on the dispersion function of Jaeckel (1972).Jaeckel's dispersion function is given by where • ϕ is a pseudo-norm defined as where R denotes rank, a(t) = ϕ t n+1 , and ϕ is a non-decreasing, square-integrable score function defined on the interval (0, 1).Assume without loss of generality that it is standardized, so that ϕ(u) du = 0 and ϕ 2 (u) du = 1.Score functions are discussed further in a later section.
The R estimator of β is defined as This estimator is a highly efficient estimator which is robust in the Y-space.A weighted version can attain 50% breakdown in the X-space at the expense of a loss in efficiency (Chang et al., 1999).

Inference
Under assumptions outlined in the previous section, it can be shown that the solution to (4) is consistent and asymptotically normal (Hettmansperger and McKean, 2011).We summarize this result as follows: βϕ ∼N β, τ 2 ϕ (X T X) −1 where τ ϕ is a scale parameter which depends on f and the score function ϕ.An estimate of τ ϕ is necessary to conduct inference and Rfit implements the consistent estimator proposed by Koul et al. (1987).Denote this estimator by τϕ .Then Wald tests and confidence regions/intervals can be calculated.Let A Wald test of the general linear hypothesis where q = dim(M).and McKean (2011), the reduced model design matrix is easily obtained using a QR-decomposition on M T .We have implemented this methodology in Rfit.
Similar to the LS reduction in sums of squares, the rank-based test is based on a reduction of dispersion from the reduced to the full model.Let D( θϕ ) and D( βϕ ) denote the reduced and full model minimum dispersions, then the test is to reject

Computation
It can be shown that Jaeckel's dispersion function ( 3) is a convex function of β (Hettmansperger and McKean, 2011).Rfit uses optim with option `BFGS' to minimize the dispersion function.We investigated other minimization methods (e.g.Nelder-Mead or CG), however the quasi-Newton method works well in terms of speed and convergence.An orthonormal basis matrix, for the space spanned by the columns The R Journal Vol.4/2, December 2012 ISSN 2073-4859 of X, is first calculated using qr which leads to better performance in examples.The default initial fit is based on an L 1 fit using quantreg (Koenker, 2011).
Computations by Rfit of rank-based estimation and associated inference are illustrated in the examples of the next section.

Rfit computations of rank-based fitting of linear models
For the general linear model (1) the package Rfit obtains the rank-based estimates and inference, as described in the previous section.In this section we illustrate this computation for two examples.The first is for a simple linear model while the second is for a multiple regression model.
Example 1 (Telephone data).We begin with a simple linear regression example using the telephone data discussed in Rousseuw and Leroy (1987) These data represent the number of telephone calls (in tens of millions) placed in Belgium over the years .The data are plotted in Figure 1.There are several noticeable outliers which are due to a mistake in the recording units for the years 1964-1969.This is a simple dataset, containing only one explanatory variable, however it allows us to easily highlight the package and also demonstrate the robustness to outliers of the procedure.The main function of the package Rfit is rfit which, as the following code segment illustrates, uses syntax similar to lm. > abline(lm(calls ~year, data = telephone), + col = 2, lty = 2) > legend("topleft", legend = c("R", "LS"), + col = 1:2, lty = 1:2) q q q q q q q q q q q q q q q q q q q q q q q q 1950 1955 1960 Further, the output is similar to that of lm and can be interpreted in the same way.The estimate of slope is 0.146 (tens of millions of calls per year) with a standard error of 0.078.The t-statistic is the ratio of the two and the p-value is calculated using a tdistribution with n − 2 degrees of freedom.Hence one could conclude that year is a marginally significant predictor of the number of telephone calls.
The overlaid fitted regression lines in the scatter plot in Figure 1 demonstrate the robustness of the Wilcoxon fit and the lack of robustness of the least squares fit.
Example 2 (Free fatty acid data).This is a data set from Morrison (1983, p. 64) (c.f.Example 3.9.4 of Hettmansperger and McKean (2011)).The response variable is level of free fatty acid in a sample of prepubescent boys.The explanatory variables are age (in months), weight (in lbs), and skin fold thickness.For this discussion, we chose the Wilcoxon (default) scores for Rfit.Based on the residual and Q-Q plots below, however, the underlying error distribution appears to be right-skewed.In a later section we analyze this data set using more appropriate (bent) scores for a right-skewed distribution.
To begin with we demonstrate the reduction in dispersion test discussed in the previous section.As the code segment shows, the syntax is similar to that of the anova function used for reduced model testing in many of the parametric packages.
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.3 0.4 0.5 0.6 0.7 0.8 Figure 2: Studentized residuals versus fitted values for the free fatty acid data.
q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q −2 − Studentized residuals for rank-based fits are calculated in a way similar to the LS studentized residuals (see Chapter 3 of Hettmansperger and McKean, 2011).We have implemented these residuals in Rfit and demonstrate their use next.These are the residuals plotted in the residual and Q-Q plots in Figure 2 and Figure 3 respectively.The code is similar to that of least squares analysis.The func-tion fitted.valuesreturns the fitted values and residuals returns the residuals from the full model fit.The function rstudent calculates the studentized residuals.

> qqnorm(residuals(fitF))
As is shown in the plots, there are several outliers and perhaps the errors are from a right skewed distribution.We revist this example in a later section.

One-way ANOVA
Suppose we want to determine the effect that a single factor A has on a response of interest over a specified population.Assume that A consists of k levels or treatments.In a completely randomized design (CRD), n subjects are randomly selected from the reference population and n i of them are randomly assigned to level i, i = 1, . . .k.Let the jth response in the ith level be denoted by Y ij , j = 1, . .., i = 1, . . ., k.We assume that the responses are independent of one another and that the distributions among levels differ by at most shifts in location.
Under these assumptions, the full model can be written as where the e ij s are iid random variables with density f (x) and the parameter µ i is a convenient location parameter for the ith level, (for example, the mean or median of the ith level).Generally, the parameters of interest are the effects (contrasts), ∆ ii = µ i − µ i , i = i , 1, . . ., k. Upon fitting the model a residual analysis should be conducted to check these model assumptions.
Observational studies can also be modeled this way.Suppose k independent samples are drawn from k different populations.If we assume further that the distributions for the different populations differ by at most a shift in locations then Model (5) is appropriate.
The analysis for this design is usually a test of the hypothesis that all the effects are 0, followed by individual comparisons of levels.The hypothesis can be written as Robust fits based on scores more appropriate than the Wilcoxon for skewed errors are discussed later.Note that the results from a call to oneway.rfitinclude the results from the call to rfit.
> anovafit <-oneway.rfit(quail$ldl,quail$treat) Which may then be used for diagnostic procedures, such as the Q-Q plot of the studentized residuals in Figure 5.
> qqnorm(rstudent(anovafit$fit)) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −2 −1 0 1 2 With a p-value of 0.0164, generally, the null hypothesis would be rejected and the inference would pass to the comparisons of interest.Finally, we note that, the LS test of the null hypothesis has p-value 0.35; hence, with the LS analysis H 0 would not be rejected.In practice, one would not proceed with comparisons with such a large p-value.Thus, for this data set the robust and LS analyses have different interpretations.

Multiple comparisons
The second stage of an analysis of a one-way design usually consists of pairwise comparisons of the treatments.The robust (1 − α)100% confidence interval to compare the ith and i th treatments is given by Often there are many comparisons of interest.For example, in the case of all pairwise comparisons there The R Journal Vol.4/2, December 2012 ISSN 2073-4859 are ( k 2 ) confidence intervals.Hence, the overall fam- ily error rate is usually of concern.Multiple comparison procedures (MCP) try to control the overall error rate to some degree.There are many MCPs from which to choose; see Chapter 4 of Hettmansperger and McKean (2011) for a review of many of these procedures from a robust perspective.In Rfit we supply a summary function that adjusts confidence intervals and use three of the most popular such procedures: protected least significant difference (none); Tukey-Kramer (tukey); and the Bonferroni (bonferroni).These methods are described in many standard statistics texts.
Example 4 (LDL cholesterol of quail, continued).For the quail data, we selected the Tukey-Kramer procedure for all six pairwise comparisons.Use the code and example output is given below.The multiple comparison part of the output is: > summary(oneway.rfit(quail$ldl,quail$treat), + method = "tukey") The Tukey-Kramer procedure declares that the Drug Compounds I and II differ significantly.

Multi-way ANOVA
In this section, we consider a k-way crossed factorial experimental design.For these designs, the Rfit function raov computes the rank-based analysis for all 2 k − 1 hypotheses including the main effects and interactions of all orders.The design may be balanced or unbalanced.For simplicity, we briefly discuss the analysis in terms of a cell mean (median) model; see Hocking (1985) for details on the traditional LS analysis and Chapter 4 of Hettmansperger and McKean (2011) for the rank-based analysis.For this paper, we illustrate Rfit using a two-way crossed factorial design, but similarly Rfit computes the rank-based analysis of a k-way design.
Let A and B denote the two factors with levels a and b, respectively.Let Y ijk denote the response for the kth replication at levels i and j of factors A and B, respectively.Then the full model can be expressed as where e ijk are iid random variables with pdf f (t).Since the effects of interest are contrasts in the µ ij 's, these parameters can be either cell means or medians, (actually any location functional suffices).Rfit implements a reduction in dispersion tests for testing all main effects and interactions.
For the two-way model, the three hypotheses of immediate interest are the main effects hypotheses and the interaction hypothesis.We have chosen Type III hypotheses which are easy to interpret even for severely unbalanced designs.Following Hocking (1985), the hypothesis matrices M can easily be computed in terms of Kronecker products.As discussed in a previous section, for these tests the drop in dispersion test statistics can easily be constructed.We have implemented this formulation in Rfit.
Example 5 (Box-Cox data).Consider the data set discussed by Box and Cox (1964).The data are the results of a 3 × 4 two-way design, where forty-eight animals were exposed to three different poisons and four different treatments.The design is balanced with four replications per cell.The response was the log survival time of the animal.An interaction plot using the cell medians is presented in Figure 6.Obviously the profiles are not parallel and interaction is present.The output below displays the Wilcoxon ANOVA table, which indicates that interaction is highly significant, p = 0.0143, confirming the profile plot.On the other hand, the LS test F statistic for interaction is 1.87 with p = 0.1123.Hence, the LS test fails to detect interaction.

Writing score functions for Rfit
As discussed earlier, we must choose a score function for rank-based fitting.For most datasets the default option of Wilcoxon scores works quite well, however, occasionally choosing a different score function can lead to a more efficient analysis.In this section we first discuss score functions in general and then illustrate how the user may create his own score function.We have placed the score functions in an object of class "scores".A "scores" object consists of two objects of type function and an optional numeric object.The functions are the score function phi and it's derivative Dphi.The derivative is necessary in estimating τ ϕ .Below is what the class for Wilcoxon scores looks like.

> wscores
An object of class "scores" Slot "phi": function(u) sqrt( 12)*(u-0.5)Slot "Dphi": function(u) rep(sqrt(12),length(u)) Slot "param": NULL Other score functions included in Rfit are listed in Table 1.A plot of the bent score functions is provided in Figure 7.Other score functions can be plotted by getting the scores using the method getScores.For example the commands u<-seq(0.01,0.99,by=0.01)plot(u,getScores(nscores,u)) graphs the normal scores.Next we illustrate how to create the score function for the bent scores.Bent scores are recommended when the errors come from a skewed distribution.

Score
An appropriate bent score function for skewed distribution with a right heavy tail is φ(u) = 4u − 1.5 if u ≤ 0.5 0.5 if u > 0.5 The following code segment defines the scores.

Figure 1 :
Figure 1: Scatter plot of the telephone data with overlaid regression lines.

Figure 3 :
Figure 3: Normal Q-Q plot of the studentized residuals for the free fatty acid data.

Figure 4 :
Figure 4: Comparison boxplots for quail data.Using Wilcoxon scores, we fit the full model.The summary of the test of hypotheses (6) as computed by the Rfit function oneway.rfitfollows.The resulting Q-Q plot (Figure 5) of the studentized residuals indicates that the random errors e ij have a distribution.> data(quail) > oneway.rfit(quail$ldl,quail$treat) Call: oneway.rfit(y= quail$ldl, g = quail$treat)

Figure 5 :
Figure 5: Normal Q-Q plot of the studentized residuals for the quail data.

Figure 6 :
Figure 6: Interaction Plot for Box-Cox Data.

Figure 7 :
Figure 7: Plots of four bent score functions.

Table 1 :
Table of available score functions.Unless otherwise noted, distribution is assumed to be symmetric.