Looking for Patterns and Applying the Method of Least Squares to Real Data

Introduction

OBJECTIVE: Minimize the sum of squared residuals to fit an arbitrary function to a set of data.

We look for patterns in real data sets including medical, socio-economic, meteorological, and pollution statistics. Does the least squares method apply only to finding lines of best fit, or does it apply to any type of curve that you wish to fit to a set of data points? We begin by analyzing data with standard fit functions and then utilize the calculus definition of best fit to find a nonstandard fit function.

Technology Guidelines

NOTE:  If you have just finished a module, restart Mathematica or Quit Kernel under the Evaluation pull-down menu before executing a new module.
TO OPEN CELLS, click on the arrow to the left of the cell or put your cursor on the right cell bracket and double click. To open all cells, highlight the bracket on the far right and select Grouping under the Cell pull-down menu and select Open All Subgroups.
TO STOP AN EXECUTION
Select the Evaluation pull-down menu and click on Abort Evaluation.
ORDER OF EXECUTION
Execute cells in the order given. Do not skip any Input cells within a given notebook.
SAVING NOTEBOOKS
You can save anytime to any directory you choose, and it is wise to save often. However, before you do your final save to reduce the size of your file, it is a good idea to delete all your output by selecting the Delete All Output selection under the Cell pull-down menu.
EXPERIENCING MAJOR PROBLEMS
Save if appropriate, then shut down Mathematica and start it up again.

Part I: Time Series Data

Thirty-Two Years of Salaries - Men Versus Women

The following data represent the mean income for men and women in the United States from 1967 to 1998. The study includes people 15 years old and over. All income is in 1998 consumer price index adjusted dollars.
SOURCE:  March Current Population Survey
PREPARED BY:
Income Statistics Branch/HHES Division,  U.S. Bureau of the Census, U.S. Department of Commerce
Washington, D.C.  20233-8500,  (301) 457-3242

menearnings={"Men's Salaries",27232,28209,29338,29290,29389,31214,31587,30495,29784,30168,30635,31180,31039,29916,29419,29180,29182,30027,30805,31956,32210,32635,33324,31978,31074,30670,32143,32887,33126,33553,34794,36315};
womenearnings={"Women's Salaries",11500,11631,11997,12195,12412,12930,12918,12868,12889,13172,13437,13313,13070,13207,13253,13758,14148,14805,15174,15729,16301,16703,17119,17085,17027,17070,17506,17846,18183,18790,19511,20462};
years={Years,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998};
data=Table[{years[[i]],menearnings[[i]],womenearnings[[i]]},{i,1,Length[years]}]//TableForm

As was probably expected, men's mean income exceeds women's mean income, but we can get a better perspective if we visualize these data through plots.

men=Table[data[[1,i,2]],{i,2,Length[years]}];
women=Table[data[[1,i,3]],{i,2,Length[years]}];
pmen:=ListPlot[men,Joined->True,PlotStyle->RGBColor[0,1,0],DisplayFunction->Identity,AxesLabel->{"Years Since 1966","Salaries"},AxesOrigin->{0,10000}];
pwomen:=ListPlot[women,Joined->True,PlotStyle->RGBColor[0,0,1],DisplayFunction->Identity,AxesLabel->{"Years Since 1966","Salaries"},AxesOrigin->{0,10000}];
Show[pmen,pwomen,DisplayFunction->$DisplayFunction]
Print["Salaries for men are in green and salaries for women are in blue."]

Neither data set follows a straight line, but both have an upward trend. One way we can examine those trends is to compute the best-fit lines for each data set. We will do this using the Fit function in Mathematica.  Note that you choose a linear fit by specifying the {x,1}.

Print["The average salary of males is approximately ",fmen=Fit[men,{x,1},x]," and"]
Print["the average salary of females is approximately ",fwomen=Fit[women,{x,1},x],","]
Print["where x represents the years since 1966"]
fitplots=Plot[{fmen,fwomen},{x,0,33},PlotStyle->{RGBColor[0,1,0],RGBColor[0,0,1]},DisplayFunction->Identity];
Show[pmen,pwomen,fitplots,DisplayFunction->$DisplayFunction]

