I built up model using an interaction term using the code:
be_model1 <-glmer(number_of_fledged_hatchlings ~ species*nestbox_quality + (1|habitat) + (1|year) , family = poisson, data = be)
The fixed effects "species" comprises two levels: Great tit (GT) and Blue tit (BT) and the fixed effect "nestbox_quality" as well, namely presumably "high" and presumably "low"
Now I have trouble calculating the fittet mean and its standard errors for the interaction term.
E. g. I want to know the fitted average number of Great tit hatchlings within low quality nest boxes and the standard errors of this fitted mean.
So far I tried the following codes:
> allEffects(be_model1, se=TRUE)
model: number_of_fledged_hatchlings ~ species * nestbox_quality
species* effect
nestbox_quality
species high low
GT 6.275469 6.869725
BT 6.547347 5.819534
which didn't give me the standard errors, as well as this code:
library(ggeffects)
> (ggpredict(be_model1, c("nestbox_quality", "species")))
# Predicted counts of number_of_fledged_hatchlings
# species = BT
nestbox_quality | Predicted | 95% CI
-------------------------------------
high | 6.28 | [5.50, 7.16]
low | 6.87 | [6.07, 7.77]
# species = GT
nestbox_quality | Predicted | 95% CI
-------------------------------------
high | 6.55 | [5.80, 7.40]
low | 5.82 | [5.21, 6.50]
Adjusted for:
* habitat = 0 (population-level)
* year = 0 (population-level)
for which I failed to get standard errors as well. And I've tried
> pred_model2<-predict(be_model1,type="response")
>
> tapply(pred_model2,list(species,nestbox_quality), mean)
high low
BT 6.457035 6.838857
GT 6.498190 5.792777
> tapply(pred_model2,list(species,nestbox_quality), std.error)
high low
BT 0.10072955 0.10153544
GT 0.08011282 0.05989635
Here I did recieve standard errors, but I am a little confused about the predicted means as they differ from the values I obtained using the previous lines of code.
You can extract this information from the effects or ggpredict package output (see below), but emmeans may be the easiest way forward.
form <- fledged ~ species*nestbox_quality + (1|habitat) + (1|year)
dd <- expand.grid(species = factor(c("GT", "BT")),
nestbox_quality = factor(c("low", "high"),
levels = c("low", "high")),
year = 2002:2005,
habitat = factor(1:10))
library(lme4)
set.seed(101)
dd$fledged <- simulate(form[-2], newdata = dd,
family = poisson,
newparams = list(beta = rep(1,4),
theta = rep(1,2)))[[1]]
m <- glmer(form, data = dd, family = poisson)
Note that emmeans gives predictions on the link (log) scale by default: use type = "response" to get count-scale predictions
library(emmeans)
emmeans(m, ~species*nestbox_quality)
species nestbox_quality emmean SE df asymp.LCL asymp.UCL
BT low 1.12 0.615 Inf -0.090 2.32
GT low 2.17 0.613 Inf 0.972 3.37
BT high 2.17 0.613 Inf 0.970 3.37
GT high 4.15 0.612 Inf 2.945 5.34
using as.data.frame() will get you the standard errors ...
as.data.frame(ggpredict(m, terms = c("species","nestbox_quality")))
x predicted std.error conf.low conf.high group
1 BT 3.051631 0.6151573 0.9139219 10.18955 low
2 BT 8.774919 0.6130930 2.6386237 29.18158 high
3 GT 8.786812 0.6130915 2.6422076 29.22104 low
4 GT 63.122520 0.6121429 19.0163563 209.52765 high
Note that ggpredict gives predictions on the count (response) scale, but does not transform the standard errors back to the count scale (arguably a bug, or at least an "infelicity")
as.data.frame() to the rescue again:
as.data.frame(Effect(c("species", "nestbox_quality"), m))
species nestbox_quality fit se lower upper
1 BT low 3.051631 1.877233 0.9139219 10.18955
2 GT low 8.786812 5.387119 2.6422076 29.22104
3 BT high 8.774919 5.379841 2.6386237 29.18158
4 GT high 63.122520 38.640005 19.0163563 209.52765
These SEs are appropriately transformed.