Robust Fitting Iteratively Reweighted Least Squares

Question

I'm interested in Mathematica's capabilities with regard to robust fitting vs. outliers. In the Statistical Model Analysis Tutorial under the Generalized Linear Models section heading, there is this sentence in a paragraph talking about Options for GeneralizedLinearModelFit:

Parameter estimates are obtained via iteratively reweighted least squares with weights obtained from the variance function of the assumed distribution.

This seems to say that GeneralizedLinearModelFit does iteratively reweighted least squares fitting. However, I can't see any options for iteratively reweighting data points, nor any indication that the weight applied to the points is anything other than 1.

For example:

$Version
(*11.2.0 for Linux x86 (64-bit) (September 2, 2017)*)

data = RandomVariate[NormalDistribution[100, 30], 200];
modData=ReplacePart[data,1->10^6];

Through@{Mean,StandardDeviation}@data
(*{99.529,30.395}*)

Through@{Mean,StandardDeviation}@modData
(*{5098.99,70703.6}*)

Through@{Median,SnDispersion}@modData
(*{99.2765,31.6749}*)

Through@{LinearModelFit,GeneralizedLinearModelFit}[data,{},x]//Normal
(*{99.529,99.529}*)

Through@{LinearModelFit,GeneralizedLinearModelFit}[modData,{},x]//Normal
(*{5098.99,5098.99}*)

GeneralizedLinearModelFit[modData,{},x,DispersionEstimatorFunction->Function[{y,yhat,w},SnDispersion[y-yhat]]]//Normal
(*5098.99*)

GeneralizedLinearModelFit[modData,{},x,DispersionEstimatorFunction->Function[{y,yhat,w},0]]//Normal
(*5098.99*)

Clearly, the general method isn't reweighting according to the variance. The variance is gigantic in the modified data set and no reweighting occurred.

Secondly, the dispersion estimator function appears to do nothing.

Lastly, I thought about modifying the Weights function, but it only takes one original data point as its argument in functional form:

Weights->func associates weight func[xi1,xi2,…,yi] with the i^(th) data element.

The weight function, as documented, does not have access to the current iteration's prediction $\hat{y}_i$ or at least the current best fit parameters/model function, which is something that would be needed for reweighting (see section 2.2 Table 1). Certainly this information could be provided via a side effect in the DownValues of a named function or a call to a cleverly constructed side function that happens to have the current state of the model, but given the construction of the options for Weight, the weight function is unlikely to be called each iteration of GeneralizedLinearModelFit anyway.

So, what am I misunderstanding? How do we get GeneralizedLinearModelFit to do iteratively reweighted least squares regression to ignore outliers?

Should we just repeatedly call the function and manually update the weights?

Answer 1

I figured out this can be done using the NormFunction option of FindFit. Using the Huber weighting function of the last reference in the question, the solution would be

modData=ReplacePart[RandomVariate[NormalDistribution[100, 30], 200],1->10^6];
huberNorm[arg:{__?NumericQ}]:=With[{abs=Abs@arg},With[{k=2 Median@abs},Norm[Min[#,k]&/@abs]]]
FindFit[modData,t,{t},{x},NormFunction->huberNorm]
(*FindFit lmnl The model t is linear in the parameters t but a nonlinear method or non-Euclidean norm was specified, so nonlinear methods will be used.*)
(*FindFit lstol The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the norm of the residual. You may need more than MachinePrecision digits of working precision to meet these tolerances.*)
(*{t->100.007}*)