Any suggestions for improving my fitting?

Question

This is not a completely Mathematica question, but I hope people who are familiar both with Mathematica and statistics could help me on this.

I have a set of data as

data = {{0, 0.046}, {40, 0.111}, {80, 0.291}, {120, 0.808}, {160, 1.742}, {200, 3.319}, {240, 5.017}, {280, 5.503}, {320, 5.897}}

The data seems to follow a logistic curve, thus for fitting, we can write:

F = NonlinearModelFit[data, A/(1 + E^(-(t - t0)/b)), {{A, 1}, {t0, 50}, {b, 10}}, t];
F["ParameterTable"]

Here is the result:

which is not bad, but I would like to improve it. Any suggestions here?

One thing the came to my head is to use inverse variances as weights. Standard deviations are given as

sd = {0.003, 0.012, 0.023, 0.056, 0.083, 0.216, 0.526, 0.366, 0.313}

thus, we have

F = NonlinearModelFit[data, A/(1 + E^(-(t - t0)/b)), {{A, 1}, {t0, 50}, {b, 10}}, t, Weights -> 1/sd^2];
F["ParameterTable"]

still, the results are not very satisfactory. Is there any other way to improve my fitting?

Answer 1

In your comments you mention that the responses are means of 3 values. If the response you're really interested in is the mean of 3 values, then we're ready to go. If the response should be the individual values and you really have a repeated measures design, then you should use something other than Mathematica. The functions LinearModelFit and NonlinearModelFit are just fine for simple models (with "simple" meaning just one error term).

To fit $k$ just replace $1/b$ with $k$. The underlying model that NonlinearModelFit thinks it's estimating is for the $i$-th observation

$$y_i=\frac{a}{e^{-k (t_i-t_0)}+1}+\epsilon_i$$

where $y_i$ is the response, $a$, $k$, and $t_0$ are parameters to be estimated, $t$ is the predictor variable and $\epsilon_i$ has a normal distribution with mean 0 and variance $\sigma^2$.

nlm1 = NonlinearModelFit[data, a/(1 + E^(-k (t - t0))), {{a, 1}, {t0, 50}, {k, 0.1}}, t];
nlm1["ParameterConfidenceIntervalTable"]

We see that your desired 0.30 value for $k$ is included in the 95% confidence interval.

A plot of the predicted response vs the residuals suggests that the error term is incorrectly structured in that the variability seems to increase with an increase in predicted response.

ListPlot[Transpose[{nlm1["PredictedResponse"], nlm1["FitResiduals"]}]]

You also considered taking the logs of the responses.

data2 = data;
data2[[All, 2]] = Log[data2[[All, 2]]];
nlm2 = NonlinearModelFit[data2, Log[a/(1 + E^(-k (t - t0)))], {{a, 1}, {t0, 50}, {k, 0.1}}, t];
nlm2["ParameterConfidenceIntervalTable"]

Note that your desired 0.30+ value for $k$ is not within a much tighter 95% confidence interval.

The residual plot doesn't display the increase in variability as before which suggests meeting the model assumptions better than the other model. In other words, this model provides a better fit to the data. It just doesn't provide a fit that coincides with your desires.

ListPlot[Transpose[{nlm2["PredictedResponse"], nlm2["FitResiduals"]}]]

The data should really be analyzed according to how it was collected. And if the 3 values that constitute the means are available, that's what should be used in an analysis. But while one can use Mathematica to do so, that would entail a bit of programming as opposed to using NonlinearModelFit.

Answer 2

If you know the error in the data is essentially zero, you are better off with Monotonic-Interpolation in a case like this. The whole point of Interpolation is that the function you get goes exactly through each data point.

data=data={{0,0.046},{40,0.111},{80,0.291},{120,0.808},{160,1.742},
{200,3.319},{240,5.017},{280,5.503},{320,5.897}};
f=ResourceFunction["MonotonicInterpolation"][data];
Plot[f[t],{t,0,320},Epilog->{AbsolutePointSize@4,Point@data}]

If you want a closed form for the above function, use the next line.

pw[t_] = ResourceFunction["InterpolatingFunctionToPiecewise"][f, t]

For more on the functions used above see documentation here and here.

Similar to MonotonicInterpolation there are ResourceFunctions here, here, and here that run faster on large datasets, but MonotonicInterpolation used above gives a smoother function.

***** UPDATE *****

Next I show how to find the maximum of the derivative of f[t] above.

Plot[f'[t],{t,0,320}]

Use a smaller PlotRange and it is clear the maximum of f'[t] is between 200 and 213. Look at the Piecewise function pw[t] defined above, and we see the interval from 200 to 213 is one piece of pw[t] and Simplify gives us that piece below.

poly=Simplify[pw[t],200<t<213]

Now we can find the maximum of the derivative of the above polynomial.

FindMaximum[D[poly,t],t]

So in the plot of f'[t] above, the maximum derivative is at {t,f'[t]}={0.04799,207.91}