Abnormally Distributed

Least Absolute Squares

Least absolute squares, also called median regression, is one way to conduct a regression analysis on data that include outliers or violate the assumption of homogeneity of variance. Median regression is just an instance of quantile regression where the quantile of interest is the 50th quantile. Please see the page Quantile regression for more information about this approach. There are, however, many more approaches. The major ones will be reviewed here. It is worth noting that least absolute squares and quantile regression more generally are L1 norm based methods. The L1 norm is any method based on minimizing the absolute value of the difference between the observed \(y\) values and predicted \(y_{pred}\) values. By contrast, the L2 norm

Weighted Least Squares

As mentioned before, heteroskedasticity increases the variance of the residuals in an uneven way. This makes the standard errors used in significance tests or confidence interval construction invalid. They may lead to the correct decision, but they may also cause a Type I or Type II error (in frequentist parlance).

One solution is to adopt a weighted least squares approach with weights to balance the influence of observations of \(x\) with larger magnitude to compensate for the variance of \(y\) increasing with \(x\). A common formula for the weights in this situation is for the \(i_{\text{th}}\) weight to be the reciprocal of the square root of the \(i_{\textit{th}}\) observation of predictor \(x\). The method of calculation is shown below.

\[ \text{VCOV}_{\text{wls}} = \underbrace{ \begin{pmatrix} \sigma^2_1&0&\cdots&0\\ 0&\sigma^2_2&\cdots&\vdots\\ \vdots&\cdots&\ddots&\vdots\\ 0&\cdots&\cdots&\sigma^2_i\\ \end{pmatrix}}_{\text{VCOV}_{\text{ols}}} ~~~~~~ \cdot ~~~~~~ \underbrace{ \begin{pmatrix} w_1&0&\cdots&0\\ 0&w_2&\cdots&\vdots\\ \vdots&\cdots&\ddots&\vdots\\ 0&\cdots&\cdots&w_i\\ \end{pmatrix}}_{\textit{w}_i ~ = ~ \frac{1}{\sqrt{x_i}}} \]

Below is a comparison of two model fits, one the ordinary least squares fit and the other the weighted least squares fit. The same simulated personality data from the quantile regression page will be used.

data = read.csv("https://tinyurl.com/ybkvf6e2", header=TRUE)

wls = lm(Neuroticism ~ Age,data, weights=1/sqrt(Age))

ols = lm(Neuroticism ~ Age,data)

Another solution is to use heteroskedastic consistent variance-covariance matrices to calculate standard error. This is done through a family of estimators called “sandwich” estimators. Although this won’t change the coefficient estimates it will result in standard errors (and subsequently p-values or confidence intervals) corrected for the heteroskedasticity. For the sake of space I will show how to obtain this correction, but won’t go into the math behind it.

vcov.rob = vcovHC(ols, method="HC4m")
sammich = coeftest(ols, vcov.rob)
kable(tidy(sammich))%>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Least Trimmed Squares

One of the simplest approaches to robust regression is least trimmed squares (LTS). The error term is the squared value of the residuals (see equation below), which like ordinary least squares is an L2 norm based method because it is based on the squared residuals, rather than absolute values.

\[ \sum_{i(k)=1}^n (y − y_{pred})^2_{i(k)} \]

The error term in LTS is minimized by finding the subset \(k\) of the data set that produces the smallest value of the error term. Alternatively the amount of trimming can be specified manually, but this may not give optimal results. However, LTS has relatively low efficiency, meaning it takes more observations to achieve stable results. This is somewhat intuitive, since trimming observations does reduce the original sample size. LTS does have the advantage of a very high breakdown point of 0.5. This means that LTS is insensitive to influence by outliers, provided that the outliers constitute less than 50 % of the data (Mount et al. 2014). See Venables and Ripley (2002) and Wilcox (2012b) for more information on least trimmed squares regression.

For an example the diabetes data from Bradley Efron’s Computer Age Statistical Inference website will be used, with some “corruptions” added to the data for an example predicting ldl cholesterol from hdl. All variables are scaled (as z-scores).

diabetes = read.csv("https://tinyurl.com/ycl5qcz5", header=TRUE)

An lts model can be fit with the ltsReg function in the robustbase package. Below is a comparison of the leasts squares and least trimmed squares fits. An LTS function is also available in the MASS package, but requires bootstrapping to obtain p-values.

