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