Re: statistical comparison of parameters from two applications of NonlinearModelFit

*To*: mathgroup at smc.vnet.net*Subject*: [mg116104] Re: statistical comparison of parameters from two applications of NonlinearModelFit*From*: Ray Koopman <koopman at sfu.ca>*Date*: Wed, 2 Feb 2011 06:08:31 -0500 (EST)*References*: <ii8s8s$2o1$1@smc.vnet.net>

On Feb 1, 3:54 am, dantimatter <goo... at dantimatter.com> wrote: > Hello Everyone, > > First, I must apologize if the following questions are silly; > I am neither a statistician nor a particularly good Mathematica > user (although I do occasionally pretend to be one). And now, > on the question: > > I've got some data that I try to fit with NonlinearModelFit: > > time1 = {0, 100, 200, 300, 400, 500, 600, 700, 800}; > data1 = {1, 0.92, 0.79, 0.81, 0.81, 0.78, 0.72, 0.75, 0.68}; > > nlm1= NonlinearModelFit[Transpose[{time1, data1}], > a1*Exp[-b1*x] + c1, {{a1, .5}, {b1, 0.001}, {c1, 0.5}}, x, > ConfidenceLevel -> .95] > > with results: > > In[11] := nlm1["BestFitParameters"] > nlm1["ParameterConfidenceIntervals"] > nlm1["ParameterErrors"] > > Out[11] = {a1 -> 0.295575, b1 -> 0.00334141, c1 -> 0.696761} > Out[12] = {{0.179807, 0.411344}, {-0.000284765, 0.00696758}, > {0.582507, 0.811014}} > Out[13] = {0.0473121, 0.00148194, 0.0466929} > > I then do the same thing on a second data set: > > time2 = {0, 100, 200, 300, 400, 500}; > data2 = {0.8, 0.73, 0.65, 0.61, 0.58, 0.57}; > > nlm2= NonlinearModelFit[Transpose[{time2, data2}], > a2*Exp[-b2*x] + c2, {{a2, .5}, {b2, 0.001}, {c2, 0.5}}, x, > ConfidenceLevel -> .95] > > In[15]:= nlm2["BestFitParameters"] > nlm2["ParameterConfidenceIntervals"] > nlm2["ParameterErrors"] > > Out[15] = {a2 -> 0.287834, b2 -> 0.00362383, c2 -> 0.516796} > Out[16] = {{0.204254, 0.371413}, {0.00121124, 0.00603641}, > {0.427572, 0.60602}} > Out[17] = {0.0262627, 0.000758091, 0.0280362} > > Based on our understanding of the real system, we actually expect > b2 <= b1, but this not what the fitting tells us (of course, there's > a bit of error in the extracted parameters). > > So what I'd really like to know is: > > (1) how can I make sense of these "ParameterConfidenceIntervals" > and "ParameterErrors" (how are they determined, and what is meant > by them, beyond their brief description in > http://reference.wolfram.com/mathematica/tutorial/StatisticalModelAnalysis.html )? > > and > > (2) is there a way to establish some kind of upper limit on the > detectability of the difference between b2 and b1? > > Any thoughts you might have would be most appreciated. > > Many thanks! The parameter estimates are the values of the parameters (a, b, c) that minimize the sum of squared differences between the observed data and the values predicted from the model. Each ParameterError is an estimate of the standard error of the corresponding parameter estimate. Each ParameterConficenceInterval is an interval that you are (100*ConfidenceLevel)% certain contains the true value of the corresponding parameter *if the model is complete (all the variables that matter are included), the form of the model function is correct, and the distributional assumptions about the errors (normality, homoscedasticity, independence) are correct*. Each confidence level is constructed as estimated_parameter +/- ( estimated_standard_error * Quantile[StudentTDistribution[degrees_of_freedom], (1 + ConfidenceLevel)/2] ), where degrees of freedom = #_of_data_points - #_of_estimated_parameters. Each estimated standard error is the square root of the product of two factors -- one that is common to all the estimates, and one that is specific to each estimate. The common factor is the mean square error of approximation: the sum of squares that the parameter estimates minimize, divided by the degrees of freedom. The specific factors are more complicated, so I'll say only that: 1. They do not directly reflect the size of the data-approximation errors; 2. They are strongly tied to the values of the independent variable (your 'time'); 3. Ceteris paribus, they are inversely proportional to the number of data points. Now we're ready to talk about detecting a difference between b1 and b2. The formal test of the hypothesis that there is no difference is a student t: (b1 - b2)/Sqrt[ASE[b1]^2 + ASE[b2]^2] = -.170 . The degrees of freedom are data dependent, but in your case the t is so small that it will never be significant, no matter what the degrees of freedom are. To achieve significance, the t would have to be at least 12 times its current value, which means you would need 144 times as many sample points as you currently have. (Note: I've cut a lot of corners here. A proper analysis would be much messier but would lead to the same conclusion -- you need a whole lot more data, on the order of 200 times as much as you have shown us.) Finally, here's some code that should reproduce some of the NonlinearModelFit output: r1 = FindFit[Transpose@{time1,data1}, a1*Exp[-b1*x]+c1, {{a1,.5},{b1,.001},{c1,.5}}, x] f1 = a1*Exp[-b1*time1]+c1; e1 = f1 - data1; df1 = Length@data1 - 3; var1 = (e1.e1/.r1)/df1 (* error variance *) cov1 = var1*Inverse[Transpose at #.#]&@Outer[D[##]&,f1,{a1,b1,c1}]/.r1; ase1 = Sqrt@Diagonal@cov1 (* ASEs *) cov1/Outer[Times,ase1,ase1] (* correlations *) ({a1,b1,c1}/.r1) + Outer[Times,ase1, (* CIs *) {-1,1}]*Quantile[StudentTDistribution[df1],.975]

**Re: MatrixPlot**

**Re: Plot artifact of NDSolve result**

**Re: statistical comparison of parameters from two applications of NonlinearModelFit**

**Re: statistical comparison of parameters from two applications of NonlinearModelFit**