summary[["r.squared"]] gives strange results
Felix Flory wrote:
I am simulating an ANOVA model and get a strange behavior from the
summary function. To be more specific: please run the following code
and see for yourself: the summary()[["r.squared"]] values of two
identical models are quite different!!
## 3 x 3 ANOVA of two factors x and z on outcome y
s.size <- 300 # the sample size
p.z <- c(0.25, 0.5, 0.25) # the probabilities of factor z
## probabilities of x given z
p.x.z <- matrix(c(40/60, 20/120, 6/60,
14/60, 80/120, 14/60,
6/60, 20/120, 40/60), 3, 3, byrow = TRUE)
## the regression coefficients according to the model.matrix X, that
## is computed later
beta <- c(140, -60, -50, -80, 120, 100, -40, 60, 50)/40
## the factor z and the factor x (in dependence of z)
z <- x <- vector(mode = "integer", length = s.size)
for(j in 1:s.size) {
z[j] <- sample(1:3, 1, prob = p.z)
x[j] <- sample(1:3, 1, prob = p.x.z[, z[j]])
}
## constructing the model.matrix X
zf <- factor(z)
contrasts(zf) <- contr.treatment(nlevels(zf), base = 3)
zm <- model.matrix(~ zf)
xf <- factor(x)
contrasts(xf) <- contr.treatment(nlevels(xf), base = 3)
xm <- model.matrix(~ xf)
X <- cbind(zm, zm * xm[,2], zm * xm[,3])
## the outcome y
y <- X %*% beta + rnorm(s.size, 0, 4)
## the two linear models
lm.v1 <- lm(y ~ X -1)
lm.v2 <- lm(y ~ zf * xf)
## which are equivalent
anova(lm.v1, lm.v2)
print(sd(model.matrix(lm.v1) %*% coef(lm.v1))^2 / sd(y)^2) -
print(sd(model.matrix(lm.v2) %*% coef(lm.v2))^2 / sd(y)^2)
## so far everything is fine but why are the following r.squared
## values quite different?
print(summary(lm.v1)[["r.squared"]]) -
print(summary(lm.v2)[["r.squared"]])
Hi, Felix,
The first model fits your data without an intercept and thus has a
different formula of R^2. The justification should be in any intro
regression text. Here is the relevant snippet from summary.lm:
mss <- if (attr(z$terms, "intercept"))
sum((f - mean(f))^2)
else sum(f^2)
rss <- sum(r^2)
<snip>
ans$r.squared <- mss/(mss + rss)
Try the following to see a direct comparison of the two methods:
lm.v1 <- lm(y ~ X - 1)
lm.v2 <- lm(y ~ zf * xf)
lm.v3 <- lm(y ~ X[, -1])
summary(lm.v1)$r.squared
summary(lm.v2)$r.squared
summary(lm.v3)$r.squared
HTH,
--sundar