Test if NonlinearModelFit parameters are different between two data sets

Question

I'm fitting the following model to different data sets using NonlinearModelFit:

model[x_]:=-Log10[A[x] + B[x]]
A[x_] := (A0)*Exp[(za)*x/25.7]
B[x_] := (B0)*Exp[-(zb)*x/25.7]
params = {A0, B0, za, zb};
cons = Thread[params > 0];
nlm=NonlinearModelFit[data, {model[x], cons}, params, x]

Using the data in the attached files I get,

For nlmfitdata1:

    Estimate Standard Error t-Statistic P-Value
A0  17.1721     3.32073     5.1712      7.64783*10^-7
B0  0.139581    0.585399    0.238437    0.811881
za  0.26739     0.091721    2.91525     0.00412296
zb  0.259632    0.390068    0.665607    0.506727

For nlmfitdata2:

   Estimate  Standard Error t-Statistic P-Value
A0  8.19873     1.3118      6.24999     3.37415*10^-9
B0  0.0750962   0.233958    0.320982    0.74863
za  0.243825    0.0783766   3.11095     0.00219787
zb  0.291387    0.292965    0.994615    0.32138

Now I want to know if there a statistical significant difference between my parameters estimates for both data sets. Is there any test in NonlinearModelFit that can address this?

Yo can get the data here: nlmfitdata1 nlmfitdata2

Answer 1

I don't think that NonlinearModelFit has any ability to do what you are looking to do. I am not the foremost statistician in the world, but I think you should be able to get an estimate of the independence of these parameters directly using the mean parameter estimates and their associated standard errors. In this case, I would suggest using a Welch's t-test to test for significant differences between the parameter estimates. For each parameter you need to calculate the t-statistic and the pooled degrees of freedom.

Start by building each model:

model[x_] := -Log10[A[x] + B[x]]
A[x_] := (A0)*Exp[(za)*x/25.7]
B[x_] := (B0)*Exp[-(zb)*x/25.7]
params = {A0, B0, za, zb};
cons = Thread[params > 0];
nlm1 = NonlinearModelFit[data1, {model[x], cons}, params, x]
nlm2 = NonlinearModelFit[data2, {model[x], cons}, params, x]

Now, extract the parameter estimates and their standard errors:

nlm1parms = Quiet@nlm1["ParameterConfidenceIntervalTableEntries"];
nlm2parms = Quiet@nlm2["ParameterConfidenceIntervalTableEntries"];

Next, calculate the t-statistics and the dofs:

tstats = Table[(nlm1parms[[i, 1]] - nlm2parms[[i, 1]])/Sqrt[
  nlm1parms[[i, 2]]^2 + nlm2parms[[i, 2]]^2], {i, 1, 
   Length@nlm1parms}]

dofs = Table[(nlm1parms[[i, 2]]^2 + nlm2parms[[i, 2]]^2)^2/(
  nlm1parms[[i, 2]]^4/(Length@data1 - Length@params) + 
   nlm2parms[[i, 2]]^4/(Length@data2 - Length@params)), {i, 1, 
   Length@nlm1parms}]

With this information, you can directly compute the p-values corresponding to the hypothesis tests related to the presence of differences between the model parameter estimates:

pvals = Table[
  2*(1 - CDF[StudentTDistribution[dofs[[i]]], Abs@tstats[[i]]]), {i, 
   1, Length@nlm1parms}]

(*{0.0128006, 0.918636, 0.845281, 0.948147}*)

Based on this analysis, it looks like your data only support the conclusion that the first parameter, A0, in the models is significantly different.

Answer 2

The first part of this answer is nearly identical to that of @Marchi. The second part is not.

When the model appropriately fits the data (meaning the model assumptions are not wildly violated by features of the data), the follow code could be used:

nlm1 = NonlinearModelFit[data1, {model[x], cons}, params, x];
nlm2 = NonlinearModelFit[data2, {model[x], cons}, params, x];
mle1 = params /. nlm1["BestFitParameters"];
mle2 = params /. nlm2["BestFitParameters"];
var1 = Diagonal[nlm1["CovarianceMatrix"]];
var2 = Diagonal[nlm2["CovarianceMatrix"]];
z = (mle1 - mle2)/(var1 + var2)^0.5;
pvalues = 2 (1 - CDF[NormalDistribution[0, 1], Abs[z]]);
TableForm[Transpose[{mle1, mle2, var1^0.5, var2^0.5, z, pvalues}],
 TableHeadings -> {{"A0", "B0", "za", "zb"},
   {"MLE1", "MLE2", "Std.Err1", "Std.Err2", "z-value", "P-value"}}]

A normal distribution is used rather than a Student's t-distribution because the associated degrees of freedom is way over 100.

But all of this is contingent on the models providing an appropriate fit: a good fit is not required but rather an appropriate fit that accounts for the structure of the data is required. A special structure in the data is evident from a plot of the two datasets:

GraphicsRow[{ListPlot[data1, PlotLabel -> "Data 1", PlotRange -> {All, {-1.4, 0.1}}],
  ListPlot[data2, PlotLabel -> "Data 2", PlotRange -> {All, {-1.4, 0.1}}]}]

We see that almost certainly that there were subjects/experimental units which had repeated measures which is not the error structure assumed by NonlinearModelFit. This is more evident if the data is a bit restructured to highlight the consecutive measurements (but no changes in the ordering of the data nor any changes in the values of the data):

GraphicsRow[{ListPlot[data1, Joined -> True, PlotLabel -> "Data 1", PlotRange -> {All, {-1.4, 0.1}}],
  ListPlot[data2, Joined -> True, PlotLabel -> "Data 2", PlotRange -> {All, {-1.4, 0.1}}]}]

An appropriate analysis needs to account for such repeated measures on various experimental units. Such regression models are called "Mixed Models" where there is more than one random term. As far as I know Mathematica does not offer any direct way of fitting such models. (Although I wish they would add such models.) One such model to fit such a structure is called a "Random coefficients model" where each subject/experimental unit has a random deviation from one or more of the coefficients.

In practice the maximum likelihood estimates of the coefficients are likely not to differ by very much even after accounting for the repeated measures. However, the associated standard errors might likely be very, very different and that affects the ability to make an appropriate assessment. So the above standard errors and P-values should be considered bogus.