2

I need to extract the posterior estimates and intervals for a random effect from my model.

For illustrative purposes, a similar dataset to the one I am using would be the ChickWeight dataset in base R.

The way I extract the posterior estimates and intervals for my fixed effects is like so:

#load package
library(lme4)

#model
m.surv<-lmer(weight ~ Time + Diet + (1|Chick), data=ChickWeight)

#load packages
library(MCMCglmm)
library(arm)

#set up for fixed effects
sm.surv<-sim(m.surv)
smfixef.surv=sm.surv@fixef
smfixef.surv=as.mcmc(smfixef.surv)

#which gives
> posterior.mode(smfixef.surv)
(Intercept)        Time       Diet2  ... 
  8.5963329   8.7034260   5.1220436  ...
> HPDinterval(smfixef.surv)
                   lower      upper
(Intercept) -0.90309142 21.3617805
Time         8.42279728  9.0306337
Diet2       -6.84371527 35.1745980
...
attr(,"Probability")
[1] 0.95
>

When I try this for the random effect (Chick), I get the following error at the second line of code:

smranef.surv=sm.surv@ranef
smranef.surv=as.mcmc(smranef.surv)

Error in mcmc.list(x) : Arguments must be mcmc objects

Any suggestions on how I can modify my code to extract these values for the random effect?

Note for other users: if the model would have been a MCMCglmm model, the posterior mode values for the MCMC output for the random effects can be extracted like so:

posterior.mode(sm.surv$VCV[,1])
HPDinterval(sm.surv$VCV[,1])
Blundering Ecologist
  • 1,199
  • 2
  • 14
  • 38
  • `str(sm.surv)` suggests `sm.surv@ranef` and `sm.surv@fixef` are rather different – Henry Nov 24 '19 at 21:02
  • That makes sense. In r, `fixef` is a function to extract fixed effects estimates and `ranef` is a function to extract the random effects estimates. – Blundering Ecologist Nov 24 '19 at 21:56
  • You may wish to utilize the `merTools` package. Alternatively, perhaps an `rstan` implementation – alexwhitworth Nov 27 '19 at 16:38
  • 1
    To expand on Henry's comment, str(smfixef.surv) shows that smfixef.surv is a 100x5 matrix. str(ramranef.surv) shows that smranef.surv is a list containing a 100x50x1 array. I was able to apply as.mcmc() when I extracted the smranef.surv=as.mcmc(smranef.surv[[1]][,,1]). I think this works because ?as.mcmc revealed that the function takes vectors or matrices, not lists or higher dimensional arrays. – gregor-fausto Nov 28 '19 at 04:09
  • @fausto.siegmund That does help it run! However, the output (and maybe this is because random effects are different from fixed effects) includes an estimate and interval for each level within the random effect, not one overall estimate and interval. Any ideas on how to properly consolidate these values so that there is one estimate and one set of intervals for the random effect? – Blundering Ecologist Nov 28 '19 at 19:00

1 Answers1

1

To extract the estimate and 95% CI for your random effects, you use the following code:

sm.surv <-sim(m.surv)

#between Chick variance
bChick <-sm@ranef$Chick[,,1]
bvar<-as.vector(apply(bChick, 1, var)) #between ind variance posterior distribution
bvar<-as.mcmc(bvar)
posterior.mode(bvar) #mode of the distribution
HPDinterval(bvar)

This will then give you:

>     posterior.mode(bvar)
     var1 
     501.24353 
>     HPDinterval(bvar)
      lower   upper
var1 412.36042 630.201
attr(,"Probability")
[1] 0.95

This means that the estimate is 501 and the lower 95% interval was 412 and the upper 95% interval was 630.

Blundering Ecologist
  • 1,199
  • 2
  • 14
  • 38