What is the origin of squaring centred data as way to model variances instead of means? I recently came across this Answer by @mpiktas wherein he suggested a transformation of $y_i \rightarrow y_i^{\prime}$
$$y_i^{\prime} = (y_i - \overline{y})^2$$
followed by fitting a model for change in mean of $y_i^{\prime}$ as a way to model a change in variance.
The same idea is also present in some versions of Levene's test (although the more robust versions use medians and absolute deviations).
I was wondering then, is there an authoritative or canonical reference for this approach to modelling variance in regression settings? If it matters, the specific application I have in mind is for analysis of irregular time series, but I plan to do the analysis using GLS or GAMs.
 A: The average of squared deviations or residuals provide a biased estimate of the variance when data are IID. The simple degrees of freedom correction removes the bias. Checking the IID assumption means using what you know about the structure of the residuals to check if the average of the residual varies between groups.
In time-series data or geospatial data, you can use variograms. In clustered data, you can calculate intraclass correlation. For sandwich based standard errors, the empirical variance derivation uses a diagonal matrix of residuals as a plug in estimate for variance under no assumptions about IID. In mixed effects modeling, the EM algorithm is used to iteratively estimate fixed effects, and use those residuals to estimate the covariance components of the residuals under some parametric assumptions.
There's nothing more to the rationale of these method than the fact that these are the values used to calculate confidence intervals and p-values under certain assumptions which are verifiable when the data are correlated.
"Robustness" is a funny concept. I think the Huber-White sandwich based standard errors are the only error estimates that are worthy of being called robust. In my simulations, it's usually contrived and data-driven approaches used to justify modeling absolute deviations, or other estimates of "scale". I find them difficult to interpret and, most importantly, less efficient at estimating what I'm interested in.
