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
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
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)  -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
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)  -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
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()