Miniproject
For The 
Mathematica  Workshop

Graph of 4 regression models on a scatter plot of bivariate data

Carol Castellon

"Experimenting With the Regression Packages"

Brief Description

This project takes an arbitrary set of bivariate data (paired data) and looks at vaious ways to "fit" a model to the data.  Models which were fitted include linear, quadratic, cubic, trigonometric, logrithmic, and a sum of non-integral-power ploynomial-like functions.
					Part 1 of this project uses the "Fit" command and allows the user to compare the resulting models on the scatter plot visually.  The statistical packages are not used in this section of the project.
					Part 2 of this project uses the Linear Regression Package and the "Regress" command to obtain a statistical analysis of the models, and uses this analysis to compare three different models. 
					Part 3 of this project uses the NonLinear Regression Package to obtain several non-linear models and alter the method used for the fit.

Experimenting with the
Regression Packages

<<Graphics`Colors`

A.					Bivariate Data and Graphs

Given this set of paired data,
(1,4)				(2,3)				(3,3)				(4,5)				(5,8)					(5,5)					(6,7)

First:   look at the scatter plot of the data, 

Then:   try to "fit" a line or curve to the data.

list={{1,4},{2,3},{3,3},{4,5},{5,8},{5,5},{6,7}}
points = ListPlot[list,PlotRange->{{-1,7},{-1,9}},
Prolog->PointSize[.015],PlotStyle->VenetianRed];

Here is the scatter plot of the paired data.

The first "fit", a polynomial of degree 0, is really the mean of the dependent variable, y.  The mean of y is 5.  A horizontal line, y=5,
"fits" the data when linear regression is not appropriate, i.e. when only a point estimate is needed.

poly0 =Fit[list,{1},x]
pt = Plot[poly0,{x,-1,7},PlotStyle->Navy,
DisplayFunction->Identity];
Show[points,pt,Prolog->PointSize[.015]];

In most cases, the line of regression, or "line of best fit" 
(using Least squares criterian) is desired.  In this example,
the "line of best fit" is y = 1.94118 + 0.823529 x

poly=Fit[list,{1,x},x]
bestfit = Plot[poly,{x,-1,7},PlotStyle->Navy,
DisplayFunction->Identity];
Show[points,bestfit,Prolog->PointSize[.015]];

A second degree polynomial also "fits" this data.  
The model is y = 4.08197 - 0.783763 x + 0.231069 x¤

poly2 =Fit[list,{1,x,x^2},x]
bestfit2 = Plot[poly2,{x,-1,7},PlotStyle->Green,
DisplayFunction->Identity];
two = Show[points,bestfit2,Prolog->PointSize[.015]];

Look at both linear and quadratic models on the scatter plot.

combine = Show[points,bestfit,bestfit2,
Prolog->PointSize[.015]];

A third degree polynomial may also "fit" the data.
The model is:  y = 8.37688 - 6.09859 x + 1.96411 x¤ - 0.164154 x‹

poly3 =Fit[list,{1,x,x^2,x^3},x]
bestfit3 = Plot[poly3,{x,-1,7},PlotStyle->HotPink,
DisplayFunction->Identity];
Show[points,bestfit3,Prolog->PointSize[.015]];

Look at al three models (linear, quadratic, cubic) on the scatter plot.

Three = Show[points,bestfit,bestfit2,bestfit3];

One thing is clear, any of the three models will make a "good" prediction within the observed values (range) of x.

Now try a trig "model".

Clear[trig,x]
trig[x_] =Fit[list,{1,Cos[(2 Pi/14) x],
Cos[4 Pi/14 x],Cos[6 Pi/14 x],Cos[8 Pi/14 x]},x]
besttrg = Plot[trig[x],{x,-1,7},PlotStyle->Purple,
DisplayFunction->Identity];
Show[points,besttrg,Prolog->PointSize[.015]];

The trig model is:

Visually, it appears the trig model is the best "fit" for this data.


Look at all models on the scatter plot.

Four = Show[points,bestfit,bestfit2,bestfit3,besttrg];

Once again, any of the models will make a "good" prediction within the observed values (range) of x.  Note the cubic and the trig models are nearly the same.  Can you explain why?

B.					Regression Package

Initialization:

<<Statistics`LinearRegression`
Off[General::spell1];

The Statistics Package may help us decide which model "fits" best. 

Recall that the linear model was y = 1.94118 + 0.823529 x

data := list
regone = Regress[data,{1,x},x,IncludeConstant->False,
		OutputList->{BestFit,PredictedResponse,FitResiduals}]

predict1 := PredictedResponse /. regone
errors1 := FitResiduals /. regone

This table compares the linear model to a smaller model consisting of only a constant.  The F-ratio and the P-value indicate the linear model is a significantly better "fit",  i.e., a linear model is better for predicting the dependent variable within the observed range of values of the independent variable.

Recall that the second degree polynomial model was
 y = 4.08197 - 0.783763 x + 0.231069 x¤
Now we will look at the table for the quadratic model.

regtwo = Regress[data,{1,x,x^2},x,IncludeConstant->False,
		OutputList->{BestFit,PredictedResponse}]

predict2 := PredictedResponse /. regone
errors2 := FitResiduals /. regone

Once again, the table provides evidence that a quadratic model is a better "fit" than a constant.  Later, we will try to find out how the linear compares to the quadratic.

Recall that the third degree polynomial was:
y = 8.37688 - 6.09859 x + 1.96411 x¤ - 0.164154 x‹

