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

[Fiona Scarff]

Thanks heaps Jarrod, that worked a treat.

Fiona

On Tue, Sep 19, 2023 at 7:30 PM Jarrod Hadfield <j.hadfield using ed.ac.uk> wrote:
>
> 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]
> 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.