How to compare model coefficients from models with different distribution family and link functions

Question

I am trying to understand if I can compare two models with different family distributions / link functions directly, or if this does not make sense mathematically.

In my example, I am measuring some aspect of acoustic noise DayL50 and trying to understand how it affects bird abundance model1 and bird foraging model2.

If we take a look at bird abundance:

library(glmmTMB)
model1<-glmmTMB(Birds ~ scale(DayL50)*scale(Med)+        
                        scale(Veg)+scale(Elev)+                
                        +scale(jDay)+   
                        (1|SITE),                           
                      data=PC, 
                      family=nbinom1(link="log"),            
                      ziformula=~1,                       
                      offset = log(p*CF_Offset))

> summary(model1)
 Family: nbinom1  ( log )
Formula: Birds ~ scale(DayL50) * scale(Med) + scale(Veg) + scale(Elev) +scale(jDay) + (1 | SITE)
Zero inflation:         ~1
Data: PC
 Offset: log(p * CF_Offset)

     AIC      BIC   logLik deviance df.resid 
  4681.0   4734.1  -2330.5   4661.0     1486 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev.
 SITE   (Intercept) 0.2634   0.5132  
Number of obs: 1496, groups:  SITE, 20

Overdispersion parameter for nbinom1 family (): 1.42 

Conditional model:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)               0.66252    0.13347   4.964 6.91e-07 ***
scale(DayL50)            -0.24456    0.06052  -4.041 5.32e-05 ***
scale(Med)                0.02232    0.05323   0.419  0.67500    
scale(Veg)                0.03547    0.05813   0.610  0.54175    
scale(Elev)              -0.22777    0.11302  -2.015  0.04387 *  
scale(jDay)               0.11326    0.03650   3.103  0.00192 ** 
scale(DayL50):scale(Med)  0.04591    0.04428   1.037  0.29979    

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.8072     0.3117  -5.798 6.72e-09 ***

There is a negative effect of DayL50 on Birds (which is simply integer counts of birds). The coefficient from this model is -0.24456, yet I understand that this is scaled. So to get the un-scaled version, I divide this by the standard deviation of the raw data to get:

> -0.24456/sd(PC$DayL50)
[1] -0.02165274

So this coefficient should be now relative to units of DayL50, rather than in sd. The goal is to make this coefficient comparable to a coefficient from the same predictor of the next model, which is a measure of bird foraging.

I put fake clay caterpillars in trees, and then scored them to see how many were predated by birds (bird bills are pretty easy to spot in clay). So the next model is binomial - raw data are 0 or 1 for a caterpillar being not attacked or attacked, respectively, by birds.

library(lme4)
model2<-glmer(Battack~scale(DayL50)*scale(Med)+scale(Elev)+scale(Veg)+scale(jDay)+
                (1|SITE), data=C, family=binomial(link="logit"))

> summary(model2)
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: Battack ~ scale(DayL50) * scale(Med) + scale(Elev) + scale(Veg) +      scale(jDay) + (1 | SITE)
   Data: C

     AIC      BIC   logLik deviance df.resid 
   566.6    603.2   -275.3    550.6      712 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-0.8293 -0.4017 -0.3356 -0.2822  5.0913 

Random effects:
 Groups Name        Variance Std.Dev.
 SITE   (Intercept) 0.1475   0.384   
Number of obs: 720, groups:  SITE, 20

Fixed effects:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)              -1.958300   0.158586 -12.348  < 2e-16 ***
scale(DayL50)            -0.465611   0.150432  -3.095  0.00197 ** 
scale(Med)                0.153074   0.157914   0.969  0.33237    
scale(Elev)              -0.084015   0.192976  -0.435  0.66330    
scale(Veg)               -0.006801   0.154334  -0.044  0.96485    
scale(jDay)               0.436336   0.172069   2.536  0.01122 *  
scale(DayL50):scale(Med)  0.195088   0.147389   1.324  0.18563    

