# [R] R^2 for multilevel models

Andy Fugard a.fugard at ed.ac.uk
Mon Aug 13 22:06:49 CEST 2007

```Hi there,

In multiple regression one way to view R^2 is as (the square of) the
correlation between original y's and the estimated y's.

Suppose you fit a multilevel model with random intercept for each
cluster.  Would it be valid to compute an R^2 by using fixed effects
plus the group intercepts to reduce the residuals?

I suspect this has been done and, given its absence from the lmer
output, refuted somewhere.  A reference would be lovely.

Code below.

Many thanks,

Andy

#################################################################

require(lme4)

# First generate some data

people = 100
obs = 4

# The random intercepts:

sub_noise = data.frame(id = 1:people, rand_int = rnorm(people,0,60))

# Merge everything

id = rep(1:people,obs)
thedata = data.frame(id)
thedata = merge(sub_noise, thedata)
thedata\$x = rep(1:obs,people)

# Setup the relationship between x and y

thedata\$y = 23*thedata\$x + 20 + thedata\$rand_int
+ rnorm(people*obs, 0, 20) # plus residuals

# Have a look

plot(y~x, data = thedata)

# Now fit a standard regression model

lm1 = lm(y ~ x, data = thedata)
summary(lm1)

# Use the model to get the R^2
predicted = coef(lm1)[1] + thedata\$x * coef(lm1)[2]
plot(thedata\$y, predicted)

# It's a noisy mess

cor(thedata\$y, predicted)^2

# intercept?

lmer1 = lmer(y ~ x + (1|id), data=thedata)
summary(lmer1)

# Get the random intercepts and stick them in a table

ran_effects = data.frame(rownames(ranef(lmer1)\$id), ranef(lmer1)\$id[1])
names(ran_effects) = c("id", "b")
ran_effects_data = merge(thedata, ran_effects)

# Now compute

predicted.ml = fixef(lmer1)[1] + ran_effects_data\$x * fixef(lmer1)[2]
+ ran_effects_data\$b
plot(thedata\$y, predicted.ml)
cor(thedata\$y, predicted.ml)^2

# Looks much nicer

```