Table of Advertising Expenditures, US, in millions of dollars (Crain, et al.). Cited in Mooney and Swift.

Andy Long
Spring, 2018

Years=Table[i,{i,1970,1989}];
Years=Years-1900;
AdExpenditures={19550,20700,23210,24980,26620,27900,33300,37440,43330,48780,53550,60430,66580,75850,87820,94750,102140,109650,118050,123930};
data=Transpose[{Years,AdExpenditures}];

lm=LinearModelFit[data,x,x,ConfidenceLevel->.95]
p1=ListPlot[data];
p2=Plot[lm[x],{x,Min[Years],Max[Years]}];
Show[p1,p2,Frame->True]
lm["RSquared"]

Let’s look at an exponential model, obtained by linearization of the problem. We’ll also compare it to the model estimated by non-linear regression.

ty=Log[AdExpenditures];
transdata=Transpose[{Years,ty}];

tm=LinearModelFit[transdata,x,x,ConfidenceLevel->.95]

p1=ListPlot[transdata];
p2=Plot[tm[x],{x,Min[Years],Max[Years]}];
Show[p1,p2,Frame->True]
tm["RSquared"]

{alpha,beta}=tm["BestFitParameters"]
CIs=tm["ParameterConfidenceIntervals"]
{a,b}={Exp[alpha],beta}
Exp[First[CIs]] (* so a is not in the center of the exponentiated CI *)
a

y=AdExpenditures;
ybar=Mean[y];
expModel[x_]:=a Exp[b x]
yhat=expModel[Years];

p1=ListPlot[data];
p2=Plot[expModel[x],{x,Min[Years],Max[Years]}];
Show[p1,p2,Frame->True]
RSquaredTrans=1-(y-yhat).(y-yhat)/((y-ybar).(y-ybar))

pres=ListPlot[Transpose[{Years,y-yhat}]]

Now I want to try NonlinearModelFit with the same model
to see how much difference there is (and to compare CIs against the transformed ones above).

Clear[aa,bb]
nlm=NonlinearModelFit[data,aa Exp[bb x],{{aa,a},{bb,b}},x]
{a,b}=nlm["BestFitParameters"]
CIs=nlm["ParameterConfidenceIntervals"]
 (* 
 Notice that these are centered on the parameter estimate. Also notice the difference in parameter values! The values for both parameters fall outside of the CI from the linearization!
 *)

p3=Plot[nlm[x],{x,Min[Years],Max[Years]}];
Show[p1,p3,Frame->True]
yhat=Table[nlm[Years[[i]]],{i,Length[Years]}];
nlm["RSquared"]
RSquaredTrans=1-(y-yhat).(y-yhat)/((y-ybar).(y-ybar))
(* so we're using different definitions of R-squared... and I can't find Mathematica's definition anywhere! :( *)

Now let’s try a power model:

ty=Log[AdExpenditures];
tx=Log[Years];
loglogdata=Transpose[{tx,ty}];

pm=LinearModelFit[loglogdata,x,x,ConfidenceLevel->.95]

p1=ListPlot[loglogdata];
p2=Plot[pm[x],{x,Min[tx],Max[tx]}];
Show[p1,p2,Frame->True]
pm["RSquared"]

{alpha,beta}=pm["BestFitParameters"]
{a,b}={Exp[alpha],beta}
y=AdExpenditures;
ybar=Mean[y]
powerModel[x_]:=a x^b
yhat=powerModel[Years];

p1=ListPlot[data];
p2=Plot[powerModel[x],{x,Min[Years],Max[Years]}];
Show[p1,p2,Frame->True]
RSquaredPower=1-(y-yhat).(y-yhat)/((y-ybar).(y-ybar))

pres=ListPlot[Transpose[{Years,y-yhat}]]

Now let’s try a logistic model (for which I’m going to use non-linear regression!):

Clear[a,b,c];nlm=NonlinearModelFit[data,c /(1+Exp[(a+b x)]),{{a ,Log[5]},{b,-2/80},{c,150000}},x]

Clear[a,b,c]
p1=ListPlot[data];
p2=Plot[nlm[x],{x,Min[Years],Max[Years]}];
Show[p1,p2,Frame->True]
ListPlot[nlm["FitResiduals"]]
params=nlm["BestFitParameters"]
params[[1]]
pa=a/.params[[1]]
pb=b/.params[[2]]
pc=c/.params[[3]]
pc
logisticModel[x_]:=pc /(1+Exp[(pa+pb x)])
yhat=logisticModel[Years];
nlm["RSquared"]
RSquaredLogistic=1-(y-yhat).(y-yhat)/((y-ybar).(y-ybar))