There is a negative effect of DayL50 on Battack (which is 0 or 1; not attacked or attacked by birds). The coefficient from this model is -0.465611, but un-scaled is:

> -0.465611/sd(C$DayL50)
[1] -0.04095002

So now I have two unscaled coefficients of the same predictor data DayL50. I would like to compare these two coefficients to understand if they differ.

The reasoning behind this is that model2 may just be another test of model1. For example, if bird abundance decreases with increasing DayL50, and attacks on caterpillars Battack also goes down to a similar degree, this may simply because birds have decreased in numbers (rather than foraging behavior having changed).

Yet I fear that there are potential problems here with comparing these directly.

1) There are two different link functions log and logit, but I am unsure as to whether or not I need to do anything differently, since they are both log transformations?

2) I am using two different packages (glmmTMB and lme4), should I be concerned with differences in how these coefficients would be calculated?

3) Model1 is zero-inflated and includes offsets. Should these things be concerning when comparing coefficients?

If this is this an appropriate comparison, below is a visualization of what I am trying to understand - are these two model coefficients different? It seems to be that this data would suggest that they are not different, since 95% CI overlap substantially.

Here is reproducible code from this plot.

df<-data.frame("model"=c("Bird abundance", "Caterpillar attacks"),"coef"=c(-0.24456,-0.465611),"se"=c(0.06052,0.150432),"sd"=c(11.29464,11.37023))

library(ggplot2)
ggplot(df, aes(x=model, y=(coef/sd))) + 
  geom_errorbar(aes(ymin=((coef-1.96*se)/sd), ymax=((coef+1.96*se)/sd), width=.2)) +
  geom_point()+
  labs(y="Partial slopes of sound pressure level (dB)\n", x="Model")+
  geom_hline(yintercept=0, linetype="dashed")+
  theme_classic()

Answer 1

You can't directly compare the estimated coefficients since the units of the response variable are not the same in both models.

See, a logistic regression will estimate a binomial probability of observing the event you modelled. So a number between 0 and one is ultimately estimated. Note also that the estimate is not linearly related to the covariates you estimate in the model. So the effect depends on the initial value.

Now, in a Negative binomial regression, the outcome is a count, a number between zero and infinity, and again, the covariates are not linearly related to the outcome, due to the link function.

But it doesn't mean you cannot compare both models. It just takes a little more effort. I would plot the partial effects of the covariate of interest (and relevant interactions) to each response variables and build my rationale from there.

Update

My idea is to help you investigate the variables effects onto the different response variables by exposing the model's mechanics. For this we will need an example. In the following code I generate the independent variables x and z. Then I generate the linear dependent portion of a logistic regression (log odds) response_binary and the negative binomial counts and theta using this SO answer.

suppressMessages(library(tidyverse))
suppressMessages(library(Hmisc))
suppressMessages(library(glmmTMB))
suppressMessages(library(broom))
suppressMessages(library(MASS))
suppressMessages(library(modelr))

N <-  500
set.seed(1)
df <- tibble(x = runif(n = N), z = runif(n = N) < 0.3)

bin_link <- function(x) 1/(1 + exp(-x))

df <- 
df %>% 
    mutate(logodds = 1.2*x + 0.2*z - 1.2*x*z + rnorm(N, 0, 0.1) , 
                 mu = 1.2*x + 2*x^2 + z + rnorm(N, 0, 0.1)
                 ) %>% 
    mutate( prob = bin_link(logodds),
                    neg_par = exp(mu),
                    response_binary = rbinom(n=N, size=1, prob=prob),
                    theta = sample(c(5,8,10, 15), replace = T, size=N),
                    response_count = rnbinom(n = N, size = theta, mu=neg_par)
                    )

The resulting distribution of response_count is in the following histogram.

df %>% ggplot(aes(response_count)) + geom_histogram(color='black', fill='skyblue') + labs(title='Count variable distribution')
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

And the scatter plots between x and response_count with z colours indicate an idea of what the model should be capable of identifying.

df %>% ggplot(aes(x = x, y = response_count, color=z)) + geom_point() + labs(title='relation between x, z and the count response')

