# What are R-structure G-structure in a glmm?

I've been using the MCMCglmm package recently. I am confused by what is referred to in the documentation as R-structure and G-structure. These seem to relate to the random effects - in particular specifying the parameters for the prior distribution on them, but the discussion in the documentation seems to assume that the reader knows what these terms are. For example:

optional list of prior specifications having 3 possible elements: R (R-structure) G (G-structure) and B (fixed effects)............ The priors for the variance structures (R and G) are lists with the expected (co)variances (V) and degree of belief parameter (nu) for the inverse-Wishart

...taken from from here.

EDIT: Please note that I have re-written the rest of the question following the comments from Stephane.

Can anyone shed light on what R-structure and G-structure are, in the context of a simple variance components model where the linear predictor is $$\beta_0 + e_{0ij} + u_{0j}$$ with $e_{0ij} \sim N(0,\sigma_{0e}^2)$ and $u_{0j} \sim N(0,\sigma_{0u}^2)$

I made the following example with some data that comes with MCMCglmm

> require(MCMCglmm)
> require(lme4)
> data(PlodiaRB)
> prior1 = list(R = list(V = 1, fix=1), G = list(G1 = list(V = 1, nu = 0.002)))
> m1 <- MCMCglmm(Pupated ~1, random = ~FSfamily, family = "categorical",
+ data = PlodiaRB, prior = prior1, verbose = FALSE)
> summary(m1)

G-structure:  ~FSfamily

post.mean l-95% CI u-95% CI eff.samp
FSfamily    0.8529   0.2951    1.455      160

R-structure:  ~units

post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

Location effects: Pupated ~ 1

post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)   -1.1630  -1.4558  -0.8119    463.1 <0.001 ***
---

> prior2 = list(R = list(V = 1, nu = 0), G = list(G1 = list(V = 1, nu = 0.002)))
> m2 <- MCMCglmm(Pupated ~1, random = ~FSfamily, family = "categorical",
+ data = PlodiaRB, prior = prior2, verbose = FALSE)
> summary(m2)

G-structure:  ~FSfamily

post.mean l-95% CI u-95% CI eff.samp
FSfamily    0.8325   0.3101    1.438    79.25

R-structure:  ~units

post.mean l-95% CI u-95% CI eff.samp
units    0.7212  0.04808    2.427    3.125

Location effects: Pupated ~ 1

post.mean l-95% CI u-95% CI eff.samp  pMCMC
(Intercept)   -1.1042  -1.5191  -0.7078    20.99 <0.001 ***
---

> m2 <- glmer(Pupated ~ 1+ (1|FSfamily), family="binomial",data=PlodiaRB)
> summary(m2)
Generalized linear mixed model fit by the Laplace approximation
Formula: Pupated ~ 1 + (1 | FSfamily)
Data: PlodiaRB
AIC  BIC logLik deviance
1020 1029   -508     1016
Random effects:
Groups   Name        Variance Std.Dev.
FSfamily (Intercept) 0.56023  0.74849
Number of obs: 874, groups: FSfamily, 49

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -0.9861     0.1344  -7.336  2.2e-13 ***


So based on the comments from Stephane I think the G structure is for $\sigma_{0u}^2$. But the comments also say that the R structure is for $\sigma_{0e}^2$ yet this does not seem to appear in the lme4 output.

Note that the results from lme4/glmer() are consistent with both examples from MCMC MCMCglmm.

So, is the R structure for $\sigma_{0e}^2$ and why doesn't this appear in the output for lme4/glmer() ?

-
With the SAS terminology (but it is possibly a more common terminology), the G matrix is the variance matrix of the random effects and the R matrix is the variance matrix of the "errors terms" (in your case perhaps it is the estimated residual variance $\sigma_{0e}^2$ ?) –  Stéphane Laurent Jul 25 '12 at 11:59
@StéphaneLaurent thank you. I wondered if it might be estimated $\sigma_{0e}^2$ but when I first learned about the generalised linear model I remember that $\sigma_{0e}^2$ isn't estimated - only "deviance" is calculated (as with lme4). Maybe I am missing something ? –  Joe King Jul 25 '12 at 12:06
I had never deeply learnt about the theory of glmmm but I don't see why the residual variance wouldn't be estimated. By the way lmer() returns the estimated residual variance. –  Stéphane Laurent Jul 25 '12 at 12:25
@Stéphane Laurent, hmmm, glmer() doesn't return estimated residual variance for some reason (only deviance, loglik, AIC, BIC). Now I am confusing all the things I learned about glmms a few weeks ago :( –  Joe King Jul 25 '12 at 12:39
maybe the sense of the residual variance is not clear when the distribution family is not the Gaussian one –  Stéphane Laurent Jul 25 '12 at 12:43

