[R] R^2 for multilevel models

[Andy Fugard]

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

# Now how about adjusting everyone's intercept towards the group
# 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