On obtaining PearsonChiSquareTest of the NonLinearModelFit

Question

Using NonlinearModelFit I found the FittedModel of the data twice: once by linearly parametrizing and then by nonlinearly parametrizing as following,

curve1 = NonlinearModelFit[data, c + a*x + b*y, {a,b,c}, {x,y}, MaxIterations->2000]
Curve2 = NonlinearModelFit[data, c + a*x^n + b*y^m, {a,b,c,n,m}, {x,y}, MaxIterations->2000]

Now, I want to perform PearsonChiSquareTest on both of these models to know the best description of the data. Commands like,

PearsonChiSquareTest[data,curve1]
PearsonChiSquareTest[data,curve2]

for example, are not working.

I request you to help me in this regard. Thanks in advance..

Answer 1

Pearson-$\chi^2$ test would show if the data comes from a certain distribution (normal by default). For a test that measures goodness of fit of a linear model you could use an F-test that would compare the model with the baseline (would be an intercept only in the linear case). For non-linear models, you should use other criteria, such as AIC or BIC. For example, if you have the data:

gp[x_, y_] := 
 3.1 - 3. x^1.2 + 2. y^2 + RandomVariate[NormalDistribution[0, 1000]]
data = Flatten[Table[{x, y, gp[x, y]}, {x, 1, 100}, {y, 1, 100}], 1];
ListPointPlot3D[data]

Fitting the data now to both models you suggested:

curve1 = NonlinearModelFit[data, c + a*x + b*y, {a, b, c}, {x, y}]
curve2 = NonlinearModelFit[data, c + a*x^n + b*y^m, {a, b, c, n, m}, {x, y}]

Row[{
  Show[Plot3D[curve1[x, y], {x, 1, 100}, {y, 1, 100}, 
    ImageSize -> Large], Graphics3D[{Blue, Point[data]}]],
  Show[Plot3D[curve2[x, y], {x, 1, 100}, {y, 1, 100}, 
    ImageSize -> Large], Graphics3D[{Blue, Point[data]}]]
  }]

Visually, we see that of course, curve2 fits better, but we want to be able to formally say it. We can always look at the properties of fitted models, say curve2:

curve2["Properties"]

We can then construct a helper function to have a nice table of the tests we want:

gf[f_, tests_] := f[#] & /@ tests
tests = {"AdjustedRSquared", "AIC", "AICc", "BIC", "RSquared"};
Grid[Prepend[Transpose@{tests, gf[curve1, tests], gf[curve2, tests]}, 
  Style[#, Bold] & /@ {"Test", "Model 1", "Model 2"}], 
 Alignment -> Left, Dividers -> {None, 2 -> True}]

Clearly, you can now say, that regardless of the criteria, the curve2 is the way to go (R-squared is higher, and both AIC and BIC are lower). Now, this won't always be the case if the models are comparable, so pay attention to that.

Now, the other thing to check is the residuals. They have to have expectation of 0, be (ideally) normally distributed, and clearly don't have any patterns in the plot that would indicate bias.

res1 = curve1["StandardizedResiduals"];
res2 = curve2["StandardizedResiduals"];
Row[{
  ListPlot[res1, PlotRange -> All, Filling -> Axis, 
   ImageSize -> Medium],
  ListPlot[res2, PlotRange -> All, Filling -> Axis, 
   ImageSize -> Medium]
  }]

We can now check the expectation of the residuals:

{Expectation[x, x \[Distributed] EmpiricalDistribution[res1]],
Expectation[x, x \[Distributed] EmpiricalDistribution[res2]]}

{0.0000368449, -4.9723*10^-7}

Obviously, the residuals of the first model are biased, while the second one looks normal-ish. Let's verify it using Pearson-$\chi^2$ test now:

h1 = PearsonChiSquareTest[res1, Automatic, "HypothesisTestData"];
h2 = PearsonChiSquareTest[res2, Automatic, "HypothesisTestData"];
{h1["TestDataTable"], h2["TestDataTable"]}

The $H_0$ in the test is that the data comes from a normal distribution, so low $p-value$ means that we can reject it. This is the clear case for the curve1, while we cannot reject this hypothesis for curve2, which means it is an unbiased estimator.