Weighted Least Squares wrapper¶

Source: R/fitWLS.R · Help: ?fitWLS

Refit a model by feasible generalised least squares, weighting each observation by the inverse of an estimated error variance.

Usage¶

fitWLS(
  model,
  data = NULL,
  var_formula = NULL,
  form = c("exponential", "power")
)

Arguments¶

model: A fitted model of class lm, or a glm with the Gaussian family and the identity link. Since 0.17.0 it refuses a MASS::rlm() fit, an M-estimate, whose residuals are not the least-squares residuals it fits its variance function to, and an mgcv::gam() fit with penalised smooth terms, which it refitted by weighted least squares on the basis of its smooths without the penalty; a gam that is the least-squares fit of its design, as one whose smooth terms are all fixed (s(x, fx = TRUE)), is refitted on that design (decisions 4, 37 and 38 of the design of the model context).
data: Optional base::data.frame holding the variables named in var_formula. It is needed only when they are not in the model frame. Since 0.16.0 its rows are matched to the observations of the fit as the tests that read data match them, data with row names of their own in any order and data that number their rows (a tibble) in the order of the fit, and other data are refused, where numbered data in another order were read by position. A variable of var_formula that neither data nor the model frame holds is refused, where it was found in the workspace. Without var_formula, data is not read, but it is matched in the same way, and other data are refused, where they were ignored.
var_formula: One-sided formula naming the variance regressors, for example ~ x. NULL, the default, uses the regressors of model.
form: Character scalar choosing the variance function: "exponential" (the default) or "power". See the section on variance functions.

Value¶

A new lm object fitted with weights. The estimated variances are attached as the "variance_model" attribute, one per row, as weights() of the result gives the weights: on a fit with na.action = na.exclude, NA at the rows the fit dropped (since 0.17.0; before, one per observation of the fit). The fitted variance function (its form, regressors and coefficients) is attached as the "variance_function" attribute, which is what performVarianceFormTest() reads, and what the bootstrap tests read to estimate the weights again in each replicate. Its regressors, the variance regressors on the scale of the form, and index, the fitted log-variance \(\hat\gamma_0 + z_i^\top \hat\gamma\), have one row, or value, per observation of the fit, also on a fit with na.action = na.exclude, whose rows they are not padded to. The estimated variances are proportional to the exponential of index, unless the variance function could not be used and the fit fell back to equal weights, which the record says with usable = FALSE and its reason; index is then the fitted log-variance that was refused, or NA. When model was fitted with weights, the attribute also holds them, as initial_weights, one per observation of the fit too.

Details¶

The variance is modelled, not read off the residuals directly. A single squared residual is a one-degree-of-freedom estimate of \(\sigma_i^2\) and far too noisy to invert: weighting by \(1/e_i^2\) hands almost all of the weight to whichever observations the initial fit happened to reproduce most closely. This function instead fits a variance function to the squared residuals and takes \(\hat\sigma_i^2\) from its fitted values, the standard feasible-GLS recipe; weights are \(1/\hat\sigma_i^2\).

Variance functions¶

form chooses the shape of the variance function and var_formula the variables \(z_i\) it depends on. Both forms are linear in the logarithm of the variance and carry their own constant.

"exponential": \(\log \sigma_i^2 = \gamma_0 + z_i^\top \gamma\), the multiplicative model of Harvey (1976), estimated by regressing \(\log e_i^2\) on \(z_i\). The log scale keeps the fitted variances positive without constraining the auxiliary regression. This is the default, and with var_formula = NULL it is what fitWLS() has computed since 0.9.0.
"power": \(\log \sigma_i^2 = \gamma_0 + \sum_j \delta_j \log z_{ij}\), that is \(\sigma_i^2 \propto \prod_j z_{ij}^{\delta_j}\), the model of Park (1966). It is the exponential form in the logarithms of the regressors, which must be strictly positive. \(\delta = 2\) is the textbook case of a standard deviation proportional to \(z\).

Other log-linear shapes are written through var_formula, for example ~ x + I(x^2) or ~ sqrt(x). The additive model \(\sigma_i^2 = \alpha_0 + z_i^\top \alpha\) is not offered: nothing keeps its fitted variances positive, and in the validation study it was not usable in a material share of samples drawn from an additive model itself (inst/validation/README.md).

Residuals that are numerically zero are floored before the logarithm is taken, with a warning.

Weights estimated this way are consistent under a correctly specified variance function, so standard errors from the returned fit are usable. They were not before 0.9.0, when the weights were the raw inverse squared residuals: the weighted residual sum of squares then collapsed towards \(n\) regardless of the data, and nominal 95% intervals covered the truth about 10% of the time.

Whether the variance function is correctly specified can be tested: pass the result to performVarianceFormTest().

A test of constant variance on the result takes its weights as known, although they were estimated from the same data, so its reference distribution is only approximate: see the validation study in inst/validation/README.md. The bootstrap tests performWildBootstrapTest(), performWhiteTestBootstrap() and bootstrapTestStatistic() allow for the estimation. Since 0.15.0 they run fitWLS() again on each replicate by default: they fit it as model was fitted, fit the recorded variance function again to its residuals, on the recorded regressors, and refit it with the weights that gives. Their argument reestimate_weights = FALSE holds the weights fixed instead, as they did before 0.15.0.

If the variance function cannot be fitted, or yields no usable variation, the function falls back to equal weights, which reduces the result to the original OLS fit.

The weighted fit keeps the design and the offset of model, whether the offset was given as an argument, lm(..., offset = o), or written in the formula, offset(o). Its coefficients, fitted values and residuals are those lm() returns for the same call with the estimated weights added. On a model fitted with weights, the refit uses the estimated weights in their place, and weights() of the result returns them. Since 0.16.0 the response, the design and the offset are read from the fit, as the tests read them, and a fit with model = FALSE whose model frame, rebuilt from the data it was fitted to, no longer holds its response, because those data changed after the fit, is refused, where it was refitted to the response of the changed data; a fit with model = FALSE whose data are gone is refitted from data that give its response and fitted values.

A glm with the Gaussian family and the identity link is a least-squares fit, with the design and residuals of the lm, and since 0.16.0 the result for it is the lm returned for that lm, prior weights included, with the call of that lm: stats::lm(), without the arguments only glm() takes, so that update() fits an lm again. Before, it kept the class and the call of the glm without its family, which summary() could not read and every test refused. Any other glm is refused.

References¶

Harvey, A. C. (1976). Estimating regression models with multiplicative heteroscedasticity. Econometrica, 44(3), 461–465. doi:10.2307/1913974

Park, R. E. (1966). Estimation with heteroscedastic error terms. Econometrica, 34(4), 888. doi:10.2307/1910108

Examples¶

data(mtcars)
m <- lm(mpg ~ wt + qsec, data = mtcars)
wls <- fitWLS(m)
summary(wls)

# Variance proportional to a power of one regressor
wls_power <- fitWLS(m, var_formula = ~ wt, form = "power")
attr(wls_power, "variance_function")$coefficients