I would prefer to post my comments below as a comment but this would not be comfortable enough. These are questions rather than an answer (simlarly to @gung I don't feel strong enough on the topic).

I am under the impression that MCMCglmm does not implement a "true" Bayesian glmmm. The true Bayesian model is described in section 2 of this paper. Similarly to the frequentist model, one has $g(E(y \mid u)) = X\beta + Zu$ and there is a prior required on the dispersion parameter $\phi_1$ in addition to the fixed parameters $\beta$ and the "G" variance of the random effect $u$.

But according to this MCMCglmm vignette, the model implemented in MCMCglmm is given by $g(E(y \mid u,e)) = X\beta + Zu + e$ , and it does not involve the dispersion parameter $\phi_1$. It is not similar to the classical frequentist model.

Therefore I would be not surprised that there is no analogue of $\sigma_e$ with glmer.

-
Thank you. Is this topic supposed to be difficult, because I am finding it quite hard ? I think I am satisfied with the meaning of R and G structure now. I am still confused about the lack of $\sigma_e$ with glmer and I am very curious about your comment that MCMCglmm is not truly Bayesian. I can't honestly say I understand all of the paper that you linked to and I am also struggling with parts of the MCMCglmm vignette, but just purely from the perspective of my example, I believe the dispersion parameter $\phi_1$should be constant (because the example is binomial). What am I missing ? –  Joe King Jul 26 '12 at 20:30
Sorry, my words were not totally appropriate. MCMCglmm is truly Bayesian, but it does not exactly implement the classical glmm (I think). In addition you have to be aware that it is difficult to set priors yielding an inference on the variance components close to the frequentist inference. –  Stéphane Laurent Jul 27 '12 at 4:59
Thanks again. In my studying I have found that I can use the default inverse-wishart distribution for variance components in MCMCglmm using a variety of parameters, and the 95% credible intervals always contain the variance value for the random effects estimate by glmer so I felt that this was reasonable, but how should I interpret this case, which may not be typical, where the result that the MCMCglmm intervals are not very sensitive to choice of prior ? Maybe I should ask a new question about this ? –  Joe King Jul 27 '12 at 6:02
Maybe you have a big sample size ? In regards to your initial question, I am under the impression that, at least for the binomial case, the glmer model is equivalent to the MCMCglmm model with $\sigma_e=0$. What happens if you set a prior on $\sigma_e$ highly concentrated at $0$ ? –  Stéphane Laurent Jul 27 '12 at 6:53
Yes, I have a quite large sample size: 50,000 observations in 225 clusters (my own data, not the example in my question). When I set a prior very concentrated near zero on $\sigma_e$, by setting V=0.01 and nu=100 then I obtain 0.25 (CI: 0.16, 0.29) for $\sigma_e$ and 0.53 (0.38, 0.73) for $\sigma_u$. When I set a less informative prior, with V=10 and nu=0.01 then I obtain 0.18 (0.12, 0.23) and 0.49 (0.34, 0.63) respectively. This compares with 0.51 from glmer. I even tried an improper flat prior, which gave 0.10 (0.08, 0.13) and 0.47 (0.25, 0.68). –  Joe King Jul 27 '12 at 20:10

I am late to the game, but a few notes. The $\mathbf{R}$ structure is the residual structure. In your case, the "structure" only has a single element (but this need not be the case). For Gaussian response variable, the residual variance, $\sigma^{2}_{e}$ is typically estimated. For binary outcomes, it is held constant. Because of how MCMCglmm is setup, you cannot fix it at zero, but it is relatively standard to fix it at $1$ (also true for a probit model). For count data (e.g., with a poisson distribution), you do not fix it and this automatically estimates an overdispersion parameter essentially.

The $\mathbf{G}$ structure is the random effects structure. Again in your case, just a random intercept, but if you had multiple random effects, they would form a variance-covariance matrix, $\mathbf{G}$.

A final note, because the residual variance is not fixed at zero, the estimates will not match those from glmer. You need to rescale them. Here is a little example (not using random effects, but it generalizes). Note how the R structure variance is fixed at 1.

# example showing how close the match is to ML without separation
m2 <- MCMCglmm(vs ~ mpg, data = mtcars, family = "categorical",
prior = list(
B = list(mu = c(0, 0), V = diag(2) * 1e10),
R = list(V = 1, fix = 1)),
nitt = 1e6, thin = 500, burnin = 10000)
summary(m2)


Here is the rescaling constant for the binomial family:

k <- ((16*sqrt(3))/(15*pi))^2


Now divide the solution by it, and get the posterior modes

posterior.mode(m2\$Sol/(sqrt(1 + k)))


Which should be fairly close to what we get from glm

summary(glm(vs ~mpg, data = mtcars, family = binomial))

-