What do these lines tell you about the rate at which salaries are increasing for men versus women? You could verify that these lines would intersect in about 183 years. Why could this be considered a worthless estimate?

You Try It: Part I

Lake Pollution

The following data sets represent weekly readings of aluminum pollutant levels in Lake Erie over a two-year period. What observations can you make about the amounts of aluminum pollutant? Does the amount of pollution follow a pattern? The first data set represents the data in the order that they were recorded. We then order those data to better analyze the distribution of the amounts of pollutant.

pollutant={59.7390,43.6978,45.4508,55.6974,43.6578,56.6556,49.9503,59.2501,59.2778,30.8085,50.9860,59.0104,57.1471,48.4696,47.4461,37.0215,52.0141,53.8839,35.8155,43.8406,58.4949,56.7300,45.6032,37.7449,45.7960,44.0860,55.0693,30.5343,57.3854,30.2706,42.1799,44.9115,40.9325,57.4171,30.2881,46.1233,36.4139,46.7199,36.3994,42.1318,54.0219,34.4747,49.2455,43.0440,58.2225,57.4492,53.2136,30.9330,48.8363,32.6732,53.6791,31.7362,35.7099,39.6966,48.3676,42.1876,53.2427,47.2459,31.4635,43.3965,44.8305,46.8731,39.5781,37.0135,37.8624,37.1184,36.5494,41.2847,48.9843,50.6743,48.5370,33.1112,54.3590,56.7471,34.9079,43.4007,34.6556,42.0995,55.4403,42.9701,56.2188,39.4222,33.5046,38.6973,59.6348,38.3582,46.5821,36.8290,35.3832,56.9938,33.0390,37.8167,51.4268,37.3356,30.3771,37.7568,43.9597,50.1189,47.6849,53.0739,34.6360,36.0122,31.8986,49.6269};
orderedpollutant=Sort[pollutant]

Plot the ordered data, and see what the plot tells you about the pattern of the data.

pollutantplot=ListPlot[orderedpollutant,Joined->True]

Look for fit functions for these data. Replace the terms in red with an appropriate expression.

pollutantfit=Fit[orderedpollutant,{(x^2),x,1},x]
fitplot=Plot[pollutantfit,{x,0,102},DisplayFunction->Identity,PlotStyle->RGBColor[1,0,0]];
Show[pollutantplot,fitplot,DisplayFunction->$DisplayFunction]

The good fit for the linear model could imply that the pollutant distribution is somewhat uniform. What would uniform mean in this context?

Part II: Looking for Cause and Effect: Temperature as a Function of Latitude

Following are data showing 56 cities in the U. S., their January temperatures ((^0)F) averaged over a thirty-year period, and their latitudes.

Off[Set::write]
city={City,MobileAL,MontgomeryAL,PhoenixAZ,LittleRockAR,Los AngelesCA,San FranciscoCA,DenverCO,NewHavenCT,WilmingtonDE,WashingtonDC,JacksonvilleFL,KeyWestFL,MiamiFL,AtlantaGA,BoiseID,ChicagoIL,IndianapolisIN,DesMoinesIA,WichitaKS,LouisvilleKY,NewOrleansLA,PortlandME,BaltimoreMD,BostonMA,DetroitMI,MinneapolisMN,St.LouisMO,HelenaMT,OmahaNE,ConcordNH,AtlanticCityNJ,AlbuquerqueNM,AlbanyNY,NewYorkNY,CharlotteNC,RaleighNC,BismarckND,CincinnatiOH,ClevelandOH,OklahomaCityOK,PortlandOR,HarrisburgPA,PhiladelphiaPA,CharlestonSC,NashvilleTN,AmarilloTX,GalvestonTX,HoustonTX,SaltLakeCityUT,BurlingtonVT,NorfolkVA,SeattleWA,SpokaneWA,MadisonWI,MilwaukeeWI,CheyenneWY};

