Quantifying uncertainty around R-squared for generalized linear mixed models

[This article was first published on Ecology in silico, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

People love $R^2$. As such, when Nakagawa and Schielzeth published an article in the journal Methods in Ecology and Evolution earlier this year, ecologists (amid increasing use of generalized linear mixed models (GLMMs)) rejoiced. Now there’s an R function that automates $R^2$ calculations for GLMMs fit with the lme4 package.

$R^2$ is usually reported as a point estimate of the variance explained by a model, using the maximum likelihood estimates of the model parameters and ignoring uncertainty around these estimates. Nakagawa and Schielzeth (2013) noted that it may be desirable to quantify the uncertainty around $R^2$ using MCMC sampling. So, here we are.

Background

$R^2$ quantifies the proportion of observed variance explained by a statistical model. When it is large (near 1), much of the variance in the data is explained by the model.

Nakagawa and Schielzeth (2013) present two $R^2$ statistics for generalized linear mixed models:

1) Marginal $R^2_{GLMM(m)}$, which represents the proportion of variance explained by fixed effects:

where $\sigma^2_f$ represents the variance in the fitted values (on a link scale) based on the fixed effects:

$\boldsymbol{X}$ is the design matrix of the fixed effects, and $\boldsymbol{\beta}$ is the vector of fixed effects estimates.

$\sum_{l=1}^{u}\sigma^2_l$ represents the sum the variance components for all of $u$ random effects. $\sigma^2_d$ is the distribution-specific variance (Nakagawa & Schielzeth 2010), and $\sigma^2_d$ represents added dispersion.

2) Conditional $R^2_{GLMM(c)}$ represents the proportion of variance explained by the fixed and random effects combined:

Point-estimation of $R^2_{GLMM}$

