Determining statistical significance of linear regression coefficient in the presence of multicollinearity

Question

Suppose I have a bunch of cities with different population sizes, and I wanted to see if there was a positive linear relationship between the number of liquor stores in a city and the number of DUIs. Where I'm determining whether this relationship is significant or not based on a t-test of the estimated regression coefficient.

Now clearly the pop. size of a city is going to be positively correlated with both the number of DUIs as well as the number of liquor stores. Thus if I run a simple linear regression on just liquor stores and see if its regression coefficient is statistically significant, I will likely run into a problem of multicollinearity, and over-estimate the effect of liquor stores on DUIs.

Which of the two methods should I use to correct for this?

1. I should divide the number of liquor stores in the city by its population in order to get a liquor store per capita value and then regress on that.

2. I should regress on both liquor stores and size, and then look to see if the liquor store coefficient is significant when controlling for size.

3. Some other method?

I honestly can't decide which seems more sensible. I vacillate between them, depending on which one I think about I'm able to convince myself that that's the right way.

On the one hand liquor stores per capita seems like the right variable to use, since DUIs are committed by individuals, but that doesn't seem very statistically rigorous. On the other hand, controlling for size seems statistically rigorous, but rather indirect. Furthermore, if I rescale after computing the liquor stores per capita variable, I get very similar regression coefficients between the two methods, but method 1 produces a smaller p-value.

Answer 1

I would regress the "DUI per capita" (Y) on "liquer stores per capita" (X) and "population size" (Z). This way your Y reflects the propensity to drunk driving of urban people, while X is the population characteristic of a given city. Z is a control variable just in case if there's size effect on Y. I don't think you are going to see multicollinearity issue in this setup.

This setup is more interesting than your model 1. Here, your base is to assume that the number of DUIs is proportional to population, while $\beta_Z$ would capture nonlinearity, e.g. people in larger cities are more prone to drunk driving. Also X reflects cultural and legal environment directly, already adjusted to size. You may end up with roughly the same X for cities of different sizes in Sough. This also allows you introduce other control variables such as Red/Blue state, Coastal/Continental etc.

Answer 2

If estimating your model with ordinary least squares, your second regression is rather problematic.

And you may want to think about how the variance of your error term varies with city size.

Regression (2) is equivalent to your regression (1) where observations are weighted by the square of the city's population:

For each city $i$, let $y_i$ be drunk driving incidents per capita, let $x_i$ be liquor stores per capita, and let $n_i$ be the city's population.

Regression (1) is: $$y_i = a + b x_i + \epsilon_i $$ If you run regression (2) without a constant, you've essentially scaled each observation of regression (1) by the population, that is, you're running:

$$ n_i y_i = a n_i + b n_i x_i + u_i $$

This is weighted least squares, and the weights you're applying are the square of the city's population. That's a lot of weight you're giving the largest cities?!

Note that if you had an observation for each individual in a city and assigned each individual the average value for the city, that would be equivalent to running a regression where you are weighting each city by population (not population squared).

Answer 3

I ran a few experiments on simulated data to see which method works best. Please read my findings below.

Lets look at two different scenarios - First where there is no direct relationship between DUI & Liquor stores & Second where we do have a direct relationship. Then examine each of the methods to see which method works best.

Case 1: No Direct relationship but both are related to the population

library(rmutil)
############
## Simulating Data

set.seed(111)  
# Simulating city populations 
popln <- rpareto(n=10000,m=10000,s=1.2)

# Simulating DUI numbers
e1 <- rnorm(10000,mean=0,sd=15)
DUI = 100 + popln * 0.04 + e1
summary(DUI)
truehist(log(DUI))

# Simulating Nbr of Liquor stores
e2 <- rnorm(100,mean=0,sd=5)
Nbr_Liquor_Stores = 20 + popln * 0.009 + e2
summary(Nbr_Liquor_Stores)
truehist(log(Nbr_Liquor_Stores))

dat <- data.frame(popln,DUI,Nbr_Liquor_Stores)

Now that the data is simulated, lets see how each of the methods fare.

## Method 0: Simple OLS
fit0 <- lm(DUI~Nbr_Liquor_Stores,data=dat)
summary(fit0)

Coefficients:
                   Estimate Std. Error  t value Pr(>|t|)    
(Intercept)       9.4353630  0.2801544    33.68   <2e-16 ***
Nbr_Liquor_Stores 4.4444207  0.0001609 27617.49   <2e-16 ***

