R/wls.R

In stevemisc: Steve's Miscellaneous Functions

#' Get Weighted Least Squares of Your OLS Model
#'
#' @description \code{wls()} takes an OLS model and re-estimates it using a weighted least squares
#' approach. Weighted least squares is often a "textbook" approach to dealing with the presence of
#' heteroskedastic standard errors, for which the weighted least squares estimates are compared
#' to the OLS estimates of uncertainty to check for consistency or potential inferential implications.
#'
#' @details The function *should* be robust to potential model specification oddities (e.g. polynomials
#' and fixed effects). It also should perform nicely in the presence of missing data, if and only
#' if `na.action = na.exclude` is supplied first to the offending OLS model supplied to the function
#' for a weighted least squares re-estimation.
#'
#' @return \code{wls()} returns a new model object that is a weighted least squares re-estimation
#' of the OLS model supplied to it.
#'
#' @author Steven V. Miller
#'
#' @param mod a fitted OLS model
#'
#' @examples
#'
#' M1 <- lm(mpg ~ ., data=mtcars)
#' M2 <- wls(M1)
#'
#' summary(M2)

wls <- function(mod) {
  if (!identical(class(mod), "lm")) {
    stop("Weighted Least Squares only makes sense in the context of OLS, and this model/object does not appear to be OLS. As of right now, this function works for just those OLS models generated by the lm() function in base R.")
  }

  resid <- resid(mod)
  fitted <- fitted(mod)

  WM <- lm(abs(resid) ~ fitted, na.action=na.exclude)

  wts <- 1/(fitted(WM)^2)

  A <- eval(mod$call$data)
  A$wts <- 1/(fitted(WM)^2)

  WLSM <- update(mod, ~., data=A, weights = wts)
  return(WLSM)


}