Example 60.2 Iteratively Reweighted Least Squares

With the NLIN procedure you can perform weighted nonlinear least-squares regression in situations where the weights are functions of the parameters. To minimize a weighted sum of squares, you assign an expression to the _WEIGHT_ variable in your PROC NLIN statements. When the _WEIGHT_ variable depends on the model parameters, the estimation technique is known as iteratively reweighted least squares (IRLS). In this situation you should employ the NOHALVE option in the PROC NLIN statement. Because the weights change from iteration to iteration, it is not reasonable to expect the weighted residual sum of squares to decrease between iterations. The NOHALVE option removes that restriction.

Examples where IRLS estimation is used include robust regression via M-estimation (Huber 1964, 1973), generalized linear models (McCullagh and Nelder 1989), and semivariogram fitting in spatial statistics (Schabenberger and Pierce 2002, Sect. 9.2). There are dedicated SAS/STAT procedures for robust regression (the ROBUSTREG procedure) and generalized linear models (the GENMOD and GLIMMIX procedures). Examples of weighted least-squares fitting of a semivariogram function can be found in Chapter 95, The VARIOGRAM Procedure.

In this example we show an application of PROC NLIN for M-estimation only to illustrate the connection between robust regression and weighted least squares. The ROBUSTREG procedure is the appropriate tool to fit these models with SAS/STAT software.

M-estimation was introduced by Huber (1964, 1973) to estimate location parameters robustly. Beaton and Tukey (1974) applied the idea of M-estimation in regression models and introduced the biweight (or bisquare) weight function. See Holland and Welsch (1977) for this and other robust methods. Consider a linear regression model of the form

In weighted least-squares estimation you seek the parameters that minimize

where is the weight associated with the th observation. The normal equations of this minimization problem can be written as

In M-estimation the corresponding equations take on the form

where is a weighing function. The Beaton-Tukey biweight, for example, can be written as

Substitution into the estimating equation for M-estimation yields weighted least-squares equations

The biweight function involves two constants, and . The scale can be fixed or estimated from the fit in the previous iteration. If is estimated, a robust estimator of scale is typically used. In this example is fixed at . A common value for the constant is .

The following DATA step creates a SAS data set of the population of the United States (in millions), recorded at 10-year intervals starting in 1790 and ending in 1990. The aim is to fit a quadratic linear model to the population over time.

   title 'U.S. Population Growth';
   data uspop;
      input pop :6.3 @@;
      retain year 1780;
      year   = year+10;
      yearsq = year*year;
      datalines;
   3929 5308 7239 9638 12866 17069 23191 31443 39818 50155
   62947 75994 91972 105710 122775 131669 151325 179323 203211
   226542 248710
   ;

The PROC NLIN code that follows fits this linear model by M-estimation and IRLS. The weight function is set to a zero or nonzero value depending on the value of the scaled residual. The NOHALVE option removes the requirement that the (weighted) residual sum of squares must decrease between iterations.

   title 'Beaton/Tukey Biweight Robust Regression using IRLS';
   proc nlin data=uspop nohalve;
      parms b0=20450.43 b1=-22.7806 b2=.0063456;
      model pop=b0+b1*year+b2*year*year;
      resid = pop-model.pop;
      sigma = 2;
      k     = 4.685;
      if abs(resid/sigma)<=k then _weight_=(1-(resid / (sigma*k))**2)**2;
      else _weight_=0;
      output out=c r=rbi;
   run;

Parameter estimates from this fit are shown in Output 60.2.1, and the computed weights at the final iteration are displayed in Output 60.2.2. The observations for 1940 and 1950 are highly discounted because of their large residuals.