# Methods or processes for discovering the major reasons for linear regression coefficients?

Other than inspecting data at random, what are some statistical or graphical processes or methods for determining the major reasons for coefficients in linear regressions?

### Example

Suppose we have the iris dataset, and regress

Sepal.Length = Sepal.Width + Petal.Length + Petal.Width + Speciesversicolor_dummy + Speciesvirginica_dummy


giving

Coefficients:

(Intercept)              Sepal.Width             Petal.Length              Petal.Width
2.1713                   0.4959                   0.8292                  -0.3152

Speciesversicolor_dummy   Speciesvirginica_dummy
-0.7236                  -1.0235


When inspecting the dummy variables, the result looks odd at first glance, since the omitted Species, setosa, appears to increase the prediction by 2.1713 (given by the intercept), and the other two species (versicolor and virginica) have negative coefficients, despite their average Sepal.Length being much higher than that of setosa.

  Species    ave_sep_len
1 setosa            5.01
2 versicolor        5.94
3 virginica         6.59


This seemed odd, until I inspected the other variables, which showed the full picture:

  Species    ave_sep_len ave_sep_wid ave_pet_len ave_pet_wid
1 setosa            5.01        3.43        1.46       0.246
2 versicolor        5.94        2.77        4.26       1.33
3 virginica         6.59        2.97        5.55       2.03


and it becomes evident that the massive difference in ave_pet_len is driving a huge amount of the difference in the estimates, so much so that the dummy variables are pushed into having counter-intuitive sign.

### Question

I learned of the main driver(s) of the coefficients through manual inspection, and the dataset was sufficiently small as to allow this, but on larger datasets this may not be viable.

What methods (statistical, graphical or otherwise) are commonly used for larger datasets?

### Note

Here is some R code to reproduce if desired

library(dplyr)

# With Species
Species_Dummies <- model.matrix(~ iris$$Species - 1) %>% as.data.frame %>% select(-1) %>% # Remove first column and let whatever it was be the 'base' model colnames<-(str_split(names(.), "\\$$") %>% sapply(function(x) { paste0(x[2], "_dummy") } )) # Give columns nice names

(fit <- lm(Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width + Speciesversicolor_dummy + Speciesvirginica_dummy, data = cbind(iris, Species_Dummies)))

# Coefficients:
#   (Intercept)              Sepal.Width             Petal.Length              Petal.Width
# 2.1713                   0.4959                   0.8292                  -0.3152
#
# Speciesversicolor_dummy   Speciesvirginica_dummy
# -0.7236                  -1.0235

iris %>%
group_by(Species) %>%
summarise(ave_sep_len = mean(Sepal.Length),
ave_sep_wid  = mean(Sepal.Width),
ave_pet_len  = mean(Petal.Length),
ave_pet_wid  = mean(Petal.Width))
#   Species    ave_sep_len ave_sep_wid ave_pet_len ave_pet_wid
# 1 setosa            5.01        3.43        1.46       0.246
# 2 versicolor        5.94        2.77        4.26       1.33
# 3 virginica         6.59        2.97        5.55       2.03



For visualizing, since you are going to fit a linear model, you can try calculating the correlations and visualizing on a heatmap:

library(pheatmap)
pheatmap(cor(model.matrix(Sepal.Length ~ .,data=iris)[,-1]))


The above should be ok for moderately sized datasets. You can also try corrplot but I find it less intuitive.

In terms of diagnostics, what you are looking at is multicollinearity, and they are a few discussions for example this on CV and this from SO. Since you are using R, I recommend also Faraway 2002's section on Scale Changes, Principal Components and Collinearity, below I just one of them, variance-inflation factor, which tells you how much the variance is inflated for each coefficient. If the coefficient is not correlated to any others (i.e orthogonal), it's vif will be close to 1, and the norm is something >30 is considered problematic, >10 is of concern:

X = model.matrix(Sepal.Length ~ .,data=iris)
library(faraway)
vif(X)
(Intercept)       Sepal.Width      Petal.Length       Petal.Width
2.002515          2.227466         23.161648         21.021401
Speciesversicolor  Speciesvirginica
20.423390         39.434378


The other often mention is conditional index, in R it is available in the perturb. I am not so familiar with that so you can check it out in this post and also the recommended reading.

• Great answer. also, I found another library that does a similar thing to pheatmap, corrplot Apr 18, 2020 at 3:00

Linear Least-Squares is not robust to outliers, which are even harder to detect in large data sets. Per Wikipedia on Robust Regression, to quote:

Certain widely used methods of regression, such as ordinary least-squares, have favourable properties if their underlying assumptions are true, but can give misleading results if those assumptions are not true; thus ordinary least squares is said to be not robust to violations of its assumptions....

In particular, least squares estimates for regression models are highly sensitive to outliers.

A more robust and iterative methodology is, for example, Least-Absolute-Deviations (LAD). Parameters can be derived by iteratively re-weighted least-squares methodology to find the smallest average of the sum of the absolute difference between each fitted and observed point. LAD is not a mean perspective, but essentially a median focus (the one-parameter model solution is, in fact, the median). In practice, with a two-parameter LAD model (at least per my experience), one observes that the solution consists of a line that passes through at least two reported points. So, this does mirror the process of graphically looking for 'best' lines passing through points.

A technical point per Wikipedia, to quote:

The least absolute deviations estimate also arises as the maximum likelihood estimate if the errors have a Laplace distribution.

If it is apparent (by testing, for example) that the error terms are not from the Laplace distribution (perhaps more like log-normal), I would suggest also exploring data transformations.