Nbr_Liquor_Stores highly significant, as expected. Although the relationship is indirect.

## Method 1: Divide Liquor Stores by population and then regress
fit1 <- lm( I(DUI/popln) ~ Nbr_Liquor_Stores, data=dat)
summary(fit1)

                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        4.981e-01  4.143e-02  12.022   <2e-16 ***
Nbr_Liquor_Stores -1.325e-05  2.380e-05  -0.557    0.578    

Nbr_Liquor_Stores has no significance. Seems to work, but lets not jump to conclusions yet.

## Method 2: Divide Liquor Stores by population and then regress
fit2 <- lm( DUI ~ Nbr_Liquor_Stores + popln, data=dat)
summary(fit2)

                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        1.003e+02  6.022e-01 166.569   <2e-16 ***
Nbr_Liquor_Stores -1.603e-02  3.042e-02  -0.527    0.598    
popln              4.014e-02  2.738e-04 146.618   <2e-16 ***

Nbr_Liquor_Stores not significant, p-value is also quite close to Method 1.

## Method 3: "DUI per capita" on "liquer stores per capita" and "population size" 
fit3 <- lm( I(DUI/popln) ~ I(Nbr_Liquor_Stores/popln) + popln, data=dat)
summary(fit3)

                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 2.841e-02  1.300e-02   2.187   0.0288 *  
I(Nbr_Liquor_Stores/popln)  4.886e+00  1.603e-02 304.867   <2e-16 ***
popln                      -8.426e-09  6.675e-08  -0.126   0.8996    

(Nbr_Liquor_Stores/popln) highly significant! Didn't expect that, maybe this method isn't the best for your problem statement.

Case 2: Direct relationship with both Population & Nbr_Liquor_Stores

### Simulating Data    

set.seed(111)  
# Simulating city populations 
popln <- rpareto(n=10000,m=10000,s=1.2)

# Simulating Nbr of Liquor stores
e2 <- rnorm(100,mean=0,sd=5)
Nbr_Liquor_Stores = 20 + popln * 0.009 + e2
summary(Nbr_Liquor_Stores)
truehist(log(Nbr_Liquor_Stores))

# Simulating DUI numbers
e1 <- rnorm(10000,mean=0,sd=15)
DUI = 100 + popln * 0.021 + Nbr_Liquor_Stores * 0.01 + e1
summary(DUI)
truehist(log(DUI))

dat <- data.frame(popln,DUI,Nbr_Liquor_Stores)

Let's see the performance of each of the methods in this scenario.

## Method 0: Simple OLS
fit0 <- lm(DUI~Nbr_Liquor_Stores,data=dat)
summary(fit0)

                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       5.244e+01  1.951e-01   268.8   <2e-16 ***
Nbr_Liquor_Stores 2.343e+00  1.121e-04 20908.9   <2e-16 ***

Expected, but not a great method to make causal inferences.

## Method 1: Divide Liquor Stores by population and then regress
fit1 <- lm( I(DUI/popln) ~ Nbr_Liquor_Stores, data=dat)
summary(fit1)

                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        4.705e-01  4.005e-02  11.747   <2e-16 ***
Nbr_Liquor_Stores -1.294e-05  2.301e-05  -0.562    0.574    

That is a surprise for me, I was expecting this method to capture the relationship but it doesn't pick it up. So this method fails in this scenario!

## Method 2: Divide Liquor Stores by population and then regress
fit2 <- lm( DUI ~ Nbr_Liquor_Stores + popln, data=dat)
summary(fit2)

                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        1.013e+02  5.945e-01 170.391   <2e-16 ***
Nbr_Liquor_Stores -5.484e-02  2.825e-02  -1.941   0.0523 .  
popln              2.158e-02  2.543e-04  84.875   <2e-16 ***

Nbr_Liquor_Stores is significant, p-value make a lot of sense. A clear winner for me.

## Method 3: "DUI per capita" on "liquer stores per capita" and "population size" 
fit3 <- lm( I(DUI/popln) ~ I(Nbr_Liquor_Stores/popln) + popln, data=dat)
summary(fit3)

                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 6.540e-02  1.485e-02   4.405 1.07e-05 ***
I(Nbr_Liquor_Stores/popln)  3.915e+00  1.553e-02 252.063  < 2e-16 ***
popln                      -2.056e-08  7.635e-08  -0.269    0.788    

TLDR; Method 2 produces most accurate p-values across different scenarios.