regtwo = Regress[data,{1,x,x^2,x^3},x,IncludeConstant->False,
		OutputList->{BestFit,PredictedResponse}]

The ANOVA Table provides a comparison of the given model to a smaller model which I reduced to the data.  The F-test compares the two models by the ratio of their mean squares.  If the value of F is large, the null hypothesis supporting the smaller model is rejected.  In all three cases above, the smaller model (model of data only) was rejected (which was expected).  The coefficient of determination, R¤, is a summary statistic that describes the relationship between the regressors and the response variable.  AdjustedSquared gives an adjusted value that you can use to compare subsequent subsets of models.  The AdjustedSquared for the above models were:
linear:  0.937295
Quadratic:  0.939489
Cubic:  0.944307
which would seem to imply that all three models are nearly the same in their ability to predict the dependent variable within the observed range of the independent variable.

We might get some insight by looking at the plot of predicted values against residuals for all three models on the same axes.

opt1 = ListPlot[Transpose[{predict1,errors1}],
		PlotStyle->Navy,
		DisplayFunction->Identity];
opt2 = ListPlot[Transpose[{predict2,errors2}],
		PlotStyle->{PointSize[0.03],Green},
		DisplayFunction->Identity];				
opt3 = ListPlot[Transpose[{predict3,errors3}],
		PlotStyle->{PointSize[0.05],Thistle},
		DisplayFunction->Identity];
Show[opt3,opt2,opt1,DisplayFunction->$DisplayFunction];

What this shows is that all three models are the same in their ability to predict the dependent variable within the observed range of the independent variable!!!

c.					Nonlinear Fit Regression Package

<<Statistics`NonlinearFit`

The idea here is to experiment with various non-linear models of any variety, to find out the parameters which would lead to the best fit of the given model.  The model must be input, ie., the computer doesn't create one!
Using the same data set, (1,4) (2,3) (3,3) (4,5) (5,8) (5,5) (6,7)
Note the exponents the package generated to "fit" the model.
You are to compare the results by observing the graphs generated, keeping in mind that any of the models will "fit" the data.

The first model tried is of the form:  y = axı

params:= NonlinearFit[data,a x^b,x,{a,b}]
Clear[nonlin,x]
nonlin[x_] = params[[1,2]] x^params[[2,2]]
try1 = Plot[nonlin[x],{x,0,7},PlotStyle->Green,
DisplayFunction->Identity];
Show[points,try1,Prolog->PointSize[.015]];

The next model tried is of the form:  y = axı + cx + d

params5:= NonlinearFit[data,a x^b + c x + d,x,{a,b,c,d}]
Clear[nonlin5,x]
nonlin5[x_] = params5[[1,2]] x^params5[[2,2]] + params5[[3,2]] x  + params5[[4,2]]
try5 = Plot[nonlin5[x],{x,0,7},PlotStyle->Green,
DisplayFunction->Identity];
Show[points,try5,Prolog->PointSize[.015]];

The next model tried is of the form y = b cos(ax) + c

params2 := NonlinearFit[data,b Cos[a x] + c,x,{a,b,c}]
Clear[nonlin2,x]
nonlin2[x_] = params2[[2,2]] Cos [params2[[1,2]] x] + params2[[3,2]]
try2 = Plot[nonlin2[x],{x,0,7},PlotStyle->Violet,
DisplayFunction->Identity];
Show[points,try2,Prolog->PointSize[.015]];

The next model tried is of the form y = a log(bx) + c

params3:= NonlinearFit[data,a Log[b x] + c,x,{a,b,c}]
Clear[nonlin3,x]
nonlin3[x_] = params3[[1,2]] Log[params3[[2,2]] x] + params3[[3,2]]
try3 = Plot[nonlin3[x],{x,.25,7},PlotStyle->Navy,
DisplayFunction->Identity];
Show[points,try3,Prolog->PointSize[.015]];

This package is reliable when a linear model y = a + bx is tried.
Recall the "line of best fit" was y = 1.94118 + 0.823529 x

params4:= NonlinearFit[data,a+b x,x,{a,b}]
Clear[nonlin4,x]
nonlin4[x_] = params4[[1,2]] + params4[[2,2]] x
try4 = Plot[nonlin4[x],{x,0,7},PlotStyle->Turquoise,
DisplayFunction->Identity];
Show[points,try4,Prolog->PointSize[.015]];

The following is an experiment with the options available in the non-linear fit regression package.  The first options tried changes the method used in the fit, using the model y = axı + cx + d again.

params6:= NonlinearFit[data,a x^b + c x + d,x,{a,b,c,d},Method->FindMinimum]
Clear[nonlin6,x]
nonlin6[x_] = params6[[1,2]] x^params6[[2,2]] + params6[[3,2]] x  + params6[[4,2]]
try6 = Plot[nonlin6[x],{x,.00001,7},PlotStyle->Magenta,
DisplayFunction->Identity];
Show[points,try6,Prolog->PointSize[.015]];

This gives very different parameters from the earlier result,
y = 2.52558 + 0.127878x + 0.195419x∕˘ﬁ‡ﬁ·∕
(Curiously, this is the only model whose parameters changed significantly when a different method was used).  Both graphs are plotted together to show the differences.

Show[points,try5,try6,Prolog->PointSize[.015]];

Conclusions:
					Ultimately, to do a good regression analysis, we really need more than seven observations.  With a sample size of seven, regression analysis using any model would work.  So it is not surprizing that any of the models "fit" the data in terms of prediction value!
					The value of this project is in producing the codes which could be used for other data sets.