R^2 for linear mixed effects models

Linear mixed effects models are a powerful technique for the analysis of ecological data, especially in the presence of nested or hierarchical variables. But unlike their purely fixed-effects cousins, they lack an obvious criterion to assess model fit.

[Updated October 13, 2015: Development of the R function has moved to my piecewiseSEM package, which can be found here under the function sem.model.fits]

In the fixed-effects world, the coefficient of determination, better known as R2, is a useful and intuitive tool for describing the predictive capacity of your model: its simply the total variance in the response explained by all the predictors in your model.

In a least squares regression, R2 is the sum of differences in the observed minus the fitted values, over the sum of the deviation of each observation from the mean of the response, all subtracted from 1:


For models estimated using maximum likelihood, this equation was generalized to the ratio of the likelihood of the intercept-only model over the full model, raised to the 2/N, and all subtracted from 1:


Both equations can be interpreted identically.

R2 has a number of useful properties. Its independent of sample size, bound (0,1), and dimensionless, which makes it ideal for comparing fits across different datasets.

I should note, however, that its a poor tool for model selection, since it almost always favors the most complex models. If the goal is to select among the best models, an information criterion approach (such as AIC or BIC) is preferred, since these indicators penalize for the number of predictors.

Unfortunately (and this is a common pitfall), the best model as determined by AIC is not synonymous with a good model. In other words, AIC provides a great test of relative model fit, but not absolute model fit. AIC values are also not comparable across models fit to different datasets. Thus, R2 and AIC both have their place in ecological statistics.

Often, researchers using mixed models report an R2 from a linear mixed model as simply the squared correlation between the fitted and observed values (see here), but this is a pseudo-R2 and is technically incorrect.

Why? Because a mixed effects model yields a variance associated with each random factor and the residual variance, so its not entirely clear which to use when calculating the R2 (the approach above ignores this issue by choosing to calculate R2 relative to only the residual variance, which omits any structure in the data imposed by the random factor).

When hierarchical models with >1 levels are considered, this problem is exacerbated since now we are dealing with variances at multiple levels, plus the residual error. Some people have discussed calculating an R2 value for each level of the random factor, but this can lead to negative R2 values when addition of predictors reduces the residual error while increasing the variance of the random component (or vice versa), even though the sum of the variance components remains unchanged.

So what have we done, historically-speaking? We’ve compared AIC values for the full and null models, or calculated all manner of pseudo-R2s, and quietly (and perhaps for some, ashamedly) went about our business.

Which is why I was so happy to see a paper with the title “A general and simple method for obtaining R2 from generalized linear mixed-effects models” by Shinichi Nakagawa and Holger Schielzeth appear in my newsfeed. Basically, they have derived two easily interpretable values of R2 that address the above issues, specifically that of negative pseudo-R2s with more predictors, while still honoring the random structure of the data.

The first is called the marginal R2 and describes the proportion of variance explained by the fixed factor(s) alone:


The fixed-effects variance is in the numerator. The denominator is the total variance explained by the model, including (in order): the fixed-effects variance, the random variance (partitioned by level l), and the last two terms add up the residual variance and are the additive dispersion component (for non-normal models) and the distribution-specific variance.

The second is the conditional R2, which describes the proportion of variance explained by both the fixed and random factors:


In this case, the numerator contains both the variance of the fixed effects, as well as the sum of random variance components for each level l of the random factor. You’ll note the denominator is identical.

This approach was later modified by Johnson (2014) to include random slope models.

While this may not be the holy grail of fit statistics (and the authors concede on this point as well), in my opinion, it does address many of the problems of pseudo-R2s and so is going to be my preferred statistic.

The authors also note (and I echo) that looking at R2 and AIC values is not a replacement for testing basic model assumptions, such as normality of predictors and non-constant variance, and should be used in tandem.

I’ve written a function in R called sem.model.fits (along with an amiable fellow named Juan Sebastian Casallas) that calculates the marginal and condition R2s, as well as AIC values for a complementary approach to model comparison. The input is either single model, or a list of models of most common classes, including models from the lme4 and nlme packages. The output is a data.frame where each row is a model in the same order as in the input list, with the model class, marginal and conditional R2s, and the AIC value as columns if the response is identical across all models.

