5
$\begingroup$

I am using a hurdle model (binomial logit for the zero counts, truncated Poisson for the positive counts) to study the effects of demographic (and other) variables on the adoption of private air-quality sensors. I want to model adoption rates, therefore I use population as an offset.

My unit of observation is the municipality. I am using the following demographics:

  • average income
  • average age
  • population density
  • average years of education
  • percentage of females (<-- my question is about this variable)

From a first simple comparison I see that municipalities with private air-quality sensors have a significantly higher percentage of females. However, the results of my regression show the complete opposite and with a weird magnitude:

Call:
hurdle(formula = private_sensors ~ lag_private_sensors + 
                 gov_monitors + density + avg_income + avg_age + 
                 fem_percent + years_edu + offset(log(k_people)), 
       data = df_2019_ita, 
       dist = "poisson", zero.dist = "binomial", link = "logit")
    
    Pearson residuals:
         Min       1Q   Median       3Q      Max 
    -1.87636 -0.14096 -0.08435 -0.05103 33.91909 
    
    Count model coefficients (truncated poisson with log link):
                         Estimate Std. Error z value Pr(>|z|)    
    (Intercept)         -17.45150    5.89024  -2.963 0.003049 ** 
    lag_private_sensors   0.02504    0.06401   0.391 0.695626    
    gov_monitors         -0.21234    0.02574  -8.249  < 2e-16 ***
    density              -0.43161    0.05392  -8.005 1.19e-15 ***
    avg_income            0.21397    0.04393   4.871 1.11e-06 ***
    avg_age               0.23298    0.06200   3.757 0.000172 ***
    fem_percent         -13.92275   14.53042  -0.958 0.337972    
    years_edu             0.61405    0.19028   3.227 0.001250 ** 
    Zero hurdle model coefficients (binomial with logit link):
                          Estimate Std. Error z value Pr(>|z|)    
    (Intercept)          -0.699734   3.866614  -0.181   0.8564    
    lag_private_sensors   0.200362   0.096645   2.073   0.0382 *  
    gov_monitors          0.161704   0.075342   2.146   0.0319 *  
    density              -0.223444   0.098359  -2.272   0.0231 *  
    avg_income            0.122489   0.025001   4.899 9.62e-07 ***
    avg_age               0.001647   0.036823   0.045   0.9643    
    fem_percent         -18.177834   8.444732  -2.153   0.0314 *  
    years_edu             0.197443   0.185116   1.067   0.2862    
    ---
    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
    
    Number of iterations in BFGS optimization: 53 
    Log-likelihood: -930.6 on 16 Df

As you can see, fem_percent has a big negative coefficient, that means that the exponentiated coefficient will be 0.0000009 for the count part and 0.000000012 for the zero part. Check all the exponentiated coefficients:

A tibble: 16 × 2
   Coeff                            Value
   <chr>                            <dbl>
 1 count_(Intercept)         0.0000000264
 2 count_lag_private_sensors 1.03        
 3 count_gov_monitors        0.809       
 4 count_density             0.649       
 5 count_avg_income          1.24        
 6 count_avg_age             1.26        
 7 count_fem_percent         0.000000898 
 8 count_years_edu           1.85        
 9 zero_(Intercept)          0.497       
10 zero_lag_private_sensors  1.22        
11 zero_gov_monitors         1.18        
12 zero_density              0.800       
13 zero_avg_income           1.13        
14 zero_avg_age              1.00        
15 zero_fem_percent          0.0000000127
16 zero_years_edu            1.22        

So, for example, increasing the average income in a municipality by one unit (1000€) is associated with an increase in the odds of having a positive adoption rate (sensors/population) by 1.13 times or 13%.

Question: How am I to interpret the coefficient on the percentage of females? There clearly is something wrong, and I also can't explain to myself why wouldn't the coefficient be positive. Does it have to do with the offset? Any help would be very appreciated because I am not that experienced with count models and in general (this is the first analysis that I do by myself).

Edit: Additional information

Correlation/Collinearity

Correlation among variables had been checked by plotting a corrplot in R and by running performance::check_collinearity(model):

enter image description here

performance::check_collinearity(modelHurdle.ita.lag)
# Check for Multicollinearity

* conditional component:

Low Correlation

                Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
 lag_private_sensors 1.05 [1.03, 1.08]         1.02      0.95     [0.93, 0.97]
        gov_monitors 1.51 [1.47, 1.56]         1.23      0.66     [0.64, 0.68]
             density 3.02 [2.91, 3.14]         1.74      0.33     [0.32, 0.34]
          avg_income 4.94 [4.75, 5.14]         2.22      0.20     [0.19, 0.21]
             avg_age 1.58 [1.53, 1.62]         1.26      0.63     [0.62, 0.65]
         fem_percent 2.35 [2.27, 2.44]         1.53      0.43     [0.41, 0.44]
           years_edu 3.34 [3.22, 3.47]         1.83      0.30     [0.29, 0.31]

* zero inflated component:

Low Correlation

                Term  VIF   VIF 95% CI Increased SE Tolerance Tolerance 95% CI
 lag_private_sensors 1.06 [1.04, 1.09]         1.03      0.94     [0.92, 0.96]
        gov_monitors 1.27 [1.24, 1.31]         1.13      0.79     [0.77, 0.81]
             density 1.27 [1.23, 1.30]         1.13      0.79     [0.77, 0.81]
          avg_income 1.66 [1.62, 1.72]         1.29      0.60     [0.58, 0.62]
             avg_age 1.17 [1.15, 1.21]         1.08      0.85     [0.83, 0.87]
         fem_percent 1.49 [1.45, 1.54]         1.22      0.67     [0.65, 0.69]
           years_edu 1.99 [1.92, 2.05]         1.41      0.50     [0.49, 0.52]

I add the results from the vif() function as asked by @EdM: I had to split the model into the two components to apply it because I got a warning message (followed inputs from this question: vif() with more than 1 set of coefficients)

(I split the model into two following a procedure from this question: Truncated Poisson vs Hurdle model)

## Hurdle part:
hurdlePart <- glm(formula = I(private_sensors>0) ~ 
    lag_private_sensors + gov_monitors + density + avg_income + 
    avg_age + fem_percent + years_edu + offset(log(k_people)), 
                  data    = df_2019_ita,
                  family  = binomial(link = "logit"))
car::vif(hurdlePart)

# Positive/truncated poisson part:
pos.poiss <- vglm(formula = private_sensors ~ lag_private_sensors + gov_monitors + density + avg_income + avg_age + fem_percent + years_edu + offset(log(k_people)), 
                  data = filter(df_2019_ita, private_sensors > 0),
                  family = pospoisson())
car::vif(pos.poiss)

> car::vif(hurdlePart)
lag_private_sensors        gov_monitors             density          avg_income             avg_age         fem_percent 
                1.1                 1.3                 1.3                 1.7                 1.2                 1.5 
          years_edu 
                2.0 


> car::vif(pos.poiss)
                    GVIF Df GVIF^(1/(2*Df))
lag_private_sensors 14.3  0             Inf
gov_monitors         1.0  1             1.0
density              1.5  1             1.2
avg_income           3.0  1             1.7
avg_age              4.9  1             2.2
fem_percent          1.6  1             1.3
years_edu            2.4  1             1.5

I noticed something interesting/frightening when running the positive part of the hurdle model:

> summary(pos.poiss)

Call:
vglm(formula = private_sensors ~ lag_private_sensors + gov_monitors + 
    density + avg_income + avg_age + fem_percent + years_edu + 
    offset(log(k_people)), family = pospoisson(), 
    data = filter(df_2019_ita, private_sensors > 0))

Coefficients: 
                    Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -17.4510     5.8903      NA       NA    
lag_private_sensors   0.0250     0.0640    0.39  0.69563    
gov_monitors         -0.2123     0.0257   -8.25  < 2e-16 ***
density              -0.4316     0.0539   -8.01  1.2e-15 ***
avg_income            0.2140     0.0439    4.87  1.1e-06 ***
avg_age               0.2330     0.0621    3.75  0.00017 ***
fem_percent         -13.9244    14.5322      NA       NA    
years_edu             0.6141     0.1903    3.23  0.00125 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Name of linear predictor: loglink(lambda) 

Log-likelihood: -276 on 179 degrees of freedom

Number of Fisher scoring iterations: 7 

Warning: Hauck-Donner effect detected in the following estimate(s):
'(Intercept)', 'gov_monitors', 'fem_percent'

I think the NAs are due to something gone wrong in the maximum likelihood, although they do not appear using the hurdle() function. Could someone please confirm or comment/why does it not show up in the hurdle function? However, this problem disappears when I use the rescaled fem_percent (see below for more details).

lag_private_sensors, fem_percent, private_sensors

private_sensors is a count variable that measures the number of sensors present in a municipality in a given year (I use 2019 for the hurdle regression). It has excess number of zeros (see histogram, that's why I'm using the hurdle model)

lag_private_sensors is a measure of the presence of sensors in nearby municipalities in the previous year (2018) and is calculated with the poly2nb() and lag.listw() functions in R.

fem_percent is the percent of females in a given municipality.

enter image description here

enter image description here

enter image description here

In case it's useful, I add the following boxplots that show the difference in characteristics between municipalities with and without sensors:

enter image description here

Rescaled variable female

I tried to rescale the variable on female as follows:

df_2019_ita$rescaled_fem_percent <- scale(df_2019_ita$fem_percent)[,1]

and I got a much more reasonable coefficient! although still negative, but at least I can try to interpret it. How tho? I would maybe say something like: if fem_percent increased by one standard deviation, the odds of having a positive sensors-count would decrease by exp(-0.30201) = 0.74 times (or 26%). Does that make sense?

> summary(modelHurdle.ita.lag.rescaled)