januarytemp={January Temperature,44,38,35,31,47,42,15,22,26,30,45,65,58,37,22,19,21,11,22,27,45,12,25,23,21,2,24,8,13,11,27,24,14,27,34,31,0,26,21,28,33,24,24,38,31,24,49,44,18,7,32,33,19,9,13,14};

latitude={Latitude,31.2,32.9,33.6,35.4,34.3,38.4,40.7,41.7,40.5,39.7,31,25,26.3,33.9,43.7,42.3,39.8,41.8,38.1,39,30.8,44.2,39.7,42.7,43.1,45.9,39.3,47.1,41.9,43.5,39.8,35.1,42.6,40.8,35.9,36.4,47.1,39.2,42.3,35.9,45.6,40.9,40.9,33.3,36.7,35.6,29.4,30.1,41.1,45,37,48.1,48.1,43.4,43.3,41.2};

city data=Table[{city[[i]],januarytemp[[i]],latitude[[i]]},{i,1,Length[city]}]//TableForm

First, do a scatterplot to see if there is a relationship between the January temperatures and latitude.

scatter=ListPlot[Table[{latitude[[i]],januarytemp[[i]]},{i,2,Length[latitude]}],AxesLabel->{"Latitude",January Temperature},PlotStyle->PointSize[0.015`]]

Since the scatterplot appears to have a somewhat linear pattern, that shows the average January temperature dropping as the latitude increases, we will investigate this relationship further. We will use Mathematica's built-in Fit function to get results. Note that you choose a linear fit by specifying the {x,1}. We can write the line of best fit as follows.

yfit=Fit[Table[{latitude[[i]],januarytemp[[i]]},{i,2,Length[latitude]}],{x,1},x]

Now we can plot this line with our scatterplot.

linearfit=Plot[yfit,{x,20,50},DisplayFunction->Identity];
Show[scatter,linearfit]

Can you identify cities that deviate from the pattern? Might you speculate on the reason for the deviation? Name one important factor that might explain at least a part of this deviation.

You Try It: Part II

Here are two data sets you can explore. You can enter your own data sets and then re-execute the computations to explore relationships.

Age and EEG

The electroencephalogram (EEG) is a device used to measure brain waves. Neurologists have found that the peak EEG frequency in children increases with age. The data below represent the results of a study of children age 2 to 16, and the EEG readings (in hertz) represent the average peak EEG frequencies for each age group. It would be reasonable to think of age as the independent variable and eeg as the dependent variable.

age={2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
eeg={5.33`,5.75`,5.8`,5.6`,6.`,5.78`,5.9`,6.23`,7.28`,7.06`,7.6`,7.45`,8.23`,8.5`,9.38`};
eegdata=Table[{age[[i]],eeg[[i]]},{i,1,Length[age]}];
TableForm[eegdata]
scatter=ListPlot[eegdata,AxesOrigin->{0,5},AxesLabel->{"age","eeg"},PlotStyle->PointSize[0.01`],PlotStyle->PointSize[0.03`]]

Alter the terms in red until you get a good fit.

eegfit=Fit[eegdata,{Exp[0.2 x],1},x]
fitplot=Plot[eegfit,{x,0,17},DisplayFunction->Identity,PlotStyle->RGBColor[1,0,0]];
Show[scatter,fitplot,DisplayFunction->$DisplayFunction]

Reliability of Construction Material

The Canadian Geotechnical Journal (Aug., 1985) reported on a study that was conducted to investigate the reliability of the use of fragmented Queenston Shale, a compaction shale, as a rockfill construction material. In particular, the researchers wanted to estimate the stress-strain relationship of the fragmented material. Their study yielded the following results, where axial strain (x) is given as percentages and deviatoric stress (y) is given in kPa.  You might want to consider a quadratic fit for these data.

strain={1.`,2.8`,4.2`,6.`,7.5`,9.`,10.5`,13.5`,16.7`,19.8`,23.`,26.`,27.5`};
stress={500,2000,2750,3500,4375,4875,5250,6000,6625,7000,7125,7000,7125};
ss=Table[{strain[[i]],stress[[i]]},{i,1,Length[strain]}];
ssplot=ListPlot[ss,AxesOrigin->{0,500},AxesLabel->{"strain","stress"},PlotStyle->PointSize[0.015`]]

