Non-linear regression on the Niderkorn’s (1872) observations on times to rigor mortis in 114 bodies. 

This dataset still provides the main reference database for the development of rigor mortis, from what I can determine.

Andy Long
Math Modeling

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}];
dataplot= ListPlot[data,PlotRange->{{0,15},{0,120}},
AxesLabel->{"Hour","Number of Bodies"},
PlotLabel->"Number of bodies in rigor mortis by a given hour"]

We’ll start by using Mathematica’s native non-linear regression routine to model this data.

I’m going to try to fit a model of the form 
				y=a(e^(-(b/(x^c))))     
to the cumulative data (don’t ask me where I got that!), which I want to turn into the probability of a body being in rigor mortis (so my graph will asymptote to 1).

Mathematica has its own routine for carrying out non-linear regression. 

If we use it without a starting guess for the parameters, it uses whatever default values it’s programmed with, and it may give junk.

(* 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[dataplot,Plot[nlm[x],{x,0,14}]]

But we need a better guess, because Mathematica’s default guesses don’t seem to be working.... 

So it’s better if we can seed it with some good guess values. We can obtain some by using a transformation first. Let's see if we can linearize this one (and the answer is "sort of"...). 

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)=(e^(-(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). That is, we're guessing that  a=120 or so.

Then we log-transform the data, so we get 

			Log[(y/120)]= -(b/(x^c))
			
Now move the negative to the other side (and everything on the right is positive). 

			-Log[(y/120)]= (b/(x^c))

This is a power model, and we know how to linearize power models; and so we log-transform again,

		Log[-Log[(y/120)]]=Log[ b]-cLog[x]

Upshot: we can get rough estimates (initial guesses) for b and c by doing linear regression on Log[-Log[(y/120)]] as a function of Log[x] , with an initial guess of a=120 as well. 

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

Here we go: linear regression on a “pseudo-linearized” problem will give us our starting guesses....

starta=120;
lm=LinearModelFit[Transpose[{Log[hours],Log[-Log[cumulative/starta]]}],x,x]
{startb,startc}=lm["BestFitParameters"]

actually I have to transform these back: the constant is Log[b], and the slope is -c, so

startb=Exp[startb]
startc=-startc

Back to our non-linear regression: we’re going to use these start values to help NonlinearModelFit do a better job:

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

And I think that you’ll agree that it does! Now we’ll take a look at some diagnostics:

nlm["ParameterTable"]
{pa,pb,pc}=nlm["BestFitParameters"]
ListPlot[nlm["FitResiduals"]]

See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2892436/
for a discussion of the mis-application of (R^2) in non-linear models. They argue for the use of no (R^2).

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)

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]
(* 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[
(* Use our current best guess for the parameters in the model to get estimates: *)
estimates=f[hours,alpha,beta,gamma];
(* Set up the design matrix *)
fpa=fa[hours,alpha,beta,gamma];
fpb=fb[hours,alpha,beta,gamma];
fpc=fc[hours,alpha,beta,gamma];
(* Linear regression of the design matrix and difference between data and the current model, to get new parameters *)
diffLinMod=LinearModelFit[
Transpose[{fpa,fpb,fpc,(cumulative-estimates)}],
{x,y,z},{x,y,z},
IncludeConstantBasis->False
];
(* And the new parameters are *)
{newa,newb,newc}=diffLinMod["BestFitParameters"];
(* Update the parameters using the new parameters *)
alpha=alpha+newa;
beta=beta+newb;
gamma=gamma+newc;
(* Print out the current state of the linear regression *)
Print[diffLinMod["RSquared"]];
Print[diffLinMod["ParameterTable"]],
(* and iterate *)
{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 model values ((OverHat[y])), as well as the residuals, and perhaps compute something like (R^2) values (I’ve not discovered how Mathematica computes  (R^2) in the non-linear case -- remember, some folks don’t think that it’s a good idea).

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! But it does seem odd....