How to tell which variable is more meaningful when modeling the relationship between several predictors and outcome variable?

Question

I'm facing a problem in which I need to figure out two things:

• which predictor, out of several relevant ones, is the most meaningful one in its effect/predictive power over a predicted variable.
• the order of meaningfulness (from most meaningful to least) across those different predictors.

As I do not have an a-priory hypothesis about this investigation, I thought I should be doing some sort of multiple regression analysis. Then, perhaps I should be extracting the terms from the model and see which one is the most meaningful. I already know that going by p-value isn't the right way. Then what is?

Example

Let's say that I want to investigate which factors affect the well-being of city residents. I sample people (residents) from both New York City and San Francisco, and ask them to rate:

1. Their general satisfaction from their city.
2. How clean their city (in their opinion)
3. How good the level of education in schools
4. How good the public transportation.

The way I see this, there are 3 relevant variables here (cleanliness, education, and transportation) that may be related to overall satisfaction. I want to model this relationship, and conclude how different NYC from SF is. For example, when it comes to overall satisfaction, is the impact of education > transportation > cleanliness in NYC, whereas in SF transportation > education > cleanliness?

Here's some toy data to demonstrate.

my_df <- structure(list(location = c("sf", "nyc", "nyc", "sf", "nyc", 
                                  "nyc", "nyc", "nyc", "nyc", "sf", "nyc", "sf", "sf", "sf", "nyc", 
                                  "sf", "sf", "nyc", "nyc", "nyc", "sf", "sf", "sf", "sf", "sf", 
                                  "nyc", "sf", "sf", "nyc", "sf", "nyc", "nyc", "nyc", "sf", "nyc", 
                                  "nyc", "nyc", "sf", "nyc", "sf", "sf", "nyc", "nyc", "nyc", "nyc", 
                                  "nyc", "nyc", "nyc", "nyc", "sf", "nyc", "nyc", "sf", "sf", "nyc", 
                                  "nyc", "nyc", "nyc", "sf", "sf", "nyc", "sf", "nyc", "nyc", "sf", 
                                  "nyc", "sf", "sf", "nyc", "nyc", "nyc", "nyc", "sf", "nyc", "nyc", 
                                  "nyc", "sf", "sf", "nyc", "nyc", "nyc", "nyc", "sf", "sf", "nyc", 
                                  "nyc", "nyc", "sf", "sf", "sf", "nyc", "sf", "sf", "sf", "nyc", 
                                  "nyc", "nyc", "nyc", "nyc", "nyc"), 
                     satisfied = c(5, 1, 7, 5, 
                                   7, 1, 5, 5, 5, 7, 7, 4, 1, 3, 5, 6, 7, 7, 6, 4, 4, 5, 6, 5, 5, 
                                   7, 5, 6, 5, 4, 7, 7, 5, 5, 4, 7, 7, 5, 6, 6, 3, 6, 5, 7, 5, 7, 
                                   6, 5, 4, 3, 6, 5, 7, 3, 5, 5, 7, 5, 6, 7, 7, 7, 7, 5, 4, 7, 6, 
                                   7, 7, 6, 6, 6, 5, 7, 5, 4, 6, 4, 7, 5, 6, 6, 5, 5, 6, 7, 6, 5, 
                                   1, 5, 2, 7, 7, 7, 7, 7, 1, 3, 7, 7), 
                     clean = c(4, 1, 
                                         7, 3, 4, 1, 6, 6, 7, 4, 5, 1, 1, 1, 4, 6, 6, 1, 6, 1, 4, 2, 2, 
                                         7, 3, 5, 2, 4, 1, 1, 4, 6, 3, 5, 1, 4, 5, 2, 5, 5, 4, 5, 4, 7, 
                                         3, 6, 5, 4, 5, 4, 5, 4, 5, 1, 5, 2, 6, 5, 7, 6, 3, 7, 5, 6, 4, 
                                         6, 6, 5, 5, 5, 4, 1, 4, 4, 5, 1, 3, 1, 2, 2, 6, 4, 3, 6, 7, 7, 
                                         5, 2, 4, 1, 3, 1, 5, 3, 5, 5, 1, 1, 5, 6), 
                     edu = c(5, 
                                             1, 7, 4, 4, 1, 6, 6, 6, 4, 5, 4, 4, 3, 3, 5, 4, 1, 1, 3, 5, 6, 
                                             5, 5, 3, 6, 2, 4, 4, 4, 6, 3, 4, 7, 1, 4, 7, 5, 6, 5, 5, 5, 4, 
                                             7, 3, 7, 6, 5, 5, 4, 5, 3, 4, 4, 4, 7, 5, 4, 6, 6, 4, 7, 4, 2, 
                                             5, 6, 6, 7, 6, 7, 3, 3, 2, 6, 6, 2, 5, 3, 6, 5, 6, 4, 4, 5, 6, 
                                             7, 3, 3, 4, 5, 4, 1, 3, 4, 4, 6, 5, 1, 4, 6), 
                     transportation = c(1, 
                                              1, 7, 5, 7, 1, 4, 6, 6, 6, 6, 1, 1, 1, 5, 5, 4, 7, 6, 6, 7, 5, 
                                              2, 7, 3, 6, 1, 4, 7, 5, 6, 7, 4, 3, 2, 6, 4, 2, 6, 5, 4, 7, 6, 
                                              7, 3, 7, 4, 4, 5, 4, 6, 3, 5, 2, 7, 3, 7, 7, 7, 6, 7, 7, 7, 5, 
                                              3, 5, 4, 7, 6, 6, 4, 2, 4, 4, 5, 6, 5, 2, 6, 2, 6, 6, 3, 4, 7, 
                                              7, 7, 4, 5, 4, 5, 3, 7, 5, 7, 7, 7, 1, 6, 6)), 
                row.names = c(NA, -100L), 
                class = c("tbl_df", "tbl", "data.frame"))


