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?


1 Answer 1


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"

ANOVA table with pure error sum of squares

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

ANOVA with total error sum of squares

(* 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.

  • $\begingroup$ I'm not sure about the principle behind it but it works. Thank you very much. $\endgroup$
    – citrate
    Apr 26, 2020 at 9:16
  • $\begingroup$ Good answer, shorter below. $\endgroup$
    – Carl
    Aug 12, 2023 at 9:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.