# How to do a lack-of-fit test by Mathematica

I want to do a lack-of-fit test to my linear fit result by Mathematica. I know the error sum of squares is $$\mathrm{SS_E}=\sum_{i=0}^n(y_i-(b_0+b_1x_i))^2$$ and the pure sum of square is $$\mathrm{SS_{E,pe}}=\sum_{i=1}^k\sum_{j=1}^{n_i}(Y_{ij}-\overline Y_{i})^2$$ , where we assume we have $$n_1+\cdots+n_k$$ sets of data and among them, we have $$x$$ values $$x_1,\ldots,x_k$$.

I've used LinearModelFit in Mathematica to do a linear regression and used ANOVATable to obtained the total sum of square and $$\mathrm{SS_E}$$.

However, how I can get pure sum of squares or how I can do a lack-of-fit test directly?

Edit: I’ve added a little bit about the statistical meaning of “lack-of-fit” and “pure” sum of squares.

When fitting a model deviations can occur (among other things) because of measurement error (labeled “pure error”) and because the model is inadequate (labeled “lack-of-fit error”). A model might be $$y=f(x)+\epsilon+\gamma$$ where $$y$$ is the response, $$f(x)$$ is some function with constants to be estimated, $$\epsilon$$ has a normal distribution with variance $$\sigma^2$$, and $$\gamma$$ is a deviation from the simple linear model. If the linear model is correct, then $$\gamma=0$$.

If there are multiple responses with the same predictor(s), then we can get two independent estimates of $$\sigma^2$$. If $$\gamma=0$$, then the ratio of the estimates will vary around 1 according to an $$F$$ distribution. If $$\gamma \neq 0$$, then the ratio should tend to be larger than 1.

The example below has multiple observations for just one of the 10 “reps”. To obtain the “lack-of-fit” and “pure error” sums of squares and associated degrees of freedom which are used to construct the F-test are outlined below. One run of the model (which might be LinearModelFit or NonlinearModelFit) obtains the pure error statistics and another run of the model obtains the lack-of-fit statistics.

The essential need is to have multiple observations with the same predictors. For more details see https://en.wikipedia.org/wiki/Lack-of-fit_sum_of_squares.

End of edit

I certainly hope that there is a better way than the following brute force approach.

The example data is taken from a SAS example so that results can be compared. I have added a "nominal" variable to the beginning of each observation to show the particular replication group (called rep below) to which the observation belongs.

data = {{1, 10.0, 1.0, 100, 140, 6.0, 37.5}, {2, 10.0, 1.0, 120, 180, 3.0, 28.5},
{3, 10.0, 2.0, 100, 180, 3.0, 40.4}, {4, 10.0, 2.0, 120, 140, 6.0, 48.2},
{5, 15.0, 1.0, 100, 180, 6.0, 50.7}, {6, 15.0, 1.0, 120, 140, 3.0, 28.9},
{7, 15.0, 2.0, 100, 140, 3.0, 43.5}, {8, 15.0, 2.0, 120, 180, 6.0, 64.5},
{9, 12.5, 1.5, 110, 160, 4.5, 39.0}, {9, 12.5, 1.5, 110, 160, 4.5, 40.3},
{9, 12.5, 1.5, 110, 160, 4.5, 38.7}, {9, 12.5, 1.5, 110, 160, 4.5, 39.7}};


Using your notation $$k=9$$ and $$n_1, n_2, n_3, n_4, n_5, n_6, n_7, n_8, n_9$$ equals '{1,1,1,1,1,1,1,1,4}.

LinearModelFit is run twice: once to get the pure error sum of squares ($$SS_{Epe}$$) and once to get the total error sum of squares ($$SS_E$$).

lm = LinearModelFit[data, {rep, x1, x2, x3, x4, x5}, {rep, x1, x2, x3, x4, x5},
NominalVariables -> rep];
lm["ANOVATable"] /. "Error" -> "Pure Error"


lm2 = LinearModelFit[data, {x1, x2, x3, x4, x5}, {rep, x1, x2, x3, x4, x5}];
lm2["ANOVATable"]


(* Get the associated sum of squares and degrees of freedom *)
(* Pure Error and degrees of freedom *)
lmANOVA = lm["ANOVATableEntries"]
{dfP, ssP} = lmANOVA[[2, {1, 2}]]

(* Error sum of squares and degrees of freedom *)
lm2ANOVA = lm2["ANOVATableEntries"]
{dfE, ssE} = lm2ANOVA[[Length[lm2ANOVA] - 1, {1, 2}]]

(* F-ratio for lack of fit test *)
f = ((ssE - ssP)/(dfE - dfP))/(ssP/dfP)
(* 22.0689 *)
pValue = 1 - CDF[FRatioDistribution[dfE - dfP, dfP], f]
(* 0.0151209 *)
`

The values are (essentially) identical to the values given in the SAS example.

• I'm not sure about the principle behind it but it works. Thank you very much. Commented Apr 26, 2020 at 9:16