library(magrittr)
library(effectsize)

my_df %>%
  lm(satisfied ~ clean*location + edu*location + transportation*location, data = .) %>%
  effectsize() %>%
  plot()
#> Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please
#> use `guide = "none"` instead.

^{Created on 2021-08-15 by the reprex package (v2.0.0)}

To compare between cities, I added interaction terms between each edu/clean/transportation variable and location variable. Then I used effectsize::effectsize() to get the estimates from the model. But what can I conclude from those estimates?

If I completely got this wrong, please advise what other path I should take for tackling this problem.

Thanks!

Answer 1

I believe you're looking for a metric of (relative) variable importance (see also this thread). Many available methods rely on the decomposition of the $R^2$ to assign ranks or relative importance to each predictor in a multiple linear regression model. A certain approach in this family is better known under the term "Dominance analysis" (see Azen et al. 2003). Azen et al. (2003) also discuss other measures of importance such as importance based on regression coefficients, based on correlations of importance based on a combination of coefficients and correlations. A general good overview of techniques based on variance decomposition can be found in the paper of Grömping (2012). These techniques are implemented in the R packages relaimpo, domir and yhat.

Here, I'm going to illustrate a method that is model-agnostic (i.e. it can be applied to a variety of model types) and has intuitive appeal: Variable importance based on permutation. The idea is very simple:

1. Decide on a performance metric that is important to you. Examples include: Root mean square error (RMSE), mean absolute error (MAE), $R^2$ etc. This also is somewhat dependent in the model type.
2. Calculate the metric on your dataset, $M_{orig}$. This serves as baseline performance metric.
3. For $i = 1, 2, \ldots, j$:
   (a) Permute the values of the predictor $X_i$ in the data set.
   (b) Recompute the metric on the permuted data and call it $M_{perm}$.
   (c) Record the difference from baseline using $imp(X_i)=M_{perm} - M_{orig}$.

Do this repeatedly, say 1000 times, and take the average of the importance values. Intuitively, the permutations break the relationship between the predictor $X_i$ and the outcome. The larger the change in the performance metric, the higher the predictors' importance. More information can be found in this chapter of an online book by Christoph Molnar.

