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

Alyssa Farmer (Most of the code attributed to Andy Long)
Spring, 2018

The linear model does estimate the ad expenditures with some accuracy, but it is not the best model. Looking at the plot of the model and the residual plot, it is obvious that the residuals are uneven and another term or a different model is needed in order to account for curvature. The value for R^2 is 0.95307, meaning that the model does explain much of the variation in the data, but we can do better.

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"]

lm["ParameterConfidenceIntervals"]

ListPlot[lm["FitResiduals"], Filling->Axis]

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.

After log transposing the ad expenditures, the linear fit to the line looks much more reasonable, although not perfect (the residuals look like they have a pattern). When we untransform to get the exponential model, the model fits the data much better. The value of R^2 is now 0.987, which is better, but we already were expecting it to be better than the linear model, so it does not mean much. However, the confidence intervals on the parameters of the transformed model do not contain zero, so a transformation does fit a positive linear trend. The residual plot of the untransformed data looks like it has non-constant variance, which is a problem.

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"]

ListPlot[tm["FitResiduals"], Filling->Axis]

{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:

The power model seems to be similar to the exponential model. We first transform both the years and the ad expenditures, and the parameter confidence intervals on the transformed linear model again did not contain zero, so the transformation was necessary and helpful to fit the linear model better. The untransformed model and the residuals look very similar to that of the exponential model, meaning that neither of the two models necessary fit the data better. However, more issues are discussed below that address the time scale on the two models.

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];
pm["ParameterConfidenceIntervals"]

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}]]

Comparing the exponential and power models

There is an issue with the power model that does not exist with the  exponential model. With the power model, when we change the starting value of t to where it represents 70 as 1970 instead, there is a big difference in the equation for the final model. The exponent on the exponential model stays the same whether we have a shift in the time scale or not (b=0.105), but with the power model, there is a difference. If the year is represented as 70, then b=8.354, but if the year is represented as 1970, then b=208.817. This is an issue  because interpretation of the equations becomes much more difficult, and predictions will also be different depending on how the time is represented. Either way, the starting value needs to be 0 at t=0, but we need the models to fit both data, and this results in a bigger exponent for one of them. The meaning of each exponent does not become very clear in either situation, which is why this model may not be the best.

As far as interpreting the parameters goes, the value of a for the exponential model is equivalent to the estimated ad expenditures in the year 1900, where the value of a for the power model is the estimated ad expenditures for the year 1901. The value of b is more difficult to interpret in each situation, but in the power model, it has to do with how quickly the estimated values increase, and a larger value for b would increase this rate. The value of b in the exponential model just determines what the “base” value is. So for every additional year, the estimated ad expenditures increases by a multiple of e^b.

Years=Years+1900;

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

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

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

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

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

The logistic model seems to be the true best model, in my opinion. It looks like the data may follow an S shape, and the logistic model accounts for this. The residuals at least look like they have constant variance, although they also look like they have a pattern to them. So again, the model is not perfect, but it is at least useful and may be best for prediction, accounting for a “leveling-off” effect in the future.

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

Years=Table[i,{i,70,89}];
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))