You can find an example of this function in action here: https://github.com/jslefche/piecewiseSEM/blob/master/README.md#get-r2-for-individual-models

I’ve since incorporated this function into my piecewiseSEM package for piecewise structural equation modeling. Since the function found here will no longer be actively developed by me, I recommend using the function in the package to minimize the potential for issues.

[Edit: Also check out this awesome tutorial for quantifying uncertainty around these estimates!]


Nakagawa, S., and H. Schielzeth. 2013. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4(2): 133-142. DOI: 10.1111/j.2041-210x.2012.00261.x

Johnson, Paul C.D. 2014. Extension of Nakagawa & Schielzeth’s R2GLMM to random slopes models. Methods in Ecology and Evolution. DOI: 10.1111/2041-210X.12225.

199 thoughts on “R^2 for linear mixed effects models

  1. Bastien Ferland-Raymond says:


    Thanks a lot for sharing your function, very usefull! However, I’ve found a little annoyance (bug) that I would like to point out. In some situation, i.e. when the “data” object doesn’t exist in the global environment, the function fails with: Error: ‘data’ not found, and some variables missing from formula environment

    This is (I think) caused by the update fonction you use at line 70 (for sure) and maybe also 90 and 133. It seems the update() function goes to look for the data in the global environment and not in the model object (“mdl”). I’ve fixed the problem by changing the command with: mdl.aic = AIC(update(mdl, method=”ML”, data=mdl@frame))

    You may want to check into that.

    I had the problem because I was running multiple lmer() on a list of dataframe with an homemade function. My model data was therefore named with the name of my function argument and not of the original dataset.

    Thanks again for the good work!

  2. srecko says:

    Hi Jon,
    Thanks for posting this, it’s really helpful. I have one question regarding the results I got. I built a model (lme) that includes a count variable as a dependent variable and 2 predictors – condition graded/non_graded, and experience – yes/no, as well as the interaction between those two. In the first case, I specified “student” as a random effect and r squared was:
    Class Family Link Marginal Conditional AIC
    lme gaussian identity 0.3964543 0.4609847 983.9265

    In the next step I introduced a “student within a course” as a random effect. That is, instead of student only, now I have student within a course. The model yielded a better fit according to AICc, but r squared is:
    Class Family Link Marginal Conditional AIC
    lme gaussian identity 0.1058472 0.9387431 966.9221

    There is a significant drop for marginal r squared, and even bigger increase for the conditional r squared. Is this even possible? 🙂 I think I know how to interpret this, but honestly, haven’t seen something like this before.


    • Jon Lefcheck says:

      Hi Srecko, I have run across this phenomenon before, and apparently it is not uncommon. Depending on how the variance is ascribed among and between the hierarchical levels of your random structure, R2 can go up or down. See here. I suppose you will have to balance the trade-off in explanatory power of your fixed variables vs. better specification of your random structure. If you are looking to make predictions, I would go with the model that ascribes more variance to your fixed effects (e.g., higher marginal R2). If you are looking to identify important predictors I would go with the model that explains the most variance (e.g., higher conditional R2). As long as you can justify either/or I think you should be ok. Cheers, Jon

  3. Helen says:

    Hi Jon, thanks very much for sharing this helpful code. I wonder what would be the appropriate way to cite my use of this code in my thesis?

  4. James Robinson says:

    Fantastic work with the functions – and thanks for sharing.

    I’m using glmers for a Gamma distribution with log link, and would love to be able to estimate an R2. Is there a reason that you didn’t include this family in your code?


    • Jon Lefcheck says:

      Hi James, other distributions were not included because they weren’t in the papers by Schielzeth & Nakagawa or Johnson. I’m not sure what the distribution-specific variance would be (or if anyone knows), so I’m afraid I can’t update the function — sorry! Cheers, Jon

  5. Johan Craeymeersch says:

    Is it possible that the marginal and conditional r squared are the same in a mixed model? I get that situation a few times when adding a spatial autocorrelation structure to the model. Without this structure I get different r squares, e.g.:
    Class Marginal Conditional AIC
    1 lme 0.0998046 0.1640069 9804.341

    After adding spatial autocorrelation structure, I get:
    Class Marginal Conditional AIC
    1 lme 0.08996415 0.08996415 9780.116

    • Paul Johnson says:

      If marginal R2 = conditional R2 then the random effect variances included in the conditional R2 must sum to zero. It would help to know how you specified the autocorrelation structure – could you post the code that fits the model? Did you apply the correlation structure to an existing random effect, or add a separate spatial random effect? Your results could be explained by the variance of the spatial random effect not being included by the function. You could check this by (1) checking for a non-zero random effect variance in the model summary and (2) stepping through the function code to find out why it isn’t being picked up.

      • Johan Craeymeersch says:

        The mixed model without spatial autocorrelation is following, where Diepte2, Mediaan, VarSalsummer2, AvgTauGolf2, MaxVel and T5_Winter are abiotic conditions, locatienr is the station name and deelgebied is geographic area.
        lme1<-lme(Lmax2~Diepte2+Mediaan+VarSalSummer2+AvgTauGolf2+MaxVel+T5_Winter, data=Biom,random=~1|as.factor(locatienr),weights=varIdent(form=~1|factor(deelgebied)))

        I added spatial autocorrelation as follows:
        tmp1E<-update(lme1,correlation=corExp(1,form=~"X-Parijs gerealiseerd"+"Y-Parijs gerealiseerd",nugget=T), control=list(msMaxIter = 100))
        where X-Parijs … and Y-Parijs are the coordinates.

      • Jon Lefcheck says:

        Thanks Paul. Johan, I can reproduce this (code below), but I’m not quite sure what is going on. I will need to dig into this more. It might be worth reaching out to Cross validated and seeing if they can figure out why the fixed = random R2 with the addition of a spatial correlation structure.

        # Read function
        r.squared = source_url("https://raw.githubusercontent.com/jslefche/rsquared.glmm/master/rsquaredglmm.R")[[1]]

        # Create dataset

        data = data.frame(
        y = rnorm(20, 0, 1),
        x = rnorm(20, 0, 1),
        location = letters[1:10],
        lat = runif(20, 0, 20),
        long = runif(20, -90, -70)

        # Add small amount to lat long for last set of observations
        # data[11:20, "lat"] = data[11:20, "lat"] + runif(10, 0, 5)
        # data[11:20, "long"] = data[11:20, "long"] + runif(10, 0, 5)

        # Build models

        mod0 = lme(y ~ 1, random = ~ 1 | location, data)

        mod1 = update(mod0, fixed = y ~ x)

        mod2 = update(mod1, correlation = corExp(form = ~ lat + long))

  6. Beth B says:

    Hello, thanks for this. I wish to calculate the conditional and marginal R2 for a linear mixed effects model. I understand that the code you provide creates a function which I can then apply to my model (‘r.squared.lme’), but I am a newbie to R and I am not which parts of the code I need to run. I have tried all and then various sections – but I just keep getting error messages when I try to apply the function. Sorry for the very basic question – but any advice on how to apply this for a beginner would be very helpful. Thank you. Beth

  7. Frank says:

    Hi John, I’m using glmms for a gaussian distribution with log link, and would love to be able to estimate an R2.
    Example: glmer(Bioma~Fungi+Par+Fungi*Parcela+(1|Point),data=camp, family=gaussian(link=”log”))

      • Maite says:

        Hi Jon and Frank. I have exactly the same problem. After apply the sem.model.fits to a list of generalised mixed models (gamma distribution, link=log), the warming message says that: “Model link ‘ log ‘ is not yet supported for the Gamma distribution”. Neither it seems to support the inverse or the indentity link function. Can you suggest any alternative to calculate R2 for glmms of the Gamma family? Many thanks for your blog!

  8. Magali Paquot (@MagaliPaquot) says:

    Dear John,
    I am a beginner in multilevel analysis but very excited about what these methods offer.
    I’m trying to use the r.squaredGLMM function on the following model but get NA values:

    > yearRI r.squaredGLMM(yearRI)
    R2m R2c
    NA NA

    I tried and deleted the method=”ML” but apparently this is not the problem. My problem is certainly very basic but I would appreciate any help! The r.squaredLR(yearRI) returns results but I am not sure I understand the function settings and its results. I’d rather use R^2 if possible.


      • Juan says:

        Hello Jon, unfortunately piecewiseSEM is not available for the latest R version 3.2.2. I would like to use the sem.model.fits to estimate a pseudo R2 for a lme model. If you have any solution to this R version I would really appreciate. Thanks in advance, JM

      • Jon Lefcheck says:

        Hi Juan, That’s odd — I built it in 3.2.2! Are you trying to install from CRAN? The package is not yet up there (working on it though). If you are still running into problems email me and we can get it sorted. As a last ditch you can always download the .R function from GitHub and use load() to bring the function into R. Cheers, Jon

      • Juan says:

        Hi Jon, thanks for your quick reply.

        I loaded the piecewiseSEM package from GitHub, however, I have some troubles to run the list function on piecewiseSEM. I ran lme models with and without the weights argument without success. The error message was: Error in chol.default((value + t(value))/2) : the leading minor of order 1 is not positive definite.

        Could you tell me what this error means?

        Thank you!

  9. Daniel Padfield says:

    Hi Jon

    I think this function is super cool and am just trying to apply it in response to some referee’s comments. I was just wondering what happens if you change the structure of the variance? My best model has a power variance structure (easily specified in nlme) to deal with heteroscedasticity but all of the R squares after i specify this in the model come out as 1 or highly close to it. I imagine that the modification of the variance structure is invalidating the way the R squared is calculated somehow?

    Any thoughts would be much appreciated.
    Daniel Padfield

    • Jon Lefcheck says:

      Hi Daniel, Yes I would imagine modeling the variance would lead to unity, although I’m at a loss for the precise mathematics. Logically it makes sense, so I would proceed with caution! Cheers, Jon

      • Paul Johnson says:

        You could pull out the variance components and calculate R^2 manually. Wherever there’s a variance that in your model is allowed to vary among observations, whether due to heteroscedasticity in the residual variance or random slopes or whatever, you should take this into account by calculating that variance for each row of your data, then take the mean to get the variance component to plug into the R^2 equation. This also applies to other statistics based on variance components, such as ICC/reliability and heritability. This approach also allows you to calculate different R^2 (or ICC etc) values for different types of individual. A potentially interesting consequence of heteroscedasticity is that the response will be more predictable (more “explained”) for some types of individual/observation than others.

      • Daniel Padfield says:

        Hi Jon and Paul.

        Thanks for your responses I will try and implement them in future when asked to do so. Managed to get the paper published without having to solve this issue, but nice to know I have an option to work it out in the future!

        Many thanks

  10. Paul says:

    Hi John,
    I’m a month or so behind everyone else here, but thanks for putting this together. I have a quick question about the AIC value in your output. Its different than the AIC output given in the summary of the model (for both mixed effects models using lme and lmer). Am I doing something wrong or missing something?

    • Jon Lefcheck says:

      Hi Paul, Great question! The output is probably different because of the model estimation employed. Its only valid to compare models using maximum-likelihood (ML) estimation, particularly when the fixed effects change. Since I imagine this function will be used to conduct model comparisons, I have automatically refit models using ML instead of the default, which is restricted maximum-likelihood (REML) for most functions in R — I believe this is noted briefly in the documentation. I probably should provide a switch for this but I think this is the safest route at the moment. Hope this solves it! If not, drop me an email. Cheers, Jon

      • Paul says:

        That’s great Jon, thanks for clearing that up, it makes total sense. Sorry for spelling your name wrong in the previous comment. Thanks again.

  11. Bob Mauck says:

    I, too, thank you for posting this function. It has been extremely helpful. However, I have had an issue that doesn’t make sense mathematically and I’m wondering if you know of a problem like this.

    I am looking at weather effects on avian reproductive success. I have a “Null” mixed effects model that contains both random and fixed effects that are unrelated to weather (#1 below). Then, I add weather effects to the model (#2 below). I use your function to look at the change in the ratio of marginal R2 to conditional R2 to understand the effect of adding weather effects to the “Null” model (#1 below). The results for 5 of 6 investigations of this type have made complete sense mathematically (and biologically).

    However, in one comparison of models, the Conditional R-squared decreases when the extra fixed effects are added to the model. Here are the results:

    Class Family Link Marginal Conditional AIC
    1 glmerMod binomial logit 0.03633955 0.2647656 1997.570
    2 glmerMod binomial logit 0.05508496 0.1541642 1617.714

    The first model has 3 random effects and 2 fixed effects. The second has 3 random effects and 6 fixed effects. Both models use the same dataset (N ~ 1100). You can consider #2 a FULL model and #1 a REDUCED model for this example.

    How can the Conditional R-sq be lowered from #1 to #2 when the only difference is that #2 has added fixed effects to model #1?

    NOTE: neither model produced any warnings when they were run using glmer from lme4.

    As I said, I have done this with many other comparisons in the same investigation and they all make sense (i.e., Conditional R-sq never lowers with the addition of variables). I have looked at my code and had others look at it and I can’t see where I made a mistake. Any ideas and thanks for the help?

    • Jon Lefcheck says:

      Hmm, that is quite strange. It would be helpful if you could show me the code you used to generate the models. Also, if you used r.squaredGLMM from the MuMIn package, do you get the same issue? Thanks, Jon

      • Bob Mauck says:


        Thanks for the quick reply. I see you are at W&M in Marine Science. You must know Jon Allen. I used to run the Bowdoin Scientific Station at Kent Island and got to know Jon there. Give him my best. I last saw him at a SICB meeting a few years ago.

        As for the models. Here is the code:

        ## make a NULL model
        jdcN_WxNull <-glmer(FS ~ 1 + J_Lay_Date + BroodSize + #### intrinsic factors – fixed

        ## add the random effects

        (1|Female_ID) + (1|Year) + (1|Study_Site),

        ### add the other parameters

        data = JDC_Baseline, family = binomial, control= glmerControl(tolPwrss=1e-5, optimizer="bobyqa") )

        ## add weather effects to the NULL model
        jdcN_WxMeanIncAll <-glmer(FS ~ 1 + J_Lay_Date + BroodSize + #### intrinsic factors – fixed

        ## entire time mean weather fixed effects
        MeanRainInc + MeanWindInc + MeanVisInc + MeanTempMaxInc +

        ## add the random effects
        (1|Female_ID) + (1|Year) + (1|Study_Site),

        ### add the other parameters
        data = JDC_Baseline, family = binomial, control= glmerControl(tolPwrss=1e-5, optimizer="bobyqa") )

        I did not know about the MuMin package. I will try that.



      • Jon Lefcheck says:

        Hi Bob, Will do. I see Bowdoin is sponsoring MBEM this year — will I see you this week??

        Your code is a bit confusing as your s-called ‘intercept only’ model still has two fixed terms: J_Lay_Date + BroodSize.

        However, the second model still should not have a lower R2.

        I made up some ‘fake’ data and ran your exact code. You can find it here:

        That gives higher values for the second model.

        So I imagine it is something to do with your data?

        If you want to drop me an email we can dig into it further.



  12. Gitu wa Mbui says:

    Hi, does anyone know whether this function canbe applied to the mixed models of the GAM family – say models fonctructed using gamm4/gamm packages?


  13. Gitu wa Mbui says:

    Can it be applied to GLMM’s of non-gaussian distribution, say GLMM fitted using binomial?…the function returns results wthout warnings or errors. However I am not sure whether the distribution has to be of Gaussian family for the R2 to be accurate

  14. Alyse Larkin says:

    Hi Jon,

    I am unsure why I am getting the error below. Others have gotten the same error in the past, but the solution was to make sure that you have fixed effects. I have a fixed effect and still get the error.

    R Code:
    l.modlist <- list(lmer(r.ratio ~ latitude + (1 + season | OTU),data=l.dataFR))

    Linear mixed model fit by maximum likelihood ['lmerMod']
    Formula: r.ratio ~ latitude + (1 + season | OTU)
    Data: l.dataFR
    AIC BIC logLik deviance df.resid
    46246.91 46284.67 -23117.45 46234.91 3994
    Random effects:
    Groups Name Std.Dev. Corr
    OTU (Intercept) 54.11
    season 11.37 -0.56
    Residual 74.95
    Number of obs: 4000, groups: OTU, 100
    Fixed Effects:
    (Intercept) latitude
    81.46 -34.22

    Error in X[, rownames(Sigma), drop = FALSE] : subscript out of bounds

    I installed piecewiseSEM from CRAN today so my version should be fairly up to date. Thank you so much for your help!



    • Raphaël says:

      Hi Jon,

      I’m getting the same error as Alyse and MDW when using the sem.model.fits function (Error in X[, rownames(Sigma), drop = FALSE] : subscript out of bounds). Did you find any solution?
      Could it be because my random effects are nested. Here’s my model’s formula:

      lmer( lBAI ~ sizeE + mixE + BA:mixE + compet + DC + DC:compet + DC:mixE + Tan + Pan + DD + DD:sizeE + DD:mixE + DD:compet + DD:DC + DD:Tan + DD:Pan + (sizeE + mixE + BA:mixE + DC +DC:compet + DC:mixE + Tan + Pan + DD | ID_PLOT/ ID_TREE) )

      R and the piecewiseSEM package are up to date.
      Thank you in advance for your help.



      • Jon Lefcheck says:

        Hi Raphael, Wow that is some model you have there! Is there a reason you need so many random slopes to vary? Try updating to the development version of `piecewiseSEM` on GitHub. If that doesn’t help, send me an email. I suspect it may have to do with the presence of Observation-level random effects. Cheers, Jon

  15. John Bruno says:

    Thank you Jon for this great post and code.

    Im trying to calculate R2 for mixed effect, glmer models, e.g.,

    modelF sem.model.fits(modelF)
    Error in FUN(X[[i]], …) :
    Some variance components equal zero. Respecify random structure!

    Is my model not compatible?

    Also tried MuMIn, r.squaredGLMM(modelF) and got:

    > r.squaredGLMM(modelF)
    Error in glmer(formula = TotalCoral * 0.01 ~ scale(log10(hii50km + 1)) + :
    fitting model with the observation-level random effect term failed. Add the term manually
    In addition: Warning message:
    In t(mm[!is.na(ff), ]) :
    error in evaluating the argument ‘x’ in selecting a method for function ‘t’: Error in mm[!is.na(ff), ] : (subscript) logical subscript too long


      • Tom Baker says:

        Hi John

        I have had the same variance components equal zero. Respecify random structure! error. Did you find the cause for this in the above circumstance and was there a solution?


      • Jon Lefcheck says:

        Hi Tom, Try updating to the latest dev version of piecewiseSEM:


        And see if the issue persists. Cheers, Jon

  16. Lia says:


    I am modelling some count data using MASS and glm.nb (mixed effect model). I came across your package and was interested in using the marginal and conditional Rsq. However when I install the package and use the sem.model.fits(model1).

    It returns to me:
    Class Family Link n R.squared
    1 negbin Negative Binomial log 56 0.7706424

    This seems to be the same value for pseudo-Rsq as I have calculated using 1-(residual deviance/null deviance). I am not getting the marginal and conditional values. Is this something to do with the MASS modelling procedure?

    Thanks in advance!


    • Jon Lefcheck says:

      Hi Lia, Since glm.nb does not have random effects, there will be no conditional R-squared, only a marginal R-squared (fixed effects only). I believe I calculate the R-squared there using a different calculation of pseudo-R-squared for non-normal models. See the ?sem.model.fits help for a description (can’t recall off the top of my head). Hope this helps! Cheers, Jon

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s