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()
this may simply because birds have decreased in numbers (rather than foraging behavior having changed)
I may be missing something, but wouldn't it be easier to include the number of birds as a covariate in the model for bird attacks? You can even add interactions if you believe the number of birds affects the relationship the other covariates have on bird attacks. $\endgroup$