Prediction intervals given by mod_results() don't work with random effects specified with an '~ inner | outer' term

[befriendabacterium]

Describe the bug When specifying a random effect with an '~ inner | outer' term, the prediction intervals outputted are identical to the confidence intervals (which are correct). You can see the issue by comparing to the predict() function output for models without and with '~ inner | outer' terms (see code below).

To Reproduce

library(metadat)

# PREPARE DATA (ISHAK 2007) -----------------------------------------------

### copy data into 'dat' and examine data
dat <- dat.ishak2007
head(dat, 5)
## Not run:
### load metafor package
library(metafor)
### create long format dataset
dat <- reshape(dat, direction="long", idvar="study", v.names=c("yi","vi"),
               varying=list(c(2,4,6,8), c(3,5,7,9)))
dat <- dat[order(study, time),]
### remove missing measurement occasions from dat.long
dat <- dat[!is.na(yi),]
rownames(dat) <- NULL
head(dat, 8)

# PLOT DATA ---------------------------------------------------------------

### plot data
with(dat, interaction.plot(time, study, yi, type="b", pch=19, lty="solid", xaxt="n",
                           legend=FALSE, xlab="Time Point", ylab="Mean Difference", bty="l"))
axis(side=1, at=1:4, lab=c("1 (3 months)", "2 (6 months)", "3 (12 months)", "4 (12+ months)"))

# MODEL 1: SIMPLE RANDOM EFFECT -------------------------------------------

### intercept-only model with study as a random effect
res1 <- rma.mv(yi, vi, random =  list(~1 | study),
              data = dat)
print(res1, digits=2)

#expected result
predict(res1)

 pred     se    ci.lb    ci.ub    pi.lb    pi.ub 
 -27.1854 0.9450 -29.0375 -25.3333 -37.8370 -16.5339 

#actual (orchaRd) result - identical to above expected, and prediction interval is different from confidence interval
mod_results1<-orchaRd::mod_results(model = res1, mod = "1", group='study')
mod_results1

name  estimate   lowerCL   upperCL   lowerPR   upperPR
1 Intrcpt -27.18543 -29.03754 -25.33333 -37.83702 -16.53384

# MODEL 2: RANDOM EFFECT WITH INNER/OUTER TERM-------------------------------------------

### intercept-only model with heteroscedastic AR(1) structure for the true (random) effects
res2 <- rma.mv(yi, vi, random =  list(~ time | study),
              struct = "CAR", data = dat)
print(res2, digits=2)

#expected result
predict(res2)

   pred     se    ci.lb    ci.ub    pi.lb    pi.ub 
 -27.1854 0.9450 -29.0375 -25.3333 -37.8370 -16.5338 

#actual (orchaRd) result - different from above expected, and prediction interval is identical to confidence interval
mod_results2<-orchaRd::mod_results(model = res2, mod = "1", group='study')
mod_results2

     name  estimate   lowerCL   upperCL   lowerPR   upperPR
1 Intrcpt -27.18543 -29.03754 -25.33333 -29.03754 -25.33333

Expected behavior

• Prediction interval should not equal confidence interval
• Predictions should be identical to those outputted by predict() - I am less sure about this as I'm not sure if the underlying statistical methods between predict() and emmeans() (used in mod_results) differ.

Desktop (please complete the following information):

• OS: Windows
• Browser: Firefox
• Version: 115.0.2

[daniel1noble]

Thanks @befriendabacterium. This is the same or similar bug to #46. Will look into this ASAP for these models

[daniel1noble]

@befriendabacterium I think I found the bug here actually. It was rather simple I think. It was a missing tau2. Can you have a look and see if this fixes things for you?

[daniel1noble]

I did a few more checks and I think this is fixed so I'm closing unless we see it pop up.

[befriendabacterium]

Great! Just double checked this on the example I provided and my own data (before and after updating orchaRd), and works a treat! Nice one, thanks.