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

[Jarrod Hadfield]

Hi Fiona,

If you set up the measurement error model explicitly (without using mev � see code below) the predict function should work. Note that the default is to have the random effects marginalised and so if you want them in the prediction interval you will need marginal = NULL or  marginal=~idh(sqrt(sig.obs)):units or marginal=~ BirdID dpending on exactly what you want.

Cheers,

Jarrod

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.
                 G2 = list(V=1, fix=1) #mev random effects
               ))

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


From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on behalf of Fiona Scarff <fiona.scarff.4 using gmail.com>
Date: Monday, 18 September 2023 at 16:29
To: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org>
Subject: [R-sig-ME] How to get prediction interval from R package MCMCglmm fitwith measurement error *and* a random effect?
This email was sent to you by someone outside the University.
You should only click on links or attachments if you are certain that the email is genuine and the content is safe.

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]<http://cran.nexr.com/web/packages/MCMCglmm/vignettes/CourseNotes.pdf%5d>
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

_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. Is e buidheann carthannais a th� ann an Oilthigh Dh�n �ideann, cl�raichte an Alba, �ireamh cl�raidh SC005336.

	[[alternative HTML version deleted]]