Re: statistical comparison of parameters from two applications of NonlinearModelFit
- To: mathgroup at smc.vnet.net
- Subject: [mg116139] Re: statistical comparison of parameters from two applications of NonlinearModelFit
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Thu, 3 Feb 2011 05:30:12 -0500 (EST)
For nonlinear models, it seems rash beyond words, to me, to assume
"normality, homoscedasticity, and independence". Hence, I doubt most of
the statistics returned mean much of anything.
The same is true in linear modeling, but you have a better chance of
testing the assumptions.
Bobby
On Wed, 02 Feb 2011 05:08:31 -0600, Ray Koopman <koopman at sfu.ca> wrote:
> 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]
>
--
DrMajorBob at yahoo.com