lts = ltsReg(formula = ldl ~ hdl, data = diabetes)

ols = tidy(lm(formula = ldl ~ hdl, data = diabetes))

### M-Estimators and Extensions

Weighted least square can also be applied to problems with regression outliers and leverage points. To do this, a set of optimal weights must be found. This is done through an iteratively weighted least squares (IWLS) approach where a family of maximum-likelihood like estimators (M-Estimators) are used to find an optimal set of weights.

Generally speaking, M-Estimators work by minimizing a particular function \(\xi\) on the residuals (which are scaled by a scale parameter \(\tau\)) (Wilcox 2012b).

\[ \sum_{i=1}^n \xi\Big(\frac{y-y_{{pred}}}{\tau}\Big)_i^2 \]

The most popular class of \(\xi\) functions for M-Estimation is the \(\psi\) family of estimators. The two most common examples are Tukey’s bisquare and Huber’s psi. The particular mathematical details will not be reviewed here, but both Tukey’s and Huber’s measures are robust measures of central tendency. For more details Wilcox (2012a), Wilcox (2012b), and Venables and Ripley (2002) give thorough coverage of the particulars.

The purpose of \(\psi\) is to find a set of weights which reduce the influence of regression outliers on the final estimates. This is done as an iterative process for different values for \(w_i\).

M-Estimators are useful for finding an optimal solution when there are a few regression outliers. They can fail more easily in the presence of bad leverage points. S-Estimators, or scale estimators, attempt to find values which minimize the scale of the residuals. This means the predictor variables must be taken into account to find the coefficients that minimize the residual variance. There are a few approaches to this but generally S-Estimators depend on IWLS. However, S-Estimators can be inefficient compared to M-Estimators.

The solution that has emerged is the Modified M-Estimators, or MM-Estimators. MM-Estimators use S-Estimators to find coefficient values which minimize the residual standard deviation of M-Estimators. The process begins by calculating a candidate set of M-Estimator solutions. The M-Estimators need a measure of scale \(\tau\), so S-Estimators find an optimal solution to the candidate set. This increases the efficiency of S-Estimation because it does not have to search a parameter space that grows with sample size. Once the optimal scale is chosen, M-Estimation finds the most likely values for the coefficients. It is capable of achieving a breakdown point of .5, and confidence intervals and p-values can be obtained with much less effort and time than bootstrapping the LTS method.

As an example the same data used for the least trimmed squares example will be used. The function rlm can be found in the MASS package.

rlm = rlm(ldl ~ hdl, data = diabetes, psi="psi.huber", scale.est="Huber")

Ridge Regression

Ridge regression is a solution to the problem of collinearity. It adds a penalty term (see equation below) to the sum of squared errors formula for each predictor. This regularizes the estimate so that there is greater tolerance to the effects of collinearity.

\[ \underbrace{\sum_{i=1}^n (y_i - \sum_{j=1}x_{ij}\beta_j)^2}_{\text{L2 norm SSE}} + \underbrace{\lambda \sum_{j=1} \beta_j^2}_\text{L2 Penalty} \]

The penalty is known as an L2 penalty because it is based on the sum of squares, rather than an absolute value based term penalty (L1 penalty). The penalty adds the product of lambda (a tuning constant used to select how strong the penalty should be) and the squared coefficient for the \(\textit{j}_{\textit{th}}\) predictor. This has the effect of adding shrinkage to each coefficient that pulls its value towards zero.

An Example

Consider a data set based on the iris data set in R. The data can be downloaded at the link below directly into R. Ridge regression is evaluated based on its out of sample predictive error. Although ridge regression under-fits the model on the data it is given, this actually results in better out of sample predictive accuracy. The training set will be what the model sees, and the test set will be used to evaluate the model.

set.seed(2018)
iris.sample = read.csv("https://tinyurl.com/yaq7grsh", header=TRUE)

train_indices = 1:round(0.80 * nrow(iris.sample))
train = as.data.frame(iris.sample[train_indices, ])
test = as.data.frame(iris.sample[-train_indices,])
test = data.frame(lapply(test, jitter, amount=.5))

Below is shown the correlation matrix for the predictors, which indicates a high level of correlation.

High levels of correlation between predictors, or collinearity results in inflated variance and potentially incorrect estimates of \(\hat{\beta}\). The level of variance inflation in a model can be checked using the vif function of the car package. Let’s consider a model predicting the sepal width.

