[R-sig-ME] How to get prediction interval from R package MCMCglmm fitwith measurement error *and* a random effect?

[Fiona Scarff]

How to get prediction interval from R package MCMCglmm fit with
measurement error *and* a random effect?

I am having trouble producing a prediction interval from a model fit
by the R package MCMCglmm, when both a random effect and a measurement
error is specified in the fit. Many thanks in advance for your time
and any suggestions.

I can reproduce the prediction interval examples in the package author
(Jarrod Hadfield)’s MCMCglmm course notes [link
http://cran.nexr.com/web/packages/MCMCglmm/vignettes/CourseNotes.pdf]
p.47, but with my own model structure I get this error:

Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'x' in selecting a method for
function 't': invalid 'times' argument

 ‘x’ and ‘t’ aren’t variables I have declared so I imagine they must
be objects called internally by the predict function. In the MCMCglmm
course notes (p.47 again), the author does note that the predict
function for MCMCglmm is currently incomplete and needs further
testing, but should be OK for simpler models. Perhaps my model
structure is an instance of a more elaborate model not yet
accommodated. Is there some other way I can generate a prediction
interval?

For context, I am studying the flying height of birds. We obtain
observations of height from satellite tags attached to birds. The
observations are recorded with large measurement error (often
reporting negative height as if the birds were flying underground).
The variance of the measurement error can be predicted based on the
number and position of available satellites. MCMCglmm has a facility
for users to specify a vector of measurement errors associated with
observations on the response variable. We also have repeat
observations on each tagged bird, so it’s desirable to also
incorporate the flying height preferences of individual birds into the
model.

Here is a reproducible example:

Session info…
MCMCglmm_2.35
R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4.1

#Declare some parameters governing the flying height of the birds, for
purposes of simulation
set.seed(756)
alt1 <- 20 # mean of true flying height, in metres above ground level
sig <- 10 #variance in true flying height
# gamma distribution scale parameter:
beta = sig/alt1
#  gamma distribution shape parameter :
alpha = alt1/beta

 #Additional variation in height amongst individual birds
ind.re <- rnorm(n= 20, mean=0, sd = 3) # ‘ind.re’ denotes individual
random effects

#Create dataframe containing simulation data
#20 tagged birds, with 100 observations each = 2000 observations total
data.sim <- data.frame(BirdID = as.factor(rep(1:20, each=100) ) ) %>%
                        mutate(., sig.obs = rep(5,2000), #std dev of
measurement error, just for this example I'm imposing an identical
error for each observation
            true.alt = #true altitude, drawn from a gamma distribution
                         rgamma(n = 2000, shape = alpha, scale = beta) +
                          ind.re[data.sim$BirdID], #individual
variation amongst birds
             obs.alt = #observed altitude
                             true.alt +
rnorm(n = 2000, mean = 0, sd = sig.obs) # imprecision due to measurement error
                        )

 #declare priors
prior <- list(B = list(mu=0, V=1e10), #for fixed effects
               R = list(V=1, nu=0.002),  #for residuals
               G = list( #random effects
                 G1 = list(V=1, nu=0.002) #first random effect, i.e. BirdID.
               ))

#fit the glmm
m <- MCMCglmm::MCMCglmm(
            fixed = obs.alt ~ 1,
            family = "gaussian",
            random = ~BirdID,
            mev = data.sim$sig.obs, #vector of measurement errors
            prior = prior,
            data = data.sim
)
#get a prediction interval
predict(m,  type = "response", interval = "prediction", level = 0.95)[1:3, ]

Note that if the model is refit without BirdID as a random effect,
(i.e. if we comment out ‘random = ~BirdID,’ above), or alternatively
without the measurement error vector (‘mev =  data.sim$sig.obs,’) then
the predict function works fine:

(Output without random effect specified)
       fit      lwr      upr
1 20.33545 5.284893 33.06830
2 20.72960 6.874275 35.33598
3 20.57229 3.561330 33.84181

 Am I doing something wrong here, or is there a workaround? Many thanks again!
Fiona Scarff
Murdoch University