And to explore the relation between x, z and the response_binary I turn to a lowess plot visualization taken from chapter 12 in Frank Harrel's Regression Modeling Strategies.

df %>% ggplot(aes(x = x, y = response_binary, group=z, color=z)) + histSpikeg(response_binary~x*z, lowess=T, data=df) + 
    labs(title='Estimated lowess for the relation between x, z and the\nproportion/probability of the binary response')

Now to explore the relationship between x + z when it comes to the response_count and the response_binary I suggest you inspect the models partial effect plots, since the coefficients can be directly compared.

First we build two simple models. A negative binomial model nb_md for the response_count variable, and a logistic regression logi_md for the response_binary.

nb_md <- glm.nb(response_count ~ x*z, data=df, link="log")

summary(nb_md)
#> 
#> Call:
#> glm.nb(formula = response_count ~ x * z, data = df, link = "log", 
#>     init.theta = 8.292271032)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -3.1667  -0.7803  -0.1371   0.5643   2.7400  
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept) -0.36614    0.08754  -4.183 2.88e-05 ***
#> x            3.37600    0.12373  27.285  < 2e-16 ***
#> zTRUE        0.98352    0.13262   7.416 1.21e-13 ***
#> x:zTRUE     -0.08763    0.19660  -0.446    0.656    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for Negative Binomial(8.2923) family taken to be 1)
#> 
#>     Null deviance: 2448.32  on 499  degrees of freedom
#> Residual deviance:  548.59  on 496  degrees of freedom
#> AIC: 2455.9
#> 
#> Number of Fisher Scoring iterations: 1
#> 
#> 
#>               Theta:  8.29 
#>           Std. Err.:  1.24 
#> 
#>  2 x log-likelihood:  -2445.94


logi_md <- glm(response_binary ~ x + z, data = df, family = 'binomial')

summary(logi_md)
#> 
#> Call:
#> glm(formula = response_binary ~ x + z, family = "binomial", data = df)
#> 
#> Deviance Residuals: 
#>    Min      1Q  Median      3Q     Max  
#> -1.639  -1.248   0.829   1.031   1.451  
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.03917    0.19276   0.203 0.838968    
#> x            1.00490    0.33105   3.035 0.002401 ** 
#> zTRUE       -0.67308    0.20188  -3.334 0.000856 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 679.64  on 499  degrees of freedom
#> Residual deviance: 659.00  on 497  degrees of freedom
#> AIC: 665
#> 
#> Number of Fisher Scoring iterations: 4

The models could be better built, in order to better capture the non-linearity we can see in the exploration plots.

However I build the partial effects plot and plot the predictions together with the standard error. From these we can have an idea of how the variable X can have a different effect over response_count and response_binary. Here are some observations we can make:

• the effect of variable x greatly increases for larger values of x when it comes to understanding it's effect on response_count, being greater when z == T

• the effect of variable x on response_binary tends to be the same, regardless of z in lower ranges. But it tends to be much larger with z == F if x is large.

  partial_effects_frame <- 
  tibble(x = seq(0,1, 0.01)) %>% 
      tidyr::crossing(z=c(T,F))
  
  
  bind_rows(
          augment(nb_md, newdata= partial_effects_frame, type.predict='response') %>% mutate(model='nb_md'),
          augment(logi_md, newdata= partial_effects_frame, type.predict='response') %>% mutate(model='logi_md')
  ) %>%  
      ggplot(aes(x, .fitted, group=z)) + 
      geom_ribbon(aes(ymin= .fitted - .se.fit, ymax= .fitted + .se.fit), alpha=0.2) +
      geom_line(aes(group=z, color=z)) + 
      facet_wrap(~model, scales = 'free', ncol=1) + 
      labs(title = 'Partial effects of the variable x interacted with z')
  

Now this is just a simple example of what I meant with a little more effort. Changing base values, exploring the plots and everything, working with tables might help you study the effect on different response variables.

^{Created on 2019-11-22 by the reprex package (v0.3.0)}