The R package vip implements this procedure (see the documentation (PDF) for more information). The following code applies the idea to your dataset. I chose the $R^2$ and the mean absolute error (MAE) as performance metrics and permute each predictor 1000 times:

library(vip)
# The model
mod <- lm(satisfied ~ clean*location + edu*location + transportation*location, data = my_df)

# Calculate permutation-based importance with r-squared as metric
set.seed(142857) # For reproducibility
p_r2 <- vip::vi(mod, method = "permute", target = "satisfied", metric = "rsquared", pred_wrapper = predict, nsim = 1000)
p_r2

  Variable       Importance  StDev
  <chr>               <dbl>  <dbl>
1 transportation     0.198  0.0492
2 clean              0.177  0.0465
3 edu                0.0462 0.0237
4 location           0.0449 0.0250

# Calculate permutation-based importance with mae as metric
p_mae <- vip::vi(mod, method = "permute", target = "satisfied", metric = "mae", pred_wrapper = predict, nsim = 1000)
p_mae

  Variable       Importance  StDev
  <chr>               <dbl>  <dbl>
1 transportation     0.166  0.0413
2 clean              0.144  0.0400
3 location           0.0396 0.0214
4 edu                0.0368 0.0219

According to the $R^2$, permuting transportation leads to the largest change in $R^2$, followed by clean. Using the mean absolute error shows a similar ordering with transportation and clean being most important while location and edu being least important.

References

Azen R, Budescu DV (2003): The Dominance Analysis Approach for Comparing Predictors in Multiple Regression. Psychological Methods 8:2, 129-148. (link)

Grömping U (2012): Estimators of relative importance in linear regression based on variance decomposition. Am Stat 61:2, 139-147. (link)

Answer 2

The following answer is some sort of an ignorant attempt to use {rms} package, following @JTH's suggestion. I got to say that this is the first time I'm using this package, and I have very minimal understanding of what I'm doing. Hence, I ask that anybody who can -- please provide feedback!

I've followed the procedure described in this chapter.