These values are far less than ideal, and warrant the use of a regularized regression. Ridge regression can be used with the glmnet package. The glmnet package allows the fitting of ridge, lasso, and elastic net (a combination of ridge and lasso) models. It requires the model be given to it as a matrix of predictors and a vector of dependent variable values, but this is easy to get using the model.matrix function.

Next, I will fit the model using k-fold cross-validation. This splits the data into \(\textit{k}\) subsets in order to fine tune the lambda parameter, which controls the amount of shrinkage. The returned model is then the model fit to the whole data set

x.mat =model.matrix(Sepal.Width  ~ ., train)[,-1]
y.var = train$Sepal.Width 

cv.ridge = cv.glmnet(x.mat,y.var, alpha=0, nlambda=100, lambda.min.ratio=0.0001) #Set alpha to zero for ridge regression

A plot showing the the candidate values of the tuning parameter \(\lambda\) can be produced. The lower dashed line is the minimum value of \(\lambda\) that gives minimum mean cross-validated mean squared error. The other dashed line is the maximum value of \(\lambda\) that gives a regularized model such that error is within one standard error of the minimum. Either one of these values can be used. The lower dashed line is stored as lambda.min in cv.ridge object. This will give the most regularized model if used to produce coefficients. The upper dashed line is stored as lambda.1se and will result in a slightly less regularized model. I prefer lambda.min.

plot(cv.ridge)

The coefficients for the best value of lambda can then be retrieved and made into a table. In a separate tab the results for the ordinary least squares model fit can be seen.

Ridge regression really shines when it comes to out of sample predictive error. Ordinary least squares can potentially over-fit the data, such that when presented with new cases it may not perform well. Ridge regression has the advantage that by under-fitting the model it generalizes to new samples more readily.

x.test =model.matrix(Sepal.Width ~ ., test)[,-1]
y.test = test$Sepal.Width

ridge.train = data.frame(predict(cv.ridge, x.test, s="lambda.1se"), y.test)
colnames(ridge.train)=c("predicted", "actual")
ridge.train$residual = ridge.train$actual-ridge.train$predicted

ols.train = data.frame(predict(lm(Sepal.Width ~ ., train), newdata=as.data.frame(x.test)), y.test)
colnames(ols.train)=c("predicted", "actual")
ols.train$residual = ols.train$actual-ols.train$predicted

Below is a comparison of the predictive performance of both models by several common measures calculated from the code above.

Modeling Heteroskedastic Variance

One benefit of Bayesian model building is constructing a model where the variance of an independent variable \(y\) depends on the value of \(x_i\) is quite easy. This falls under the general category of distributional Bayesian modeling (Umlauf and Kneib 2018). Distributional modeling attempts to estimate model parameters such as the variance from the data, rather than treating them as nuisance parameters, and is a useful alternative to quantile regression in situations like this (Koenker, Leorato, and Peracchi 2013). The brms package can easily implement this model.

formula = bf(Neuroticism ~ Age, sigma ~ Age)
prior = get_prior(formula, data)
prior$prior[1] = "student_t(3, 0, 10)"
fit.hetvar = brm(formula, data=data, prior=prior, iter = 4000, cores = 4)

prior = get_prior(Neuroticism ~ Age, data)
prior$prior[1] = "student_t(3,0,10)"
fit.normal = brm(Neuroticism ~ Age, data, prior=prior, iter = 4000, cores = 4)

Student T Model

A Bayesian alternative to MM-Estimators or Least Trimmed Squares is to use the Student’s t model as the family function in a regression model (Gelman et al. 2014). This assumes that the observations \(y_i\) are drawn from a t-distribution, which has longer tails stretching past the center of the distribution’s mass. This is useful to represent situations where the bulk of the data are from an essentially gaussian distribution but contain some data points a distance away from the rest. Unlike MM-Estimators or Least Trimmed Squares, this directly models the outlier-generating process. The model gives the most weight to the bulk of the distribution, giving less attention to the tails. This results in a lower standard error and decreases the chances that rarely occuring observations mislead you about the relationship of the predictors to the dependent variable. Furthermore, the normality parameter is estimated from the data to indicate how “robust” the model needs to be.

For example, consider the ldl variable from the corrupted diabetes data set. The outlying observations create long, thin, tails of the distribution.