Alter the terms in red until you get a good fit.

ssfit=Fit[ss,{Log[x],1},x]
fitplot=Plot[ssfit,{x,1,30},DisplayFunction->Identity,PlotStyle->RGBColor[1,0,0]];
Show[ssplot,fitplot,DisplayFunction->$DisplayFunction]

Part III: Finding a Least Squares Curve Using Calculus

Data were collected concerning a specific genetic characteristic that researchers believed followed a particular exponential pattern. After observing a scatterplot of the data, they looked for a fit function of the form y=a (e^bx)+c. Because this is not a standard fit function, it is necessary to use the minimization techniques from multivariable calculus to determine a, b, and c from scratch. In all there were 2733 pieces of data collected, and those data were placed into eight categories. In the set defined below, the first entry represents the frequency, the second the categorical variable x, and the third the observed genetic measurement y.

Clear[x,y,a,b,c]
dataset={{579,1,38.08},{1021,2,29.70},{607,3,25.42},{324,4,23.15},{120,5,21.79},{46,6,20.91},{17,7,19.37},{9,8,19.36}};
Print["freq. count observations"]
Print[TableForm[dataset]]
y[x_]:=a Exp[b x]+c

We begin our process by writing the formula for the sum of the squares of the vertical distances between the curve of best fit and the observed data. This will be a function of a, b, and c, so we will find the first three partial derivatives and solve for the values of a, b, and c when we set those partial derivatives equal to 0.

Clear[a,b,c,f,fa,fb,fc]
g[a_,b_,c_]:=Sum[dataset[[i,1]]*(((dataset[[i,3]]- a Exp[dataset[[i,2]]*b]-c))^2),{i,1,Length[dataset]}]
ga[a_,b_,c_]=D[g[a,b,c],a];
gb[a_,b_,c_]=D[g[a,b,c],b];
gc[a_,b_,c_]=D[g[a,b,c],c];

Now we will set up our equations to solve. Due to previous studies, we had an idea of a range of values for the parameters a, b and c. For that reason, we can use the command FindRoot and give some initial starting values, rather than call upon the Solve or NSolve command.

\[MathematicaIcon] About Mathematica

Clear[a,b,c]
oursolution=FindRoot[
			{ga[a,b,c]==0,gb[a,b,c]==0,gc[a,b,c]==0},
			{a,40.},{b,-1.},{c,20.}]

Show y[x] for those values of a, b, and c.

yhat[x_]=y[x]/.oursolution

Look at the results graphically.

p1=Plot[yhat[x],{x,0,12},PlotRange->{15,55},DisplayFunction->Identity];
pts=Table[{dataset[[i,2]],dataset[[i,3]]},{i,1,Length[dataset]}];
p2=ListPlot[pts,PlotStyle->PointSize[0.01`],DisplayFunction->Identity];
Show[p1,p2,DisplayFunction->$DisplayFunction]

This looks like a good fit. Compute the sum of the squares of the errors.

yest=Table[yhat[dataset[[i,2]]],{i,1,Length[dataset]}];
yact=Table[dataset[[i,3]],{i,1,Length[dataset]}];
wts=Table[dataset[[i,1]],{i,1,Length[dataset]}];
Print["The sum of the squares of the vertical deviations for the 2723 data points is ",sse=wts.(((yact-yest))^2)]
Print["The sum of the squares of the vertical deviations written as percentages for the 2723 data points is ",sse=wts.((((yact-yest))^2)/yest)]

This last entry is important for running a statistical test on how good a fit this model is. Because of our calculus work, we know that this is the best fit model for a function of this type.

\[MathematicaIcon] About Mathematica

The disadvantage of using the Solve or NSolve commands here is that they will look for all solutions, both real and complex. This can take a long time (10 to 15 minutes) because there are many complex solutions to this system of equations. The FindRoot command will only look for real solutions if the equations are real and the starting points are real. FindRoot works somewhat like Newton's method. You will notice that it yields only approximate solutions, and it may not converge to a solution at all. 
 Go back.