my_df <- structure(list(location = c("sf", "nyc", "nyc", "sf", "nyc", 
                                     "nyc", "nyc", "nyc", "nyc", "sf", "nyc", "sf", "sf", "sf", "nyc", 
                                     "sf", "sf", "nyc", "nyc", "nyc", "sf", "sf", "sf", "sf", "sf", 
                                     "nyc", "sf", "sf", "nyc", "sf", "nyc", "nyc", "nyc", "sf", "nyc", 
                                     "nyc", "nyc", "sf", "nyc", "sf", "sf", "nyc", "nyc", "nyc", "nyc", 
                                     "nyc", "nyc", "nyc", "nyc", "sf", "nyc", "nyc", "sf", "sf", "nyc", 
                                     "nyc", "nyc", "nyc", "sf", "sf", "nyc", "sf", "nyc", "nyc", "sf", 
                                     "nyc", "sf", "sf", "nyc", "nyc", "nyc", "nyc", "sf", "nyc", "nyc", 
                                     "nyc", "sf", "sf", "nyc", "nyc", "nyc", "nyc", "sf", "sf", "nyc", 
                                     "nyc", "nyc", "sf", "sf", "sf", "nyc", "sf", "sf", "sf", "nyc", 
                                     "nyc", "nyc", "nyc", "nyc", "nyc"), 
                        satisfied = c(5, 1, 7, 5, 
                                      7, 1, 5, 5, 5, 7, 7, 4, 1, 3, 5, 6, 7, 7, 6, 4, 4, 5, 6, 5, 5, 
                                      7, 5, 6, 5, 4, 7, 7, 5, 5, 4, 7, 7, 5, 6, 6, 3, 6, 5, 7, 5, 7, 
                                      6, 5, 4, 3, 6, 5, 7, 3, 5, 5, 7, 5, 6, 7, 7, 7, 7, 5, 4, 7, 6, 
                                      7, 7, 6, 6, 6, 5, 7, 5, 4, 6, 4, 7, 5, 6, 6, 5, 5, 6, 7, 6, 5, 
                                      1, 5, 2, 7, 7, 7, 7, 7, 1, 3, 7, 7), 
                        clean = c(4, 1, 
                                  7, 3, 4, 1, 6, 6, 7, 4, 5, 1, 1, 1, 4, 6, 6, 1, 6, 1, 4, 2, 2, 
                                  7, 3, 5, 2, 4, 1, 1, 4, 6, 3, 5, 1, 4, 5, 2, 5, 5, 4, 5, 4, 7, 
                                  3, 6, 5, 4, 5, 4, 5, 4, 5, 1, 5, 2, 6, 5, 7, 6, 3, 7, 5, 6, 4, 
                                  6, 6, 5, 5, 5, 4, 1, 4, 4, 5, 1, 3, 1, 2, 2, 6, 4, 3, 6, 7, 7, 
                                  5, 2, 4, 1, 3, 1, 5, 3, 5, 5, 1, 1, 5, 6), 
                        edu = c(5, 
                                1, 7, 4, 4, 1, 6, 6, 6, 4, 5, 4, 4, 3, 3, 5, 4, 1, 1, 3, 5, 6, 
                                5, 5, 3, 6, 2, 4, 4, 4, 6, 3, 4, 7, 1, 4, 7, 5, 6, 5, 5, 5, 4, 
                                7, 3, 7, 6, 5, 5, 4, 5, 3, 4, 4, 4, 7, 5, 4, 6, 6, 4, 7, 4, 2, 
                                5, 6, 6, 7, 6, 7, 3, 3, 2, 6, 6, 2, 5, 3, 6, 5, 6, 4, 4, 5, 6, 
                                7, 3, 3, 4, 5, 4, 1, 3, 4, 4, 6, 5, 1, 4, 6), 
                        transportation = c(1, 
                                           1, 7, 5, 7, 1, 4, 6, 6, 6, 6, 1, 1, 1, 5, 5, 4, 7, 6, 6, 7, 5, 
                                           2, 7, 3, 6, 1, 4, 7, 5, 6, 7, 4, 3, 2, 6, 4, 2, 6, 5, 4, 7, 6, 
                                           7, 3, 7, 4, 4, 5, 4, 6, 3, 5, 2, 7, 3, 7, 7, 7, 6, 7, 7, 7, 5, 
                                           3, 5, 4, 7, 6, 6, 4, 2, 4, 4, 5, 6, 5, 2, 6, 2, 6, 6, 3, 4, 7, 
                                           7, 7, 4, 5, 4, 5, 3, 7, 5, 7, 7, 7, 1, 6, 6)), 
                   row.names = c(NA, -100L), 
                   class = c("tbl_df", "tbl", "data.frame"))

library(rms, warn.conflicts = FALSE)
#> Loading required package: Hmisc
#> Loading required package: lattice
#> Loading required package: survival
#> Loading required package: Formula
#> Loading required package: ggplot2
#> 
#> Attaching package: 'Hmisc'
#> The following objects are masked from 'package:base':
#> 
#>     format.pval, units
#> Loading required package: SparseM
#> 
#> Attaching package: 'SparseM'
#> The following object is masked from 'package:base':
#> 
#>     backsolve

model_fit <- rms::ols(satisfied ~ clean*location + edu*location + transportation*location, data = my_df)

my_datadist <- rms::datadist(my_df)     ## apparently, we need these two lines 
options(datadist = "my_datadist")  ## otherwise we get an error with `summary(model_fit)`
                                   ## I learned it from here: https://stackoverflow.com/a/41378930/6105259 

plot(summary(model_fit))

^{Created on 2021-08-16 by the reprex package (v2.0.0)}

---

As far as I understand this output, we see the effect each predictor carries on the outcome variable, with 95% CI. Thus, for example, we can conclude that clean has a greater effect over satisfied than edu has.