The brms package can easily accomodate this situation in a Student’s t regression model. The results can be compared to the OLS, LTS, and MM-Estimator results from before.

prior = get_prior(ldl ~ hdl, data = diabetes, family=student)
prior$prior[1] = "student_t(3,0,5)"
fit.robust = brm(ldl ~ hdl, diabetes, prior=prior, family=student, iter = 4000, cores = 4)

Skew Normal Model

Data may sometimes not quite follow a normal distribution, but lack outliers. When there are more observations on one side of the Gaussian curve than another the result is a skewed normal distribution (Arellano-Valle et al. 2008; Alhamide, Ibrahim, and Alodat 2015). If the observations are skewed badly enough one might wish to account for this skewness in a regression model. This is fully achievable by using a skew normal probability density function in lieu of a normal p.d.f. for the likelihood.

For example, consider the bmi variable in the diabetes data set used in previous examples.

The distribution is clearly skewed. For an example, let’s fit a model of bmi with age and hdl as predictors.

prior = get_prior(bmi ~ age + hdl, data = diabetes, family=skew_normal)
prior$prior[1] = "student_t(7,0,5)"
prior$prior[2] = "student_t(3,0,10)"
fit.skew = brm(bmi ~ age + hdl, diabetes, family=skew_normal, iter = 4000, cores = 4)

prior = get_prior(bmi ~ age + hdl, data = diabetes)
prior$prior[1] = "student_t(3,0,10)"
fit.gauss= brm(bmi ~ age + hdl, data = diabetes, iter = 4000, cores = 4)

When compared to the Gaussian model, the skew normal model’s coefficients are of moderately different strengths. If the sample size were smaller or the coefficients smaller in magnitude, these differences could make the difference between declaring a coefficient to be credibly non-zero and not making such a claim. If the data in the dependent variable are skewed, the skew normal model is the safest bet. If the skewness is not enough to matter the coefficients will come out the same as a typical normal model.

Regularized Regresssion

The ridge penalty is equivalent to setting a normally distributed (Gaussian) prior of \(Normal~\sim~(\mu=0, ~ \sigma = \frac{1}{2\cdot p}\cdot I(\theta)~)\), where \(I\) is the expected fisher information.

The Fisher information is estimate is greater than or equal to the reciprocal of the fisher information (a rule called the Cramer-Rao bound). The standard deviation of the estimate (or prior in this case) is then, at least, the square root of the reciprocal of the fisher information.

\[ \text{VAR}({\hat\theta}) \geq \frac{1}{I(\theta)} \\ \text{sd}({\hat\theta}) \geq \sqrt{\frac{1}{I(\theta)}} \]

Standardizing variables to a mean of zero and standard deviation of 1 makes it an easy matter to decide on the prior standard deviation for the coefficients by simply multiplying \(1/2p\) by 1 for the prior sd; \(Normal~\sim~(\mu=0, ~ \sigma = \frac{1}{2\cdot p}\cdot 1)\).

In this example there are three predictors, making the prior equal to \(Normal~\sim~(0, ~ \frac{1}{6})\). Now that the prior has been obtained the model can be fit in brms.

prior = get_prior(Sepal.Width ~ ., train)
prior$prior[1] = "normal(0, 0.1666667)"
fit.ridge = brm(Sepal.Width ~ ., train, prior=prior, iter=4000, cores = 4)

## Compiling the C++ model

## Start sampling

Notice that the estimates are nearly identical to those obtained using lambda$min in the glmnet ridge model. This is no coincidence. The prior was designed to reflect the properties of the ridge estimator after all! So why use this over glmnet? The benefit of the Bayesian approach is that full posterior distributions are obtained for each parameter, allowing statistical inference about which parameters are credibly non-zero. There is no widely accepted standard for statistical inference from a ridge regression model fit in glmnet. Classical frequentist hypothesis tests largely rely on unbiased estimators to obtain test statistics that follow specific sampling distributions, but Bayesian inference proceeds by updating prior beliefs with the data to obtain a posterior inference.

The model’s predictions can also be compared against a test sample to evaluate its performance, just as with any other regression model. Below is the mean square error, which is comparable (if not slightly better) than what was obtained using glmnet.

pred = as.data.frame(predict(fit.ridge, newdata=x.test))
mean((pred$Estimate-y.test)^2)

## [1] 0.6780166