This is an example of non-linear regression on the Niederkorn’s (1872) observations on 114 bodies.

This dataset provides the main reference database for the development of rigor mortis. I include an animation which shows how one might guess some good starting parameters, and then turn things over to the non-linear regression routine.

Andy Long
Spring, 2018

Clear["Global`*"] (* This removes any definitions that may be in place already; like if you've already been doing some Mathematica, and talking about other functions named f, x-values named a and b, etc. *)

Here is the data. It’s available from
http://www.zoology.ubc.ca/~bio300/homework/hw1_answers.pdf
(although that’s not where I originally got it).

bodies={2,14,31,14,20,11,7,4,7,1,1,2} ;
(* bodies transitioning to rigor mortis at the associated hour: *)
hours=Table[i,{i,2,13}]; (* hours 2 through 13 *)

We compute the cumulative, and will model that -- that is, the number of bodies in rigor mortis by a given hour:

cumulative=Accumulate[bodies];
data=Table[{hours[[i]],cumulative[[i]]},{i,1,12}];
p1=ListPlot[data,PlotRange->{{0,15},{0,120}}]

Now we’ll get our good starting guess -- which we’ll need for non-linear regression -- in two ways: by playing with the parameters of our model, and then by a sort of pseudo linearization. 

Once we have it, we’ll use Mathematica’s native non-linear regression routine to find the parameters.

(* I'm suggesting a model of the form a Exp[-b/t^p] for this data -- I think that it will be nice and flat near the origin, and it will asymptote as I want it to do. This animation allows me to play with the parameters, in order to get a good starting guess. *)
Clear[a,b,c]
Manipulate[
Module[
{},
p2=Plot[a Exp[-b/t^p],{t,0,15}];
Show[p1,p2]
]
,{{a,120},100,150,1}
,{{b,10},0,50,.01}
,{{p,2},.5,3,.01}
]

Let’s try Mathematica’s native non-linear regression routine to model this data, just for fun. Mathematica has its own routine for carrying out non-linear regression. If we use it without a starting guess for the parameters, it gives junk; so we seed it with some values that we obtained by using a transformation first.

(* this is kind of funny: if we don't give NonlinearModelFit any guidance, check out the "fit"! *)
nlm=NonlinearModelFit[data,a Exp[-(b/x^c)],{a,b,c},x]
Show[p1,Plot[nlm[x],{x,0,14}],Frame->True]

Let’s make a good guess:

I’m going to try to fit a model of the form 
y=aExp[-(b/x^c)]
to the cumulative data (don’t ask me where I got that!), which I want to turn into a probability (so my graph will asymptote to 1). Let me guess that the asymptote is around a=120.
 
We divide the cumulative by 120 to start as a preliminary guess (slightly more than the total number of bodies), which we might guess is the asymptote from the ListPlot above of the data,
y/120=Exp[-(b/x^c)]
because we don’t believe that we’ve hit “1” by the last body. There’s a body somewhere that would have taken longer to go into rigor mortis.

Then we log-transform the data, so we get 
Log[y/120]= -b/x^c
-Log[y/120]= b/x^c
and then we log-transform again,
Log[-Log[y/120]]=Log[ b]-cLog[x]
So I can get rough estimates for b and c by doing linear regression on 
Log[-Log[y/120]] as a function of Log[x] , and then use a=120 as well. 

Good lord, that’s a little ugly. Sorry for all the horrible wickedness....

It is a good guess, as you will see:

starta=120;
lm=LinearModelFit[Transpose[{Log[hours],Log[-Log[cumulative/starta]]}],x,x]
{startb,startc}=lm["BestFitParameters"]
(* actually I have to transform these back: 
	My constant is Log[b], and the slope is -c,
so
*)
startb=Exp[startb]
startc=-startc
(* What did I tell you? It's a good guess! *)
Show[Plot[starta Exp[-startb/t^startc],{t,0,15},Frame->True], p1]

So here we go: we’re going to use these start values to help NonlinearModelFit do a better job -- and I think that you’ll agree that it does!

nlm=NonlinearModelFit[data,a Exp[-(b/x^c)],{{a,starta},{b,startb},{c,startc}},x]
{pa,pb,pc}=nlm["BestFitParameters"]
a=a/.pa
Normal[%]
(* Now isn't this a little better?:*)
Show[Plot[nlm[x],{x,0,14},Frame->True],ListPlot[data]]

Plot[nlm'[x]/a,{x,0,14}] (* here's the PDF -- Probability Density Function : *)

Now we’ll take a look at some diagnostics:

nlm["ParameterTable"]
ListPlot[nlm["FitResiduals"]]

yhat=Table[nlm[i], {i,hours}];
ybar=Mean[cumulative];
residuals=cumulative-yhat;
(* Notice that my rsquared calculation doesn't quite match up with theirs...*)
rsquared=1-residuals.residuals/((cumulative-ybar).(cumulative-ybar))
nlm["RSquared"]
n=Length[hours];
p=3; (* number of parameters estimated *)
rsquaredadj=1-(n-1)/(n-p) (1-rsquared)
(* See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2892436/
for a discussion of the mis-application of r-squared in non-linear models. They argue for the use of no r-squared. *)

Now we’ll do it by linear regression! But we’ll have to do it iteratively, by successively better approximation.

Our objective will be to match the results of Mathematica’s routine, NonlinearModelFit, including estimates and t-values.

Clear[f,fa,fb,fc,a,b,c]
(* This is our model, which I write as a multivariate function: *)
f[x_,a_,b_,c_]:=a Exp[-(b/x^c)]
(* We'll need the partial derivatives of f: *)
fa[x_,a_,b_,c_]=D[f[x,a,b,c],a]
fb[x_,a_,b_,c_]=D[f[x,a,b,c],b]
fc[x_,a_,b_,c_]=D[f[x,a,b,c],c]

Okay, now we iterate:

runs=6;
alpha=starta;
beta=startb;
gamma=startc;
Do[
estimates=f[hours,alpha,beta,gamma];
fpa=fa[hours,alpha,beta,gamma];
fpb=fb[hours,alpha,beta,gamma];
fpc=fc[hours,alpha,beta,gamma];
(* The IncludeConstantBasis->False command means that we don't want an intercept term: *)
diffLinMod=LinearModelFit[Transpose[{fpa,fpb,fpc,(cumulative-estimates)}],{x,y,z},{x,y,z},IncludeConstantBasis->False];
{newa,newb,newc}=diffLinMod["BestFitParameters"];
alpha=alpha+newa;
beta=beta+newb;
gamma=gamma+newc;
Print[diffLinMod["ParameterTable"]],
{i,runs}
]
alpha
beta
gamma

Notice that the parameter estimates (alpha, beta, gamma) and their standard errors (from the table) are exactly the same as those given by NonlinearModelFit.

From these we can compute the confidence intervals, and the t-statistics. And once we have the model, we can look at the estimates (yhats), as well as the residuals, and compute R^2 values.

Look at those P-values! Isn’t that funny? That we’re successful when the P-values are 1, and the estimates are so close to zero. Those estimates are the differences between the “true” parameters and our “guess” parameters -- so we want them to go to zero!