• Does this make sense?
• Originally, I was interested in the interaction between each edu/clean/transportation and location, as I want to learn how the relationships between predictors and outcome change between cities. But here, as far as I can see from this output, the interaction isn't reflected in terms of effect on satisfied.

---

UPDATE

---

Following @JTH's comment, I'm adding another plot:

plot(anova(model_fit))

The info underlying this plot is here:

library(dplyr, warn.conflicts = FALSE)
library(tibble)

anova(model_fit) %>% 
  as_tibble(rownames = "factor") %>%
  mutate(across(3:6, round, 4))
#> # A tibble: 14 x 6
#>    factor                          d.f.     `Partial SS` MS      F       P      
#>    <chr>                           <anov.r> <anov.rms>   <anov.> <anov.> <anov.>
#>  1 "clean  (Factor+Higher Order F~  2        11.8209      5.9105 3.4587  0.0356 
#>  2 " All Interactions"              1         0.0000      0.0000 0.0000  0.9969 
#>  3 "location  (Factor+Higher Orde~  4         4.8955      1.2239 0.7162  0.5830 
#>  4 " All Interactions"              3         4.8821      1.6274 0.9523  0.4188 
#>  5 "edu  (Factor+Higher Order Fac~  2         4.0844      2.0422 1.1951  0.3073 
#>  6 " All Interactions"              1         3.1038      3.1038 1.8163  0.1811 
#>  7 "transportation  (Factor+Highe~  2        15.3207      7.6604 4.4828  0.0139 
#>  8 " All Interactions"              1         0.0476      0.0476 0.0279  0.8678 
#>  9 "clean * location  (Factor+Hig~  1         0.0000      0.0000 0.0000  0.9969 
#> 10 "location * edu  (Factor+Highe~  1         3.1038      3.1038 1.8163  0.1811 
#> 11 "location * transportation  (F~  1         0.0476      0.0476 0.0279  0.8678 
#> 12 "TOTAL INTERACTION"              3         4.8821      1.6274 0.9523  0.4188 
#> 13 "TOTAL"                          7        90.9761     12.9966 7.6055  0.0000 
#> 14 "ERROR"                         92       157.2139      1.7088     NA      NA

---

UPDATE 2

---

Addressing @EdM's comment, here is a print of anova(model_fit) without converting it to a tibble.

anova(model_fit) %>% 
  round(., 4)
#>                 Analysis of Variance          Response: satisfied 
#> 
#>  Factor                                                   d.f. Partial SS
#>  clean  (Factor+Higher Order Factors)                      2    11.8209  
#>   All Interactions                                         1     0.0000  
#>  location  (Factor+Higher Order Factors)                   4     4.8955  
#>   All Interactions                                         3     4.8821  
#>  edu  (Factor+Higher Order Factors)                        2     4.0844  
#>   All Interactions                                         1     3.1038  
#>  transportation  (Factor+Higher Order Factors)             2    15.3207  
#>   All Interactions                                         1     0.0476  
#>  clean * location  (Factor+Higher Order Factors)           1     0.0000  
#>  location * edu  (Factor+Higher Order Factors)             1     3.1038  
#>  location * transportation  (Factor+Higher Order Factors)  1     0.0476  
#>  TOTAL INTERACTION                                         3     4.8821  
#>  REGRESSION                                                7    90.9761  
#>  ERROR                                                    92   157.2139  
#>  MS      F    P     
#>   5.9105 3.46 0.0356
#>   0.0000 0.00 0.9969
#>   1.2239 0.72 0.5830
#>   1.6274 0.95 0.4188
#>   2.0422 1.20 0.3073
#>   3.1038 1.82 0.1811
#>   7.6604 4.48 0.0139
#>   0.0476 0.03 0.8678
#>   0.0000 0.00 0.9969
#>   3.1038 1.82 0.1811
#>   0.0476 0.03 0.8678
#>   1.6274 0.95 0.4188
#>  12.9966 7.61 <.0001
#>   1.7088

Answer 3

