There are various ways of handling this issue. emmeans in particular has some Delta-method machinery that can be used. However, in this case it's not actually very accurate.
Let's dig in a bit further.
m <- fixef(fit_lme4)
s <- c(sqrt(VarCorr(fit_lme4)$ID))
par(las=1, bty = "l")
cc <- curve(plogis(x), from = m-3*s, to = m+3*s, lwd = 2,
ylab = "probability/scaled density")
nn <- dnorm(cc$x, m, s)/dnorm(m, m, s)
polygon(c(cc$x, rev(cc$x)), c(nn, rep(0, length(nn))),
col = adjustcolor("black",alpha = 0.1),
border = NA)
abline(v=m, lty = 2)
abline(v=0.772501, lty=2, col = 2) ## see below
abline(v=mean(dd$outcome), lty=2, col = 4)
The figure shows the inverse-link function (logistic or inverse-logit); the estimated distribution of log-odds across IDs; the mean of the random effects distribution on the link/logit scale (black vertical dashed line); the mean of the distribution on the probability scale (red); and the mean of the observed response values (blue). The last two aren't identical (0.73 vs 0.77), but they're indistinguishable at this scale.
Bias correction via delta method according to emmeans vignette:
emmeans(fit_lme4, ~1, type = "response", bias.adjust = TRUE,
sigma = s)
## 1 prob SE df asymp.LCL asymp.UCL
## overall 0.969 0.0488 Inf 0.362 0.999
This is actually not very good because the range of the random effects is so wide that a quadratic approximation to the inverse-link function is bad ...
A more accurate estimate/correction:
logitnorm::momentsLogitnorm(m, s)
## mean var
## 0.7722507 0.1451088
Doing the same thing by brute force/Monte Carlo:
mean(plogis(rnorm(1e6, m, s)))
## [1] 0.7725701
There are various ways of handling this issue. emmeans in particular has some Delta-method machinery that can be used.
There are various ways of handling this issue. emmeans in particular has some Delta-method machinery that can be used. However, in this case it's not actually very accurate.
Let's dig in a bit further.
m <- fixef(fit_lme4)
s <- c(sqrt(VarCorr(fit_lme4)$ID))
par(las=1, bty = "l")
cc <- curve(plogis(x), from = m-3*s, to = m+3*s, lwd = 2,
ylab = "probability/scaled density")
nn <- dnorm(cc$x, m, s)/dnorm(m, m, s)
polygon(c(cc$x, rev(cc$x)), c(nn, rep(0, length(nn))),
col = adjustcolor("black",alpha = 0.1),
border = NA)
abline(v=m, lty = 2)
abline(v=0.772501, lty=2, col = 2) ## see below
abline(v=mean(dd$outcome), lty=2, col = 4)
The figure shows the inverse-link function (logistic or inverse-logit); the estimated distribution of log-odds across IDs; the mean of the random effects distribution on the link/logit scale (black vertical dashed line); the mean of the distribution on the probability scale (red); and the mean of the observed response values (blue). The last two aren't identical (0.73 vs 0.77), but they're indistinguishable at this scale.
Bias correction via delta method according to emmeans vignette:
emmeans(fit_lme4, ~1, type = "response", bias.adjust = TRUE,
sigma = s)
## 1 prob SE df asymp.LCL asymp.UCL
## overall 0.969 0.0488 Inf 0.362 0.999
This is actually not very good because the range of the random effects is so wide that a quadratic approximation to the inverse-link function is bad ...
A more accurate estimate/correction:
logitnorm::momentsLogitnorm(m, s)
## mean var
## 0.7722507 0.1451088
Doing the same thing by brute force/Monte Carlo:
mean(plogis(rnorm(1e6, m, s)))
## [1] 0.7725701