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

Question

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

Answer 1

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).