As all your variables are numeric and limited to the same range of values, directly comparing the absolute values of the coefficients (as you did in your answer) is one possible way to estimate effect size.

It might be, however, that for one variable the answers are only in the range 1-3, while they are for another variable in the range 1-5. If both variables are equally important, the second will have a smaller coefficient due to its wider range. To overcome this problem, you could compute standardized coefficients, aka "beta coefficients". There are convenience functions for their computation in addon packages, but with base R you can achieve the same by standardizing your data before doing the model fit:

my_df.scaled <- scale(my_df)

In addition, you should also check, whether all coefficients are significantly different from zero.

Answer 4

I'd look at the residuals for determining the effectiveness of a model. Plots of residuals versus other quantities are used to find failures of assumptions. The most common plot, especially useful in simple regression, is the plot of residuals versus the fitted values. A null plot would indicate no failure of assumptions. Curvature might indicate the fitted mean function is incorrect. Residuals that seem to increase or decrease in average magnitude with the fitted values might indicate nonconstant residual variance. A few relatively large residuals may be indicative of outliers, cases for which the model is somehow inappropriate.

Assumptions of a linear model are:

1. Linear relationship: There exists a linear relationship between the independent variable, x, and the dependent variable, y.

2. Independence: The residuals are independent. In particular, there is no correlation between consecutive residuals.

3. Homoscedasticity: The residuals have constant variance at every level of x.

4. Normality: The residuals of the model are normally distributed.

If one or more of these assumptions are violated, then the results of our linear regression may be unreliable or even misleading.

Using the given data, I build a simple linear regression model given below.

# required libraries
library(caret)
library(ggplot2)
library(magrittr)

# split the train dataset into train and test set
set.seed(2021)
index <- createDataPartition(my_df$satisfied, p = 0.7, list = FALSE)
df_train <- my_df[index, ]
df_test  <- my_df[-index, ]

# Model building 
lm_model <- lm(satisfied ~., data = my_df)
# Make predictions and compute the R2, RMSE and MAE
predictions <- lm_model %>% predict(df_test)
data.frame( R2 = R2(predictions, df_test$satisfied),
            RMSE = RMSE(predictions, df_test$satisfied),
            MAE = MAE(predictions, df_test$satisfied))
         R2      RMSE       MAE
1 0.4513587 0.9956456 0.8224509
# Residuals = Observed - Predicted
# compute residuals
residualVals <- df_test$satisfied - predictions
df.1 <- data.frame(df_test$satisfied, predictions, 
                   residualVals)
colnames(df.1)<- c("observed","predicted","residuals")
head(df.1)
  observed predicted   residuals
1        5  4.443428  0.55657221
2        7  6.930498  0.06950225
3        4  3.626674  0.37332568
4        3  3.538958 -0.53895795
5        5  5.083349 -0.08334872
6        7  6.097199  0.90280068
ggplot(data = df.1, aes(x=predicted, y=residuals))+
  geom_point()+
  xlab("Predicted values for satsified")+
  ylab("Residuals")+
  ggtitle("Residual plot")+
  theme_bw()

Discussion

To understand the strength/weakness of a model, relying on a single metric is problematic. Visualization of model fit, particularly residual plots in the context of linear regression model, are critical to understanding whether the model is fit for purpose.

When the outcome is a number, the most common method for characterizing a model’s predictive capabilities is to use the root mean squared error (RMSE). This metric is a function of the model residuals, which are the observed values minus the model predictions. The mean squared error (MSE) is calculated by squaring the residuals and summing them. The RMSE is then calculated by taking the square root of the MSE so that it is in the same units as the original data. The value is usually interpreted as either how far (on average) the residuals are from zero or as the average distance between the observed values and the model predictions. Another common metric is the coefficient of determination, commonly written as R2. This value can be interpreted as the proportion of the information in the data that is explained by the model. Thus, an R2 value of 0.45 in the above model, implies that the model can explain less than half of the variation in the outcome/dependent/response variable, satisfied variation. Simply put, the above model is not good. It should be noted that R2 is a measure of correlation and not accuracy.