Effect sizes in mixed models

[wcornwell]

This is a pretty general problem for all linear models, so someone's got to have thought about this. Be good to find current guidance...

[wcornwell]

see https://www.theanalysisfactor.com/effect-size/ and http://www.theanalysisfactor.com/calculate-effect-size/

I propose we try Eta Squared and Partial Eta Squared

thoughts @johnWilshire ?

[wcornwell]

some more issues for mixed models discussed here: https://www.journalofcognition.org/articles/10.5334/joc.10/

[wcornwell]

seems like solution is in Shinichi's paper (see equation 27)(https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/j.2041-210x.2012.00261.x)

but will take some time to work this out

[wcornwell]

someone needs to turn this paper into an R library...

[coreytcallaghan]

I've used his R2 calculation before: MuMIn::r.squaredGLMM I think is based on Shinichi's paper.

An R2 is helpful for each model. I previously calculated R2 for the top models and then averaged across the top models. Which we could do.

But, you're saying that if we pull out the equation 27 part of it, then that might represent the 'effect size' for each fixed effect?

I'll look at the function and see if I can't get anywhere.

[wcornwell]

yes, that equation is for the total variance explained by all fixed effects, but if we look at the components, they should be the variance explained by each fixed effect.

[wcornwell]

shinichi says this should work with phylo models, so seems promising

[wcornwell]

should actually be pretty straightforward to modify this from Shinichi's supplementary material:

#A general and simple method for obtaining R2 from generalized linear mixed-effects models
#Shinichi Nakagawa1,2 and Holger Schielzeth3
#1 National Centre of Growth and Development, Department of Zoology, University of Otago, Dunedin, New Zealand
#2 Department of Behavioral Ecology and Evolutionary Genetics, Max Planck Institute for Ornithology, Seewiesen, Germany
#3 Department of Evolutionary Biology, Bielefeld University, Bielefeld, Germany
#Running head: Variance explained by GLMMs
#Correspondence:
#S. Nakagawa; Department of Zoology, University of Otago, 340 Great King Street, Dunedin, 9054, New Zealand
#Tel:   +64 (0)3 479 5046
#Fax:   +64 (0)3 479 7584
#e-mail: shinichi.nakagawa@otago.ac.nz 

####################################################
# A. Preparation
####################################################
# Note that data generation appears below the analysis section.
# You can use the simulated data table from the supplementary files to reproduce exactly the same results as presented in the paper.

# Set the work directy that is used for rading/saving data tables
# setwd("/Users/R2")

# load R required packages
# If this is done for the first time, it might need to first download and install the package
# install.package("arm")
library(arm)
# install.package("lme4")
library(lme4)

####################################################
# B. Analysis
####################################################

# 1. Analysis of body size (Gaussian mixed models)
#---------------------------------------------------

# Clear memory
rm(list = ls())

# Read body length data (Gaussian, available for both sexes)
Data <- read.csv("BeetlesBody.csv")

# Fit null model without fixed effects (but including all random effects)
m0 <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data)

# Fit alternative model including fixed and all random effects
mF <- lmer(BodyL ~ Sex + Treatment + Condition + (1 | Population) + (1 | Container), data = Data)

# View model fits for both models
summary(m0)
summary(mF)

# Extraction of fitted value for the alternative model
# fixef() extracts coefficents for fixed effects
# mF@X returns fixed effect design matrix
Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3] + fixef(mF)[4] * mF@X[, 4]

# Calculation of the variance in fitted values
VarF <- var(Fixed)

# An alternative way for getting the same result
VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))

# R2GLMM(m) - marginal R2GLMM
# Equ. 26, 29 and 30
# VarCorr() extracts variance components
# attr(VarCorr(lmer.model),'sc')^2 extracts the residual variance
VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + attr(VarCorr(mF), "sc")^2)

[wcornwell]

in the help file for MuMIn::r.squaredGLMM they say they basically copied this anyway, but we have to re-write to take a phylo model object...

[coreytcallaghan]

Awesome.

Sounds like you have made some progress!

What would the difference be in order to take the phylo model?

Or, will the current function not work with the phylolm model ouput?

Currently away for the weekend, but will get back into this next week.

• Corey

On Thu, Mar 29, 2018, 5:01 PM Will Cornwell notifications@github.com wrote:

in the help file for MuMIn::r.squaredGLMM they say they basically copied this anyway, but we have to re-write to take a phylo model...

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/cornwell-lab-unsw/bird_urbanness/issues/17#issuecomment-377130908, or mute the thread https://github.com/notifications/unsubscribe-auth/Aa0iJpv9z0UlNABFs3y40XWCJLolzfAKks5tjHi1gaJpZM4S99Qg .

[wcornwell]

in principle this is sorted out. Had to move to lme4verse to use Shinichi's approach because we needed an estimate of variance explained by random factors and I was not brave enough to try to hack that out of one of the dedicated phylo packages--too much risk of me messing up either a sqrt or a ^2.

Turns out that bolker lab is currently working on sorting out phylo applications of lme4 so I used some of the approach from https://bbolker.github.io/mixedmodels-misc/notes/phylog.html and https://github.com/bbolker/mixedmodels-misc/blob/master/notes/phylocomp.rmd . if your curious the hard part is in the modify_phylo_retrms function.

all that is to say that clutch size explains about 16% of variation in urbanness. Will generalize once #14 is sorted out.

[wcornwell]

added complication, via shinichi, https://www.biorxiv.org/content/early/2017/07/26/144170

[coreytcallaghan]

Hmmm..... I haven't read it fully yet, but seems interesting. What exactly is the complication?

It seems like most of his discussion is focused on using R2 to compare among models or model sets, but we are using R2 to relate the relative effect sizes of parameters?

Could we switch to using the package?

I wouldn't classify what we are trying to do as 'hypothesis testing', so I guess that R2lr wouldn't be the best for us - as the author hints at in the last paragraph.

[wcornwell]

Yeah, if we go back to phylolm instead of lme4 we can use that package.

It's about calculating the percentage variance explained by a random variable, when the random variable has a known covariance (the phylogeny). Shinichi thought his formula would work in that case, but Ives (I think) is politely disagreeing.

I would like to make a bar plot in the end showing the variation in urbanness explained by phylogeny versus traits. And if we sort this out, can also do it for lots of other projects as well...

[coreytcallaghan]

Ah.... okay. It makes a little more sense now.

Have to really dig into it to fully grasp it. But, I see what you're saying.

Then if you know the variation explained by phylogeny, the traits follows from that I reckon.

[wcornwell]

rr2 looks like it's super buggy. but can probably pull some code out and get it to work.

[wcornwell]

seems like we're standardizing predictors and looking at coefficients