Misbehavior of `summary.lm` [duplicate]

Question

I found a misbehavior of the function summary.lm that might be a bug.

Data

d = structure(
  list(
    Treatment = structure(
      c(1L, 1L, 3L, 1L, 2L, 1L, 3L, 2L, 1L, 1L, 2L, 3L, 2L, 2L, 3L, 3L, 3L, 1L, 2L, 3L, 1L, 3L, 2L, 3L, 2L), 
      .Label = c("A", "B", "C"), class = "factor"), 
    Elapsed = c(108L, 110L, 108L, 90L, 100L, 105L, 103L, 120L, 120L, 119L, 45L, 119L, 100L, 80L, 70L, 120L, 112L, 45L, 103L, 85L, 120L, 110L, 110L, 120L, 120L), 
    Correct = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE)),
  row.names = c(NA, 25L), 
  class = "data.frame")

MWE

I perform an OLS regression of Elapsed vs. Treatment and Correct that gives rise to a degenerate model:

mdl = lm(Elapsed ~ Treatment:Correct + 0,data=d)

summary(mdl)

Then I obtain the following summary:

## 
## Call:
## lm(formula = Elapsed ~ Treatment:Correct + 0, data = d)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.571  -8.000   6.429  10.429  20.429 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## TreatmentA:CorrectFALSE  120.000     22.051   5.442 2.99e-05 ***
## TreatmentB:CorrectFALSE  120.000     15.592   7.696 2.97e-07 ***
## TreatmentC:CorrectFALSE  111.000     15.592   7.119 9.06e-07 ***
## TreatmentA:CorrectTRUE    99.571      8.335  11.947 2.79e-10 ***
## TreatmentB:CorrectTRUE    89.667      9.002   9.960 5.61e-09 ***
## TreatmentC:CorrectTRUE   103.571      8.335  12.427 1.43e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.05 on 19 degrees of freedom
## Multiple R-squared:  0.9658, Adjusted R-squared:  0.9549 
## F-statistic: 89.31 on 6 and 19 DF,  p-value: 6.789e-13

The coefficients computed by the lm() are correct. Actually they can be computed, e,g., with with(d,aggregate(Elapsed,list(Treatment,Correct),mean)).

Unfortunately the R-squared value is completely wrong as well as the p-values.

For instance, if I compute the R-squared manually

SSy = sum( (d$Elapsed - mean(d$Elapsed))^2)
SSerr = sum( mdl$residuals^2 )
Rsq = (SSy-SSerr) / SSy

I obtain Rsq = 0.1853978, which is correct if I look at the plot

Any thought on that?

Answer 1

The reason for the difference is that you have + 0 in your formula, setting an intercept of 0. This means it is comparing the model not to one with y equal to mean(Elapsed), but one with y equal to 0.

The equivalent R squared calculation would be:

SSyIntercept = sum( (d$Elapsed)^2)
(SSyIntercept-SSerr) / SSyIntercept
# [1] 0.9657585

If you fit the model with

mdl = lm(Elapsed ~ Treatment:Correct,data=d)

Then you'll find the same R-squared as in your calculation.