Keeling Analysis: non-linear regression

Andy Long
Spring, 2020

Throughout the non-linear version of R^2 is useless. I propose a “better” version, which suggests that the non-linear model with variable amplitude has the highest adjusted R^2 of 0.998882 (which just barely beats the R^2 of the linear model: R^2 of 0.998863). 

That model suggests that the amplitude is, indeed, varying (based on the parameters of the linear function for amplitude being distinctly non-zero, per the confidence intervals or p-values.

The period is, interestingly enough, not 1, based on the 95% confidence intervals (they do not contain 1 as a possibility).

Here’s one way you might Import the CO2 data: straight off the web.

dat=Import["ftp://aftp.cmdl.noaa.gov/products/trends/co2/co2_mm_mlo.txt","Table"];
DeleteCases[dat,{__String},{-1}];
dat[[96,1]]
dat[[71]]
(* First data line in the file: *)
Table[dat[[73,i]],{i,1,7}] 
dim=Dimensions[dat][[1]]
DecimalDate=Table[dat[[i,3]]-1958,{i,74,dim}] ;
CO2data=Table[dat[[i,5]],{i,74,dim}] ;
KeelingData=Transpose[{DecimalDate,CO2data}];
p1=ListPlot[KeelingData]

We’ll do some linear regression to get starting parameters:

Clear[a,b,c]
lm=LinearModelFit[
KeelingData,{x,x^2,Sin[2Pi x],Cos[2Pi x]},x]

{astart,bstart,cstart,d,e}=lm["BestFitParameters"]
ampstart=Sqrt[d^2+e^2]
x0start=-ArcTan[e/d]/(2 Pi)
lm["ParameterTable"]
lm["ANOVATable"]
lm["ParameterConfidenceIntervalTable"]
Show[p1,Plot[lm[x],{x,0,Max[DecimalDate]}]]

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit, and that produced by LinearModelFit:

SSTuncorrected=(CO2data).(CO2data)
fits=Table[lm[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
SSTotal=(CO2data-Mean[CO2data]).(CO2data-Mean[CO2data])
SSReg=(fits-Mean[CO2data]).(fits-Mean[CO2data])
SSReg/SSTotal
lm["RSquared"]

Now we’ll do non-linear regression to get the parameters:

Clear[a,b,c,amp,x0]
{astart,bstart,cstart,ampstart,x0start}
nlm=NonlinearModelFit[
KeelingData,
a+b x+c x^2 +amp Sin[2 Pi (x-x0)/Tval],
{{a ,astart},{b,bstart},{c,cstart},
{amp,ampstart},{x0,x0start},{Tval,1}},
x]
{aend,bend,cend,ampend,x0end,Tvalend}==nlm["BestFitParameters"]
nlm["ParameterTable"]
nlm["ANOVATable"]
nlm["ParameterConfidenceIntervalTable"]

What we notice is that the Tval CI does not contain 1. Hence 1 is not an acceptable value (at 95% confidence).

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit:

fits=Table[nlm[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
nlm["RSquared"]
nlm["AdjustedRSquared"]

(* This seems like the fairer comparison to the linear regression R^2. 
However it doesn't take into account the additional parameters 
(six for the nlm fit, and only five for the linear model): 
the more parameters you use, the better fit they should provide. *)
SSReg=(fits-Mean[CO2data]).(fits-Mean[CO2data])
ourSQ=SSReg/SSTotal

(* https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2 
Rsqadj=1-(1-Rsq)* (n-1)/(n-p-1)
*)
n=Length[KeelingData]
p=6
Rsqadj=1-(1-ourSQ)* (n-1)/(n-p-1)

ListPlot[Transpose[{DecimalDate,nlm["FitResiduals"]}]]
Sqrt[Mean[nlm["FitResiduals"]^2]]

Show[p1,Plot[nlm[x],{x,0,Max[DecimalDate]}]]

Not much change. Now let’s allow the amplitude of the sinusoidal variation be a linear function, and see if it’s significantly changing over time:

Clear[a,b,c,amp,x0,m]
{astart,bstart,cstart,ampstart,x0start}
nlmFinal=NonlinearModelFit[
KeelingData,
a+b x+c x^2 +(amp+m x) Sin[2 Pi (x-x0)/Tval],
{{a ,astart},{b,bstart},{c,cstart},
{amp,ampstart},{x0,x0start},{m,0},{Tval,1}},
x]
{aend,bend,cend,ampend,x0end, mend,Tvalend}=nlmFinal["BestFitParameters"]
nlmFinal["ParameterTable"]
nlmFinal["ANOVATable"]
nlmFinal["ParameterConfidenceIntervalTable"]

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit:

fits=Table[nlmFinal[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
nlmFinal["RSquared"]
nlmFinal["AdjustedRSquared"]

(* This seems like the fairer comparison to the linear regression R^2. 
However it doesn't take into account the additional parameters (seven, versus six for the nlm fit, and only five for the linear model): 
the more parameters you use, the better fit they should provide. *)
SSReg=(fits-Mean[CO2data]).(fits-Mean[CO2data])
ourSQ=SSReg/SSTotal
(* https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2 
Rsqadj=1-(1-Rsq)* (n-1)/(n-p-1)
*)
n=Length[KeelingData]
p=7
Rsqadj=1-(1-ourSQ)* (n-1)/(n-p-1)

Show[p1,Plot[nlmFinal[x],{x,0,Max[DecimalDate]}]]

phase=x0/.x0end
phaseInDaysFromNewYear=phase*365.25
periodFinal=1+period/.periodend
periodFinal= periodFinal*365.25

ListPlot[Transpose[{DecimalDate,nlmFinal["FitResiduals"]}]]
Sqrt[Mean[nlmFinal["FitResiduals"]^2]]

Let’s compare the models: A manipulate command let’s us see how the models are doing in their fit to the data:

Manipulate[
Module[{},
Show[
Plot[
{nlmFinal[x],nlm[x],lm[x]},
{x,start,start+yearspan},
PlotStyle->{Blue,{Orange,Dashed},Gray}]
,p1]
],
{start,0,Max[DecimalDate]},
{yearspan,1,10}
]

I didn’t ask you for the exponential model, but Flerlage, et al. did a nice job in their report, and concluded that this was the best model:

“Mathematica found that the best power model fit for the data is fit by the nonlinear model: CO2=257.184 +56.2922 e^ (0.0163503 * (years since 1958)) + 2.83537 * sin(2π ((years since 1958)-1.05748)] )”

Clear[a,b,c,amp,x0]
{astart,bstart,cstart,ampstart,x0start}={257.184,56.2922,0.0163503,2.83537,.05748}
nlmExpo=NonlinearModelFit[
KeelingData,
a+b Exp[c x] +amp Sin[2 Pi (x-x0)],
{{a ,astart},{b,bstart},{c,cstart},
{amp,ampstart},{x0,x0start}},
x]
{aend,bend,cend,ampend,x0end}=nlmExpo["BestFitParameters"]
nlmExpo["ParameterTable"]
nlmExpo["ANOVATable"]
nlmExpo["ParameterConfidenceIntervalTable"]

Show[p1,Plot[nlmExpo[x],{x,0,Max[DecimalDate]}]]
phase=x0/.x0end
phaseInDaysFromNewYear=phase*365.25

Examining the RSquared as calculated by the non-linear regression function NonlinearModelFit:

fits=Table[nlmExpo[x],{x,DecimalDate}];
SSModel=(fits).(fits)
SSModel/SSTuncorrected
nlmExpo["RSquared"]
(* This seems like the fairer comparison to the linear regression R^2. 
However it doesn't take into account the additional parameters: 
the more parameters you use, the better fit they should provide. *)
SSReg=(fits-Mean[CO2data]).(fits-Mean[CO2data])
ourSQ=SSReg/SSTotal
(* https://en.wikipedia.org/wiki/Coefficient_of_determination#Adjusted_R2 
Rsqadj=1-(1-Rsq)* (n-1)/(n-p-1)
*)
n=Length[KeelingData]
p=5
Rsqadj=1-(1-ourSQ)* (n-1)/(n-p-1)

Plot[nlmFinal[x]-nlmExpo[x],{x,0,Max[DecimalDate]}]
ListPlot[Transpose[{DecimalDate,nlmExpo["FitResiduals"]}]]

root=FindRoot[nlmExpo[x]==450,{x,60}]
1958+x/.root
root=FindRoot[nlmFinal[x]==450,{x,60}]
1958+x/.root
root=FindRoot[nlm[x]==450,{x,60}]
1958+x/.root
root=FindRoot[lm[x]==450,{x,60}] 
1958+x/.root