How to specify a vector full of means with degrees of freedom for a Lack of Fit F-Test in R

Question

Currently I'm working through Applied Linear Models $5^{th}$ ed - by Kutner, et al. A question I'm working on is asking me to perform an F-Test for Lack of Fit on my linear model. The linear model is a simple linear model of one parameter nothing too troublesome.

To perform the test one has to assess the difference between the full model and the reduced model. At this current junction the authors have stated to take the full model as $\hat{\mu_{j}} = \bar{Y_{j}}$. Specifically the screenshot below says the following:

The reduced model would be the simple linear model:

I have no problem being able to do this manually within R, by computing the necessary values where need be as I've done for other questions. But I'm trying to improve my R skill set and this is where my problem lies.

I have done some reading to other answers related to this and model comparison can be done directly in the anova() function. But I'm having issues stating my full model correctly to be able to leverage the anova() function. I thought about computing a "vector of means" for the subgroups of data (which I display here just for completeness)

But I'm going to run into the problem of the anova() function most likely not being able to compute the degrees of freedom correctly. My data set is very small and this seems like the sort of situation that would show up all the time. With huge data sets I wouldn't see it being feasible to compute things manually so surely there has to be a way for me to phrase my Full Model properly to allow for the computation of means from the subgroups of replicates. But how do I do so? is the question of the day.

Answer 1

First, you should note that for the reduced model (in this context), $\bar{Y}_j$ is the average of $Y$ values where $X=x_j$. That is, we relate to the $X$ variable (in this case, hours) as a category rather than to its actual value. This is simply a one-way ANOVA model. The full model, where we do relate to the nominal value of $X$, is the simple linear regression (just as you said).

It's a bit confusing in R because the anova function is for comparing models, but an ANOVA model is built using the aov function. In any case, this is what you're looking for:

df <- data.frame(
  concentration = c(.07, .09, .08, .16, .17, .21, .49, .58, .53, 1.22, 1.15, 1.07, 2.84, 2.57, 3.1),
  hours = c(9, 9, 9, 7, 7, 7, 5, 5, 5, 3, 3, 3, 1, 1, 1)
)
full <- lm(concentration~hours, data = df)
reduced <- aov(concentration~as.factor(hours), data = df)

the models are

> summary(full)

Call:
lm(formula = concentration ~ hours, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.5333 -0.4043 -0.1373  0.4157  0.8487 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.5753     0.2487  10.354 1.20e-07 ***
hours        -0.3240     0.0433  -7.483 4.61e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4743 on 13 degrees of freedom
Multiple R-squared:  0.8116,    Adjusted R-squared:  0.7971 
F-statistic: 55.99 on 1 and 13 DF,  p-value: 4.611e-06

> summary(reduced)
                 Df Sum Sq Mean Sq F value   Pr(>F)    
as.factor(hours)  4 15.364   3.841     244 6.38e-10 ***
Residuals        10  0.157   0.016                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

and the output is:

> anova(full,reduced)
Analysis of Variance Table

Model 1: concentration ~ hours
Model 2: concentration ~ as.factor(hours)
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1     13 2.9247                                  
2     10 0.1574  3    2.7673 58.603 1.194e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1