Call:
hurdle(formula = private_sensors ~ lag_private_sensors + 
    gov_monitors + density + avg_income + avg_age + 
    rescaled_fem_percent + 
    years_edu + offset(log(k_people)), data = df_2019_ita, 
    dist = "poisson", zero.dist = "binomial", link = "logit")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.8764 -0.1410 -0.0844 -0.0510 33.9191 

Count model coefficients (truncated poisson with log link):
                     Estimate Std. Error z value Pr(>|z|)    
(Intercept)          -24.4694     3.3084   -7.40  1.4e-13 ***
lag_private_sensors    0.0250     0.0640    0.39  0.69563    
gov_monitors          -0.2123     0.0257   -8.25  < 2e-16 ***
density               -0.4316     0.0539   -8.01  1.2e-15 ***
avg_income             0.2140     0.0439    4.87  1.1e-06 ***
avg_age                0.2330     0.0620    3.76  0.00017 ***
rescaled_fem_percent  -0.2313     0.2414   -0.96  0.33791    
years_edu              0.6141     0.1903    3.23  0.00125 ** 
Zero hurdle model coefficients (binomial with logit link):
                     Estimate Std. Error z value   Pr(>|z|)    
(Intercept)          -9.86207    2.31451   -4.26 0.00002035 ***
lag_private_sensors   0.20036    0.09665    2.07      0.038 *  
gov_monitors          0.16170    0.07534    2.15      0.032 *  
density              -0.22344    0.09836   -2.27      0.023 *  
avg_income            0.12249    0.02500    4.90 0.00000096 ***
avg_age               0.00165    0.03682    0.04      0.964    
rescaled_fem_percent -0.30201    0.14030   -2.15      0.031 *  
years_edu             0.19744    0.18512    1.07      0.286    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Number of iterations in BFGS optimization: 30 
Log-likelihood: -931 on 16 Df
$\endgroup$
7
  • $\begingroup$ Looks like a problem with collinearity. Have you checked that? $\endgroup$
    – Peter Flom
    Commented Jan 28 at 12:36
  • $\begingroup$ Hi Peter, thank you for your comment. I will edit my post to include the checks I did, although I am not sure if they are sufficient. $\endgroup$
    – Giovanna
    Commented Jan 28 at 13:30
  • $\begingroup$ I guess it's not collinearity. $\endgroup$
    – Peter Flom
    Commented Jan 28 at 13:38
  • $\begingroup$ ahah:D Another thing I checked was how to interpret compositional data: "fem_percent" is a percentage and thus the fact that fem_percent increases means that male_percent has to decrease, and interpreting this is apparently not always straightforward. But I don't know if this could be the cause of such a weird coefficient. $\endgroup$
    – Giovanna
    Commented Jan 28 at 13:45
  • 1
    $\begingroup$ Could you please edit the question to show the distributions of fem_percent and private_sensors values and to explain the lag_private_sensors predictor? Also, your correlation plot seems to include the outcome private_sensors , not the predictor lag_private_sensors. That plot suggests a very low two-way correlation between fem_percent and private_sensors. I wonder whether there might be some numerical precision problem leading to the extreme magnitudes of the intercept and fem_percent values. Rescaling predictors (e..g, to fem_fraction) might help. $\endgroup$
    – EdM
    Commented Jan 28 at 15:22

1 Answer 1

7
$\begingroup$

The weird magnitudes of the original regression coefficients evidently came from numerical difficulties in maximum-likelihood estimation.

It often makes sense to center and scale the values of continuous predictors in a model like this, then re-adjust the intercept and coefficients to account for the centering and scaling. For example, the coxph() function in R both centers and scales behind the scenes, to avoid the types of problems that you found.

There are two major reasons.

First, with usual R coding, the intercept is the estimated value (here, in log-odds for the hurdle or log scale for the Poisson part) when continuous predictors have 0 values and categorical predictors are at reference values. You already have a very low probability of having sensors when your continuous predictors are at typical values, most of which are probably far from 0. To estimate the intercept you then must extrapolate the results for those predictors down to 0.

Second, a coefficient estimate is inversely proportional to the measurement scale of the corresponding predictor. In models like this, a predictor with numerically high-magnitude values will thus tend to have a numerically low-magnitude regression coefficient.

Both those issues can push you to the numerical precision limits of the computer. Centering predictors addresses the first issue, scaling the second.

Your interpretation of the rescaled_fem_percent coefficient in the last binomial model, after centering and scaling fem_percent, is correct with one caveat: it's for the situation when all other predictor values are held constant. The negative value of the coefficient in the final model, despite your observation that "municipalities with private air-quality sensors have a significantly higher percentage of females," probably comes from correlations of fem_percent with other predictors (e.g.,avg_income and years_edu)that are themselves associated with having private sensors.

To get the rescaled_fem_percent coefficient expressed in terms of its original units, you divide the coefficient (in log-odds scale) for rescaled_fem_percent by the standard deviation of fem_percent and then exponentiate to get the odds per unit change in the original scale.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.