# Multiple regression with categorical predictor at two levels in lm()

I have no idea where I came across this code, but I am modeling biological activity of plant tannins in response to an environmental variable (yes,no format) and two continuous variables. The format I was using is:

mod=lm(response~treat/continuousvar1+treat/continuousvar2-1)


My output looks like this, and I think that it is giving me two multiple regressions, one in each environmental treatment (with yes and no rows giving intercepts) and both continuous predictors slopes within that regression. But, if that's the case, what is up with reporting the overall model p value and R^2? If I do this with a single continuous variable, I get back the same slopes and intercepts as when I split the data into yes and no treatments and fit a regression line. The big questions being: is it doing what I think it is, and what can I make of the overall model p?

Call:
lm(formula = stemnit ~ treat/pbo + treat/nitrogenprcnt - 1)

Residuals:
Min        1Q    Median        3Q       Max
-0.199408 -0.024005 -0.003304  0.039079  0.149653

Coefficients:
Estimate Std. Error t value Pr(>|t|)
treatNo                 0.84716    0.37564   2.255   0.0478 *
treatYes                0.61179    0.50882   1.202   0.2569
treatNo:pbo             0.06525    0.53332   0.122   0.9050
treatYes:pbo            0.39890    0.42716   0.934   0.3724
treatNo:nitrogenprcnt   0.01604    0.33802   0.047   0.9631
treatYes:nitrogenprcnt  0.22960    0.46016   0.499   0.6286

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1123 on 10 degrees of freedom
Multiple R-squared:  0.9909,    Adjusted R-squared:  0.9854
F-statistic:   181 on 6 and 10 DF,  p-value: 1.309e-09

• A separate issue might be lurking here. The high adjusted $R^2$ and low overall p-value are startling in light of the nonsignificant coefficients: that suggests there is a severe collinearity problem. (Have you run collinearity diagnostics?) It also appears you are attempting to fit a six-parameter model to just 16 observations, which would be many more than so few observations could reliably support.
– whuber
Aug 4, 2016 at 22:27
• The weirdly high r squared and low p was resolved by the first answer (it was a coding issue). I am very confident in these data, because there area actually 10 individuals averaged out in each of the 16 experimental units. Aug 5, 2016 at 0:35

The estimated parameters (plus standard errors and test statistics) as well as the fitted values and residuals all do what you expect them to do. But the $R^2$ and $p$-value do not. The reason for this is subtle but is, in fact, indicated in ?summary.lm:

r.squared: R^2, the 'fraction of variance explained by the model',

R^2 = 1 - Sum(R[i]^2) / Sum((y[i]- y*)^2),

where y* is the mean of y[i] if there is an intercept and
zero otherwise.


Note the last sentence. If the lm object was fitted with an intercept, then the reference model against which the residual sum of squares is compared is y ~ 1 otherwise it is y ~ 0. For regression lines that are actually forced through the origin this makes sense but for models with no overall intercept but group-specific intercepts (like yours) I personally find this convention confusing. But I guess in general it is not trivial to check whether the constant regressor is in $\mathrm{Im}(X)$ of the regressor matrix.

In your case I would simply use the formula stemnit ~ treat/pbo + treat/nitrogenprcnt without the -1. Then you will get the $R^2$ etc. you expect while having the same fitted values etc. and the nested (group-specific) coding of the slopes. Only the intercepts will have a treatment coding rather a nested coding, i.e., (Intercept) = 0.84716 and treatYes = -0.23537).

• That fixed it! I thought that -1 was required for to gain separate intercepts within a treatment, but I now recall learning that in lmer or nlme. Aug 5, 2016 at 0:31