Here, I’ll follow the example of an overdispersed Poisson GLMM provided in the supplement to Nakagawa & Schielzeth, available here. This is their most complicated example, and the simpler ones ought to be relatively straightforward for those that are interested in normal or binomial GLMMs.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
<span class="line"><span class="c1"># First, simulate data (code from Nakagawa & Schielzeth 2013): </span>
</span><span class="line"><span class="c1"># 12 different populations n = 960</span>
</span><span class="line">Population <span class="o"><-</span> gl<span class="p">(</span><span class="m">12</span><span class="p">,</span> <span class="m">80</span><span class="p">,</span> <span class="m">960</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># 120 containers (8 individuals in each container)</span>
</span><span class="line">Container <span class="o"><-</span> gl<span class="p">(</span><span class="m">120</span><span class="p">,</span> <span class="m">8</span><span class="p">,</span> <span class="m">960</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Sex of the individuals. Uni-sex within each container (individuals are</span>
</span><span class="line"><span class="c1"># sorted at the pupa stage)</span>
</span><span class="line">Sex <span class="o"><-</span> factor<span class="p">(</span>rep<span class="p">(</span>rep<span class="p">(</span>c<span class="p">(</span><span class="s">"Female"</span><span class="p">,</span> <span class="s">"Male"</span><span class="p">),</span> each <span class="o">=</span> <span class="m">8</span><span class="p">),</span> <span class="m">60</span><span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Habitat at the collection site: dry or wet soil (four indiviudal from</span>
</span><span class="line"><span class="c1"># each Habitat in each container)</span>
</span><span class="line">Habitat <span class="o"><-</span> factor<span class="p">(</span>rep<span class="p">(</span>rep<span class="p">(</span>c<span class="p">(</span><span class="s">"dry"</span><span class="p">,</span> <span class="s">"wet"</span><span class="p">),</span> each <span class="o">=</span> <span class="m">4</span><span class="p">),</span> <span class="m">120</span><span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Food treatment at the larval stage: special food ('Exp') or standard</span>
</span><span class="line"><span class="c1"># food ('Cont')</span>
</span><span class="line">Treatment <span class="o"><-</span> factor<span class="p">(</span>rep<span class="p">(</span>c<span class="p">(</span><span class="s">"Cont"</span><span class="p">,</span> <span class="s">"Exp"</span><span class="p">),</span> <span class="m">480</span><span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Data combined in a dataframe</span>
</span><span class="line">Data <span class="o"><-</span> data.frame<span class="p">(</span>Population <span class="o">=</span> Population<span class="p">,</span>
</span><span class="line">    Container <span class="o">=</span> Container<span class="p">,</span> Sex <span class="o">=</span> Sex<span class="p">,</span>
</span><span class="line">    Habitat <span class="o">=</span> Habitat<span class="p">,</span> Treatment <span class="o">=</span> Treatment<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Subset the design matrix (only females express colour morphs)</span>
</span><span class="line">DataF <span class="o"><-</span> Data<span class="p">[</span>Data<span class="o">$</span>Sex <span class="o">==</span> <span class="s">"Female"</span><span class="p">,</span> <span class="p">]</span>
</span><span class="line">
</span><span class="line"><span class="c1"># random effects</span>
</span><span class="line">PopulationE <span class="o"><-</span> rnorm<span class="p">(</span><span class="m">12</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> sqrt<span class="p">(</span><span class="m">0.4</span><span class="p">))</span>
</span><span class="line">ContainerE <span class="o"><-</span> rnorm<span class="p">(</span><span class="m">120</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> sqrt<span class="p">(</span><span class="m">0.05</span><span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># generation of response values on link scale (!) based on fixed effects,</span>
</span><span class="line"><span class="c1"># random effects and residual errors</span>
</span><span class="line">EggLink <span class="o"><-</span> with<span class="p">(</span>DataF<span class="p">,</span>
</span><span class="line">                  <span class="m">1.1</span> <span class="o">+</span>
</span><span class="line">                  <span class="m">0.5</span> <span class="o">*</span> <span class="p">(</span>as.numeric<span class="p">(</span>Treatment<span class="p">)</span> <span class="o">-</span> <span class="m">1</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">                  <span class="m">0.1</span> <span class="o">*</span> <span class="p">(</span>as.numeric<span class="p">(</span>Habitat<span class="p">)</span> <span class="o">-</span> <span class="m">1</span><span class="p">)</span> <span class="o">+</span>
</span><span class="line">                  PopulationE<span class="p">[</span>Population<span class="p">]</span> <span class="o">+</span>
</span><span class="line">                  ContainerE<span class="p">[</span>Container<span class="p">]</span> <span class="o">+</span>
</span><span class="line">                  rnorm<span class="p">(</span><span class="m">480</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> sqrt<span class="p">(</span><span class="m">0.1</span><span class="p">)))</span>  <span class="c1"># adds overdispersion</span>
</span><span class="line">
</span><span class="line"><span class="c1"># data generation (on data scale!) based on Poisson distribution</span>
</span><span class="line">DataF<span class="o">$</span>Egg <span class="o"><-</span> rpois<span class="p">(</span>length<span class="p">(</span>EggLink<span class="p">),</span> exp<span class="p">(</span>EggLink<span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># save data (to current work directory)</span>
</span><span class="line">write.csv<span class="p">(</span>DataF<span class="p">,</span> file <span class="o">=</span> <span class="s">"BeetlesFemale.csv"</span><span class="p">,</span> row.names <span class="o">=</span> <span class="k-Variable">F</span><span class="p">)</span>
</span>

Having simulated a dataset, calculate the $R^2$ point-estimates, using the lme4 package to fit the model.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
<span class="line"><span class="c1">## Now, calculate R-squared (code from Nakagawa & Schielzeth 2013)</span>
</span><span class="line">library<span class="p">(</span>arm<span class="p">)</span>
</span><span class="line">library<span class="p">(</span>lme4<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Clear memory</span>
</span><span class="line">rm<span class="p">(</span>list <span class="o">=</span> ls<span class="p">())</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Read fecundity data (Poisson, available for females only)</span>
</span><span class="line">Data <span class="o"><-</span> read.csv<span class="p">(</span><span class="s">"BeetlesFemale.csv"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Creating a dummy variable that allows estimating additive dispersion in</span>
</span><span class="line"><span class="c1"># lmer This triggers a warning message when fitting the model</span>
</span><span class="line">Unit <span class="o"><-</span> factor<span class="p">(</span><span class="m">1</span><span class="o">:</span>length<span class="p">(</span>Data<span class="o">$</span>Egg<span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Fit null model without fixed effects (but including all random effects)</span>
</span><span class="line">m0 <span class="o"><-</span> lmer<span class="p">(</span>Egg <span class="o">~</span> <span class="m">1</span> <span class="o">+</span> <span class="p">(</span><span class="m">1</span> <span class="o">|</span> Population<span class="p">)</span> <span class="o">+</span> <span class="p">(</span><span class="m">1</span> <span class="o">|</span> Container<span class="p">)</span> <span class="o">+</span> <span class="p">(</span><span class="m">1</span> <span class="o">|</span> Unit<span class="p">),</span>
</span><span class="line">    family <span class="o">=</span> <span class="s">"poisson"</span><span class="p">,</span> data <span class="o">=</span> Data<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Fit alternative model including fixed and all random effects</span>
</span><span class="line">mF <span class="o"><-</span> lmer<span class="p">(</span>Egg <span class="o">~</span> Treatment <span class="o">+</span> Habitat <span class="o">+</span> <span class="p">(</span><span class="m">1</span> <span class="o">|</span> Population<span class="p">)</span> <span class="o">+</span> <span class="p">(</span><span class="m">1</span> <span class="o">|</span> Container<span class="p">)</span> <span class="o">+</span>
</span><span class="line">    <span class="p">(</span><span class="m">1</span> <span class="o">|</span> Unit<span class="p">),</span> family <span class="o">=</span> <span class="s">"poisson"</span><span class="p">,</span> data <span class="o">=</span> Data<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># View model fits for both models</span>
</span><span class="line">summary<span class="p">(</span>m0<span class="p">)</span>
</span><span class="line">summary<span class="p">(</span>mF<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Extraction of fitted value for the alternative model fixef() extracts</span>
</span><span class="line"><span class="c1"># coefficents for fixed effects [email protected] returns fixed effect design matrix</span>
</span><span class="line">Fixed <span class="o"><-</span> fixef<span class="p">(</span>mF<span class="p">)[</span><span class="m">2</span><span class="p">]</span> <span class="o">*</span> mF<span class="o">@</span>X<span class="p">[,</span> <span class="m">2</span><span class="p">]</span> <span class="o">+</span> fixef<span class="p">(</span>mF<span class="p">)[</span><span class="m">3</span><span class="p">]</span> <span class="o">*</span> mF<span class="o">@</span>X<span class="p">[,</span> <span class="m">3</span><span class="p">]</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Calculation of the variance in fitted values</span>
</span><span class="line">VarF <span class="o"><-</span> var<span class="p">(</span>Fixed<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># An alternative way for getting the same result</span>
</span><span class="line">VarF <span class="o"><-</span> var<span class="p">(</span>as.vector<span class="p">(</span>fixef<span class="p">(</span>mF<span class="p">)</span> <span class="o">%*%</span> t<span class="p">(</span>mF<span class="o">@</span>X<span class="p">)))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># R2GLMM(m) - marginal R2GLMM see Equ. 29 and 30 and Table 2 fixef(m0)</span>
</span><span class="line"><span class="c1"># returns the estimate for the intercept of null model</span>
</span><span class="line">R2m <span class="o"><-</span> VarF<span class="o">/</span><span class="p">(</span>VarF <span class="o">+</span> VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Container<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span>
</span><span class="line">               VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Population<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span> VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Unit<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span>
</span><span class="line">                log<span class="p">(</span><span class="m">1</span> <span class="o">+</span> <span class="m">1</span><span class="o">/</span>exp<span class="p">(</span>as.numeric<span class="p">(</span>fixef<span class="p">(</span>m0<span class="p">))))</span>
</span><span class="line">            <span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># R2GLMM(c) - conditional R2GLMM for full model Equ. XXX, XXX</span>
</span><span class="line">R2c <span class="o"><-</span> <span class="p">(</span>VarF <span class="o">+</span> VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Container<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span> VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Population<span class="p">[</span><span class="m">1</span><span class="p">])</span><span class="o">/</span>
</span><span class="line">         <span class="p">(</span>VarF <span class="o">+</span> VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Container<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span> VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Population<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span>
</span><span class="line">           VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Unit<span class="p">[</span><span class="m">1</span><span class="p">]</span> <span class="o">+</span> log<span class="p">(</span><span class="m">1</span> <span class="o">+</span> <span class="m">1</span><span class="o">/</span>exp<span class="p">(</span>as.numeric<span class="p">(</span>fixef<span class="p">(</span>m0<span class="p">))))</span>
</span><span class="line">         <span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Print marginal and conditional R-squared values</span>
</span><span class="line">cbind<span class="p">(</span>R2m<span class="p">,</span> R2c<span class="p">)</span>
</span>

Having stored our point estimates, we can now turn to Bayesian methods instead, and generate $R^2$ posteriors.

Posterior uncertainty in $R^2_{GLMM}$

We need to fit two models in order to get the needed parameters for $R^2_{GLMM}$. First, a model that includes all random effects, but only an intercept fixed effect is fit to estimate the distribution specific variance $\sigma^2_d$. Second, we fit a model that includes all random and all fixed effects to estimate the remaining variance components.

First I’ll clean up the data that we’ll feed to JAGS:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
<span class="line"><span class="c1"># Prepare the data</span>
</span><span class="line">jags_d <span class="o"><-</span> as.list<span class="p">(</span>Data<span class="p">)[</span><span class="o">-</span>c<span class="p">(</span><span class="m">2</span><span class="p">,</span> <span class="m">3</span><span class="p">)]</span>  <span class="c1"># redefine container, don't need sex</span>
</span><span class="line">jags_d<span class="o">$</span>nobs <span class="o"><-</span> nrow<span class="p">(</span>Data<span class="p">)</span>
</span><span class="line">jags_d<span class="o">$</span>npop <span class="o"><-</span> length<span class="p">(</span>unique<span class="p">(</span>jags_d<span class="o">$</span>Population<span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># renumber containers from 1:ncontainer for ease of indexing</span>
</span><span class="line">jags_d<span class="o">$</span>Container <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> nrow<span class="p">(</span>Data<span class="p">))</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nrow<span class="p">(</span>Data<span class="p">))</span> <span class="p">{</span>
</span><span class="line">  jags_d<span class="o">$</span>Container<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> which<span class="p">(</span>unique<span class="p">(</span>Data<span class="o">$</span>Container<span class="p">)</span> <span class="o">==</span> Data<span class="o">$</span>Container<span class="p">[</span>i<span class="p">])</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">jags_d<span class="o">$</span>ncont <span class="o"><-</span> length<span class="p">(</span>unique<span class="p">(</span>jags_d<span class="o">$</span>Container<span class="p">))</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Convert binary factors to 0's and 1's</span>
</span><span class="line">jags_d<span class="o">$</span>Habitat <span class="o"><-</span> ifelse<span class="p">(</span>jags_d<span class="o">$</span>Habitat <span class="o">==</span> <span class="s">"dry"</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">)</span>
</span><span class="line">jags_d<span class="o">$</span>Treatment <span class="o"><-</span> ifelse<span class="p">(</span>jags_d<span class="o">$</span>Treatment <span class="o">==</span> <span class="s">"Cont"</span><span class="p">,</span> <span class="m">0</span><span class="p">,</span> <span class="m">1</span><span class="p">)</span>
</span><span class="line">str<span class="p">(</span>jags_d<span class="p">)</span>
</span>

Then, fitting the intercept model:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
<span class="line"><span class="c1"># intercept model statement:</span>
</span><span class="line">cat<span class="p">(</span><span class="s">"</span>
</span><span class="line"><span class="s">model{</span>
</span><span class="line"><span class="s">  # priors on precisions (inverse variances)</span>
</span><span class="line"><span class="s">  tau.pop ~ dgamma(0.01, 0.01)</span>
</span><span class="line"><span class="s">  sd.pop <- sqrt(1/tau.pop)</span>
</span><span class="line"><span class="s">  tau.cont ~ dgamma(0.01, 0.01)</span>
</span><span class="line"><span class="s">  sd.cont <- sqrt(1/tau.cont)</span>
</span><span class="line"><span class="s">  tau.unit ~ dgamma(0.01, 0.01)</span>
</span><span class="line"><span class="s">  sd.unit <- sqrt(1/tau.unit)</span>
</span><span class="line"><span class="s">  # prior on intercept</span>
</span><span class="line"><span class="s">  alpha ~ dnorm(0, 0.01)</span>
</span><span class="line">
</span><span class="line"><span class="s">  # random effect of container</span>
</span><span class="line"><span class="s">  for (i in 1:ncont){</span>
</span><span class="line"><span class="s">    cont[i] ~ dnorm(0, tau.cont)</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line">
</span><span class="line"><span class="s">  # random effect of population</span>
</span><span class="line"><span class="s">  for (i in 1:npop){</span>
</span><span class="line"><span class="s">    pop[i] ~ dnorm(0, tau.pop)</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line">
</span><span class="line"><span class="s">  # likelihood</span>
</span><span class="line"><span class="s">  for (i in 1:nobs){</span>
</span><span class="line"><span class="s">    Egg[i] ~ dpois(mu[i])</span>
</span><span class="line"><span class="s">    log(mu[i]) <- cont[Container[i]] + pop[Population[i]] + unit[i]</span>
</span><span class="line"><span class="s">    unit[i] ~ dnorm(alpha, tau.unit) </span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line"><span class="s">}</span>
</span><span class="line"><span class="s">    "</span><span class="p">,</span> fill<span class="o">=</span><span class="k-Variable">T</span><span class="p">,</span> file<span class="o">=</span><span class="s">"pois_intercept.txt"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">nstore <span class="o"><-</span> <span class="m">2000</span>
</span><span class="line">nthin <span class="o"><-</span> <span class="m">20</span>
</span><span class="line">ni <span class="o"><-</span> nstore<span class="o">*</span>nthin
</span><span class="line">require<span class="p">(</span>rjags<span class="p">)</span>
</span><span class="line">int_mod <span class="o"><-</span> jags.model<span class="p">(</span><span class="s">"pois_intercept.txt"</span><span class="p">,</span>
</span><span class="line">                      data<span class="o">=</span>jags_d<span class="p">[</span><span class="o">-</span>c<span class="p">(</span><span class="m">2</span><span class="p">,</span> <span class="m">3</span><span class="p">)],</span> <span class="c1"># exclude unused data </span>
</span><span class="line">                      n.chains<span class="o">=</span><span class="m">3</span><span class="p">,</span>
</span><span class="line">                      n.adapt<span class="o">=</span><span class="m">5000</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">vars <span class="o"><-</span> c<span class="p">(</span><span class="s">"sd.pop"</span><span class="p">,</span> <span class="s">"sd.cont"</span><span class="p">,</span> <span class="s">"sd.unit"</span><span class="p">,</span> <span class="s">"alpha"</span><span class="p">)</span>
</span><span class="line">int_out <span class="o"><-</span> coda.samples<span class="p">(</span>int_mod<span class="p">,</span> n.iter<span class="o">=</span>ni<span class="p">,</span> thin<span class="o">=</span>nthin<span class="p">,</span>
</span><span class="line">                        variable.names<span class="o">=</span>vars<span class="p">)</span>
</span>

Then, fit the full mixed-model with all fixed and random effects:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
<span class="line"><span class="c1"># covariate model statement:</span>
</span><span class="line">cat<span class="p">(</span><span class="s">"</span>
</span><span class="line"><span class="s">model{</span>
</span><span class="line"><span class="s">  # priors on precisions (inverse variances)</span>
</span><span class="line"><span class="s">  tau.pop ~ dgamma(0.01, 0.01)</span>
</span><span class="line"><span class="s">  sd.pop <- sqrt(1/tau.pop)</span>
</span><span class="line"><span class="s">  tau.cont ~ dgamma(0.01, 0.01)</span>
</span><span class="line"><span class="s">  sd.cont <- sqrt(1/tau.cont)</span>
</span><span class="line"><span class="s">  tau.unit ~ dgamma(0.01, 0.01)</span>
</span><span class="line"><span class="s">  sd.unit <- sqrt(1/tau.unit)</span>
</span><span class="line"><span class="s">  # priors on coefficients</span>
</span><span class="line"><span class="s">  alpha ~ dnorm(0, 0.01)</span>
</span><span class="line"><span class="s">  beta1 ~ dnorm(0, 0.01)</span>
</span><span class="line"><span class="s">  beta2 ~ dnorm(0, 0.01)</span>
</span><span class="line">
</span><span class="line"><span class="s">  # random effect of container</span>
</span><span class="line"><span class="s">  for (i in 1:ncont){</span>
</span><span class="line"><span class="s">    cont[i] ~ dnorm(0, tau.cont)</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line">
</span><span class="line"><span class="s">  # random effect of population</span>
</span><span class="line"><span class="s">  for (i in 1:npop){</span>
</span><span class="line"><span class="s">    pop[i] ~ dnorm(0, tau.pop)</span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line">
</span><span class="line"><span class="s">  # likelihood</span>
</span><span class="line"><span class="s">  for (i in 1:nobs){</span>
</span><span class="line"><span class="s">    Egg[i] ~ dpois(mu[i])</span>
</span><span class="line"><span class="s">    log(mu[i]) <- cont[Container[i]] + pop[Population[i]] + unit[i]</span>
</span><span class="line"><span class="s">    mu_f[i] <- alpha + beta1 * Treatment[i] + beta2 * Habitat[i]</span>
</span><span class="line"><span class="s">    unit[i] ~ dnorm(mu_f[i], tau.unit) </span>
</span><span class="line"><span class="s">  }</span>
</span><span class="line"><span class="s">}</span>
</span><span class="line"><span class="s">    "</span><span class="p">,</span> fill<span class="o">=</span><span class="k-Variable">T</span><span class="p">,</span> file<span class="o">=</span><span class="s">"pois_cov.txt"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">cov_mod <span class="o"><-</span> jags.model<span class="p">(</span><span class="s">"pois_cov.txt"</span><span class="p">,</span>
</span><span class="line">                      data<span class="o">=</span>jags_d<span class="p">,</span>
</span><span class="line">                      n.chains<span class="o">=</span><span class="m">3</span><span class="p">,</span>
</span><span class="line">                      n.adapt<span class="o">=</span><span class="m">5000</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">vars2 <span class="o"><-</span> c<span class="p">(</span><span class="s">"sd.pop"</span><span class="p">,</span> <span class="s">"sd.cont"</span><span class="p">,</span> <span class="s">"sd.unit"</span><span class="p">,</span> <span class="s">"alpha"</span><span class="p">,</span> <span class="s">"beta1"</span><span class="p">,</span> <span class="s">"beta2"</span><span class="p">)</span>
</span><span class="line">cov_out <span class="o"><-</span> coda.samples<span class="p">(</span>cov_mod<span class="p">,</span> n.iter<span class="o">=</span>ni<span class="p">,</span> thin<span class="o">=</span>nthin<span class="p">,</span>
</span><span class="line">                        variable.names<span class="o">=</span>vars2<span class="p">)</span>
</span>

For every MCMC draw, we can calculate $R^2_{GLMM}$, generating posteriors for both the marginal and conditional values.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
<span class="line"><span class="c1"># Step 1: variance in expected values (using fixed effects only)</span>
</span><span class="line">require<span class="p">(</span>ggmcmc<span class="p">)</span>
</span><span class="line">d_int <span class="o"><-</span> ggs<span class="p">(</span>int_out<span class="p">)</span>
</span><span class="line">d_cov <span class="o"><-</span> ggs<span class="p">(</span>cov_out<span class="p">)</span>
</span><span class="line">
</span><span class="line">alpha_cov <span class="o"><-</span> subset<span class="p">(</span>d_cov<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"alpha"</span><span class="p">)</span><span class="o">$</span>value
</span><span class="line">alpha_int <span class="o"><-</span> subset<span class="p">(</span>d_int<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"alpha"</span><span class="p">)</span><span class="o">$</span>value
</span><span class="line">b1_cov <span class="o"><-</span> subset<span class="p">(</span>d_cov<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"beta1"</span><span class="p">)</span><span class="o">$</span>value
</span><span class="line">b2_cov <span class="o"><-</span> subset<span class="p">(</span>d_cov<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"beta2"</span><span class="p">)</span><span class="o">$</span>value
</span><span class="line">
</span><span class="line">Xmat <span class="o"><-</span> cbind<span class="p">(</span>rep<span class="p">(</span><span class="m">1</span><span class="p">,</span> jags_d<span class="o">$</span>nobs<span class="p">),</span> jags_d<span class="o">$</span>Treatment<span class="p">,</span> jags_d<span class="o">$</span>Habitat<span class="p">)</span>
</span><span class="line">beta_mat <span class="o"><-</span> cbind<span class="p">(</span>alpha_cov<span class="p">,</span> b1_cov<span class="p">,</span> b2_cov<span class="p">)</span>
</span><span class="line">
</span><span class="line">fixed_expect <span class="o"><-</span> array<span class="p">(</span>dim <span class="o">=</span> c<span class="p">(</span>nstore<span class="p">,</span> jags_d<span class="o">$</span>nobs<span class="p">))</span>
</span><span class="line">varF <span class="o"><-</span> rep<span class="p">(</span><span class="kc">NA</span><span class="p">,</span> nstore<span class="p">)</span>
</span><span class="line"><span class="kr">for</span> <span class="p">(</span>i <span class="kr">in</span> <span class="m">1</span><span class="o">:</span>nstore<span class="p">)</span> <span class="p">{</span>
</span><span class="line">    fixed_expect<span class="p">[</span>i<span class="p">,</span> <span class="p">]</span> <span class="o"><-</span> beta_mat<span class="p">[</span>i<span class="p">,</span> <span class="p">]</span> <span class="o">%*%</span> t<span class="p">(</span>Xmat<span class="p">)</span>
</span><span class="line">    varF<span class="p">[</span>i<span class="p">]</span> <span class="o"><-</span> var<span class="p">(</span>fixed_expect<span class="p">[</span>i<span class="p">,</span> <span class="p">])</span>
</span><span class="line"><span class="p">}</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Step 2: calculate remaining variance components </span>
</span><span class="line"><span class="c1"># among container variance</span>
</span><span class="line">varCont <span class="o"><-</span> subset<span class="p">(</span>d_cov<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"sd.cont"</span><span class="p">)</span><span class="o">$</span>value<span class="o">^</span><span class="m">2</span>
</span><span class="line"><span class="c1"># among population variance</span>
</span><span class="line">varPop <span class="o"><-</span> subset<span class="p">(</span>d_cov<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"sd.pop"</span><span class="p">)</span><span class="o">$</span>value<span class="o">^</span><span class="m">2</span>
</span><span class="line"><span class="c1"># overdispersion variance</span>
</span><span class="line">varUnit <span class="o"><-</span> subset<span class="p">(</span>d_cov<span class="p">,</span> Parameter <span class="o">==</span> <span class="s">"sd.unit"</span><span class="p">)</span><span class="o">$</span>value<span class="o">^</span><span class="m">2</span>
</span><span class="line"><span class="c1"># distribution variance (Table 2, Nakagawa & Schielzeth 2013)</span>
</span><span class="line">varDist <span class="o"><-</span> log<span class="p">(</span><span class="m">1</span><span class="o">/</span>exp<span class="p">(</span>alpha_int<span class="p">)</span> <span class="o">+</span> <span class="m">1</span><span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># Finally, calculate posterior R-squared values </span>
</span><span class="line"><span class="c1"># marginal</span>
</span><span class="line">postR2m <span class="o"><-</span> varF<span class="o">/</span><span class="p">(</span>varF <span class="o">+</span> varCont <span class="o">+</span> varPop <span class="o">+</span> varUnit <span class="o">+</span> varDist<span class="p">)</span>
</span><span class="line"><span class="c1"># conditional</span>
</span><span class="line">postR2c <span class="o"><-</span> <span class="p">(</span>varF <span class="o">+</span> varCont <span class="o">+</span> varPop<span class="p">)</span><span class="o">/</span>
</span><span class="line">             <span class="p">(</span>varF <span class="o">+</span> varCont <span class="o">+</span> varPop <span class="o">+</span> varUnit <span class="o">+</span> varDist<span class="p">)</span>
</span><span class="line">
</span><span class="line"><span class="c1"># compare posterior R-squared values to point estimates</span>
</span><span class="line">par<span class="p">(</span>mfrow <span class="o">=</span> c<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">2</span><span class="p">))</span>
</span><span class="line">hist<span class="p">(</span>postR2m<span class="p">,</span> main <span class="o">=</span> <span class="s">"Marginal R-squared"</span><span class="p">,</span>
</span><span class="line">        ylab <span class="o">=</span> <span class="s">"Posterior density"</span><span class="p">,</span>
</span><span class="line">        xlab <span class="o">=</span> <span class="kc">NULL</span><span class="p">,</span> breaks <span class="o">=</span> <span class="m">20</span><span class="p">)</span>
</span><span class="line">abline<span class="p">(</span>v <span class="o">=</span> R2m<span class="p">,</span> col <span class="o">=</span> <span class="s">"blue"</span><span class="p">,</span> lwd <span class="o">=</span> <span class="m">4</span><span class="p">)</span>
</span><span class="line">hist<span class="p">(</span>postR2c<span class="p">,</span> main <span class="o">=</span> <span class="s">"Conditional R-squared"</span><span class="p">,</span>
</span><span class="line">        ylab <span class="o">=</span> <span class="s">"Posterior density"</span><span class="p">,</span>
</span><span class="line">        xlab <span class="o">=</span> <span class="kc">NULL</span><span class="p">,</span> breaks <span class="o">=</span> <span class="m">25</span><span class="p">)</span>
</span><span class="line">abline<span class="p">(</span>v <span class="o">=</span> R2c<span class="p">,</span> col <span class="o">=</span> <span class="s">"blue"</span><span class="p">,</span> lwd <span class="o">=</span> <span class="m">4</span><span class="p">)</span>
</span>

This plot shows the posterior $R^2_{GLMM}$ distributions for both the marginal and conditional cases, with the point estimates generated with lmer shown as vertical blue lines. Personally, I find it to be a bit more informative and intuitive to think of $R^2$ as a probability distribution that integrates uncertainty in its component parameters. That said, it is unconventional to represent $R^2$ in this way, which could compromise the ease with which this handy statistic can be explained to the uninitiated (e.g. first year biology undergraduates). But, being a derived parameter, those wishing to generate a posterior can do so relatively easily.

Aside: correspondence between parameter estimates

Some may be wondering whether the parameter estimates generated with lme4 are comparable to those generated using JAGS. Having used vague priors, we would expect them to be similar. We can plot the Bayesian credible intervals (in blue), with the previous point estimates (as open black circles):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
<span class="line">par<span class="p">(</span>mfrow <span class="o">=</span> c<span class="p">(</span><span class="m">1</span><span class="p">,</span> <span class="m">2</span><span class="p">))</span>
</span><span class="line">require<span class="p">(</span>mcmcplots<span class="p">)</span>
</span><span class="line">caterplot<span class="p">(</span>int_out<span class="p">,</span> style <span class="o">=</span> <span class="s">"plain"</span><span class="p">)</span>
</span><span class="line">caterpoints<span class="p">(</span>c<span class="p">(</span>m0<span class="o">@</span>fixef<span class="p">,</span>
</span><span class="line">              attr<span class="p">(</span>VarCorr<span class="p">(</span>m0<span class="p">)</span><span class="o">$</span>Container<span class="p">,</span> <span class="s">"stddev"</span><span class="p">),</span>
</span><span class="line">              attr<span class="p">(</span>VarCorr<span class="p">(</span>m0<span class="p">)</span><span class="o">$</span>Population<span class="p">,</span> <span class="s">"stddev"</span><span class="p">),</span>
</span><span class="line">              attr<span class="p">(</span>VarCorr<span class="p">(</span>m0<span class="p">)</span><span class="o">$</span>Unit<span class="p">,</span> <span class="s">"stddev"</span><span class="p">)))</span>
</span><span class="line">title<span class="p">(</span><span class="s">"Intercept model BCIs & point estimates"</span><span class="p">)</span>
</span><span class="line">
</span><span class="line">caterplot<span class="p">(</span>cov_out<span class="p">,</span> style <span class="o">=</span> <span class="s">"plain"</span><span class="p">)</span>
</span><span class="line">caterpoints<span class="p">(</span>c<span class="p">(</span>mF<span class="o">@</span>fixef<span class="p">,</span>
</span><span class="line">              attr<span class="p">(</span>VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Container<span class="p">,</span> <span class="s">"stddev"</span><span class="p">),</span>
</span><span class="line">              attr<span class="p">(</span>VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Population<span class="p">,</span> <span class="s">"stddev"</span><span class="p">),</span>
</span><span class="line">              attr<span class="p">(</span>VarCorr<span class="p">(</span>mF<span class="p">)</span><span class="o">$</span>Unit<span class="p">,</span> <span class="s">"stddev"</span><span class="p">)))</span>
</span><span class="line">title<span class="p">(</span><span class="s">"Covariate model BCIs & point estimates"</span><span class="p">)</span>
</span>

To leave a comment for the author, please follow the link and comment on their blog: Ecology in silico.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)