--- title: "Feasible GLS under multiplicative heteroskedasticity" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Feasible GLS under multiplicative heteroskedasticity} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5 ) options(hcinfer.use_emoji = FALSE) ``` Feasible generalized least squares (FGLS) complements the heteroskedasticity-consistent (HC) estimators in `hcinfer`. The HC estimators keep the ordinary least squares (OLS) coefficients and only robustify their covariance, assuming nothing about the form of the variance function. Feasible GLS instead models the conditional variance as an exponential function of observed regressors, re-weights the regression accordingly, and can be more efficient than OLS when that variance model is adequate. The `gls_mult()` function offers two estimators of this model, Harvey's two-step procedure and Gaussian maximum likelihood, and the maximum likelihood fit supplies a proper `logLik()` together with `AIC()` and `BIC()`. The method follows Harvey (1976) and Cribari-Neto and Pereira (2019). ## HC inference versus feasible GLS The two approaches answer different questions and suit different situations. Calling `hcinfer()` or `vcov_hc()` returns the OLS point estimates together with a heteroskedasticity-consistent sandwich covariance that is valid under heteroskedasticity of unknown form, so it makes no commitment to how the variance depends on the regressors. Calling `gls_mult()` instead estimates an explicit variance model, returns re-weighted coefficients that can be more efficient than OLS, and reports a model-based covariance that is valid when the variance model is correctly specified. As a practical rule, the HC estimators are the natural choice when you do not wish to model the variance and want inference that is agnostic about its form, whereas feasible GLS is attractive when a multiplicative variance model is plausible and the efficiency gain matters. ## The multiplicative variance model Feasible GLS in `hcinfer` is built on a linear mean model paired with a multiplicative model for the conditional variance. The mean model is $$ y = X \beta + e, \qquad E(e_t) = 0, \qquad \operatorname{Var}(e_t) = \sigma_t^2, $$ where $y$ is the $n \times 1$ response, $X$ is the $n \times p$ mean design matrix, and $\beta$ is the vector of mean coefficients. The variance is modelled as an exponential function of a second set of regressors, $$ \sigma_t^2 = \exp(\eta_t), \qquad \eta_t = z_t^\top \gamma, $$ so that the log-variance $\eta_t$ is linear in the dispersion regressors $z_t$, the rows of a dispersion design matrix $Z$, with coefficients $\gamma$. Because the response has response units, $\exp(\eta_t)$ has squared response units and $\exp(\eta_t/2)$ is the conditional standard deviation in response units. The exponential link keeps the fitted variances positive for any value of $\gamma$. The `variance` argument selects $Z$: a one-sided formula such as `~ income_scaled` builds $Z$ from the named terms, while the default `variance = NULL` sets $Z = X$ and reuses the mean design. This is a model for the variance, not an assumption that the residuals or leverage complements follow any particular distribution. ## Two estimators Both estimators target the same mean and variance model but reach it by different routes. The two-step estimator is a single pass through an auxiliary regression, whereas maximum likelihood locally optimises the Gaussian profile likelihood and is accepted only at a stationary solution. Maximum likelihood is the default, and it is also the fit that supports the information criteria. ### Harvey two-step The two-step estimator of Harvey (1976) begins from the OLS residuals $\hat e_t$ and regresses $\log(\hat e_t^2)$ on the dispersion regressors $Z$. Because the expectation of $\log \varepsilon^2$ under normality is not zero, the intercept of this auxiliary regression is bias-corrected by adding $c = -\operatorname{digamma}(1/2) - \log 2 \approx 1.2704$, the bias of $\log \varepsilon^2$ under normality. The fitted log-variances yield the weights $\hat w_t = \exp(-z_t^\top \hat\gamma)$, and a single weighted least squares pass with these weights delivers the mean coefficients, with reported dispersion covariance $\tfrac{\pi^2}{2}(Z^\top Z)^{-1}$. ### Maximum likelihood Maximum likelihood optimises the Gaussian log-likelihood of the same model. The mean coefficients $\beta$ are profiled out by weighted least squares at each candidate $\gamma$, so a BFGS optimiser searches only over the dispersion parameters $\gamma$, and the reported asymptotic dispersion covariance is $2(Z^\top Z)^{-1}$. Under the model and standard regularity conditions, maximum likelihood is asymptotically efficient. For an accepted locally optimised stationary fit, `logLik()` returns the likelihood evaluated at that solution, which enables `AIC()` and `BIC()` without claiming that the stationarity check has certified a global maximum. ## Data and model The examples use the `PublicSchools2` data, which have complete observations for the 50 U.S. states and the District of Columbia and include a `south` indicator that the original `PublicSchools` data lack. The `expenditure` response is annual expenditure per student enrolled in K-12 public schools for 2025, measured in U.S. dollars. Income is divided by ten thousand so that the estimated coefficients are on a readable scale. The mean model regresses expenditure on rescaled income and the regional indicator, and printing the fitted `lm` object shows the OLS starting point that `gls_mult()` will refine. ```{r} library(hcinfer) schools <- PublicSchools2 schools$income_scaled <- schools$income / 10000 fit <- lm(expenditure ~ income_scaled + south, data = schools) fit ``` ## Maximum likelihood fit With no `estimator` argument, `gls_mult()` fits the model by maximum likelihood, its default estimator, and optimizes the profile likelihood with the algorithm chosen by `method` (BFGS by default). The printed object gives a compact view of the mean and dispersion coefficients together with the convergence diagnostics. ```{r} gls_fit <- gls_mult(fit) gls_fit ``` ```{r} summary(gls_fit) ``` The summary reports fitted conditional variances in squared U.S. dollars, or USD squared, because `expenditure` is measured in USD and $\exp(\eta_t)$ is a variance. On the standard-deviation scale, the corresponding values range from about USD `r formatC(min(sqrt(gls_fit$fitted_variances)), digits = 0, format = "f", big.mark = ",")` to USD `r formatC(max(sqrt(gls_fit$fitted_variances)), digits = 0, format = "f", big.mark = ",")` and are comparable to the OLS residual standard deviation of about USD `r formatC(stats::sigma(fit), digits = 0, format = "f", big.mark = ",")`. The million-scale variance entries are therefore expected for this response scale and are not evidence of an exponentiation error. The same report gives the mean and dispersion blocks, normal Wald tests, and the likelihood evaluated at the accepted local solution together with AIC and BIC. ## Extracting the fit The usual extractor functions work on the fitted object and, where relevant, take a `model` argument to choose between the mean and dispersion blocks. ```{r} coef(gls_fit) # mean coefficients coef(gls_fit, model = "dispersion") # log-variance coefficients ``` ```{r} vcov(gls_fit) # mean covariance vcov(gls_fit, model = "dispersion") # dispersion covariance ``` ```{r} tests(gls_fit) confint(gls_fit) ``` ```{r} head(fitted(gls_fit)) head(gls_fit$fitted_variances) # estimated conditional variances head(gls_fit$weights) # GLS weights exp(-eta) ``` By default `coef()` and `vcov()` return the mean block, while `model = "dispersion"` selects the log-variance coefficients and their covariance. The `tests()` and `confint()` methods report normal Wald inference on the mean coefficients, using the model-based standard errors. The fitted object also stores the conditional variances in `fitted_variances` and the GLS weights $\exp(-\eta_t)$ in `weights`, so the estimated variance structure is available for diagnostics or plotting. These weights have inverse-squared-response units, and multiplying all of them by one positive constant leaves the weighted least squares coefficients unchanged, so only their relative scale affects those coefficients. ## Harvey two-step Setting `estimator = "two_step"` fits the same model with Harvey's non-iterative estimator. ```{r} two_step <- gls_mult(fit, estimator = "two_step") coef(two_step) ``` ```{r} data.frame( term = names(coef(gls_fit)), ml = coef(gls_fit), two_step = coef(two_step) ) ``` The two-step estimator uses one auxiliary regression, whereas maximum likelihood locally optimises the profile likelihood and applies the stationarity guard. In this example, the income slopes are similar, but the intercept and `south` coefficients differ more noticeably, so the table demonstrates estimator dependence rather than equivalence. Information criteria are defined only for the accepted maximum likelihood fit because the two-step point does not optimise the likelihood and therefore carries no comparable `logLik()`. ## Custom variance models The `variance` argument decouples the dispersion model from the mean model, so the analyst can let the variance depend on a different, and usually smaller, set of regressors. ```{r} gls_income <- gls_mult(fit, variance = ~ income_scaled) coef(gls_income, model = "dispersion") ``` Here the variance is modelled as a function of income alone, while the mean model still includes both income and the regional indicator. The default `variance = NULL` instead reuses the full mean design as $Z$, which is convenient but not always the most parsimonious choice for the variance. ## Model selection with AIC and BIC Because each maximum likelihood fit carries a proper `logLik()` with $p + q$ degrees of freedom, where $p$ counts the mean coefficients and $q$ the dispersion coefficients, the base `AIC()` and `BIC()` generics compare competing variance specifications directly. The three fits below share the same mean model but differ in their variance model: the full dispersion model on income and region, a reduced model on income alone, and the special case `variance = ~ 1`, which forces a constant variance and so reduces the fit to the homoskedastic Gaussian model. ```{r} full <- gls_mult(fit) # variance ~ income + south income_only <- gls_mult(fit, variance = ~ income_scaled) homoskedastic <- gls_mult(fit, variance = ~ 1) AIC(full, income_only, homoskedastic) BIC(full, income_only, homoskedastic) ``` The specification with the smallest criterion is preferred, so these tables let the data adjudicate between richer and sparser variance models on the same footing. Because they rest on the likelihood, the criteria are defined only for maximum likelihood fits. Meaningful comparisons also require the competing fits to have reached comparable likelihood solutions. ## Comparison with HC standard errors It is instructive to place the feasible GLS standard errors next to the OLS and HC standard errors for the same model. ```{r} data.frame( term = names(coef(fit)), ols = sqrt(diag(vcov(fit))), fgls_ml = sqrt(diag(vcov(gls_fit))), hcbeta = sqrt(diag(vcov(hcinfer(fit, type = "hcbeta")))), hc3 = sqrt(diag(vcov(hcinfer(fit, type = "hc3")))) ) ``` The feasible GLS standard errors come from a re-weighted fit under an explicit variance model, whereas the HC standard errors, here HC$\beta$ and HC3, keep the OLS coefficients and only robustify their covariance without modelling the variance. Agreement between the reported standard errors is descriptive and does not by itself establish that the variance model is correctly specified. Systematic differences can motivate inspection of the variance specification and influential high-leverage observations, but they do not identify either explanation on their own. ## Practical guidance Several practical points help in routine use. Maximum likelihood is the recommended default because it is efficient under the model and provides information criteria, while the two-step estimator is a useful non-iterative alternative. The mean covariance that `gls_mult()` reports is a model-based plug-in, valid when the variance model is correctly specified rather than an HC sandwich, so inference that should be agnostic about the variance form still belongs to `hcinfer()`. Scaling the dispersion regressors sensibly can improve numerical behaviour, because extreme scaling may require optimiser tuning through `control`, such as a tighter `reltol`. The convergence code and scale-invariant score check establish only a locally optimised stationary fit; they do not prove that the accepted solution is a global maximum. For the underlying theory of the HC estimators see `vignette("hcinfer-methodology", package = "hcinfer")`, and for the full argument list see `?gls_mult`.