A Mathematica Introduction to Simple Linear Regression

The author describes an instructional laboratory approach to the study of least squares estimates.

by Evan Fisher

The combination of graphical and analytical capabilities of Mathematica makes it of great pedagogical potential in undergraduate statistics classes. It is of particular value in the standard calculus-based mathematical statistics course. For most students, this is their first exposure to probability and statistics. At this level, the mathematical development of the theory of statistics has a  deterministic quality that can obscure the sample variation that the theory is intended to model. Mathematica can provide an immediate visual connection between this sample variation and the underlying theory.

	 At various points throughout the semester, I incorporate computer laboratory sessions combining computer simulations and graphics to motivate and illustrate theoretical concepts in the mathematical statistics course. They are designed to supplement the formal development of the subject in the class and in the textbook (in the two most recent years, the textbooks by Rice[2] and Larsen and Marx[1] were used). A laboratory class session is a regular class meeting that takes place in the department of mathematics' computer laboratory. In the example described here, students work through a Mathematica notebook at their own machines. However, the pace at which this takes place is controlled through a combination of  discussion among students and between students and the instructor. These forced breaks from the computer activities highlight and explain their connection to the theory.

	This article describes the content of a laboratory session focusing on least squares regression under the normal simple linear regression model.	We consider the model of the form

(Subscript[y, i])=A + B (Subscript[x, i])+ (Subscript[ε, i])

where (Subscript[ε, i]) ∼ N(0, (σ^2)) are independent for i = 1, 2, 3, ... n and the (Subscript[x, i])'s are assumed fixed. Prior to the laboratory session, students have been introduced to the least squares method of fitting a line to a set of points and have derived, in the usual calculus-based manner, the least squares estimates for A and B. These are denoted here by a and b, respectively, and in one of the standard formulations are expressed as

\(TraditionalForm\`a = y\&_ - b\ x\&_\)

and

\(TraditionalForm
    \`b = \(∑\+\(i = 1\)\%n\((x\_i - x\&_)\)\ 
            \((y\_i - y\&_)\)\)\/\(∑\+\(i = 1\)\%n\((x\_i - x\&_)\)\^2
            \)\)

The laboratory session has two primary purposes. One is to emphasize the nondeterministic aspect of the linear model through random observations that are generated by the model and in the resulting least squares line which is a function of the observed data. The computer simulations are also designed to demonstrate that a and b are statistics as well as to  suggest characteristics of them as estimators of A and B.

The Laboratory Session

The least squares line and the regression line

We utilize one of the built-in Mathematica packages that permits the generation of random observations from the normal distribution.

<<Statistics`NormalDistribution`

The following function allows the user to generate a sample of observations from the model with any set of parameters. The list of fixed values of the independent variable x at which each response y will be observed is named "indep".

\(sample[A_, B_, sigma_, indep_List] := 
	
      Table[{indep[[i]], 
			
          Random[NormalDistribution[
              A + B\ indep[[i]], 
              sigma]]}, 
		{i, Length[indep]}]\)

In order to demonstrate a typical laboratory session, assume data is to be observed from the linear model

(Subscript[y, i]) = 3 (Subscript[x, 1]) + (Subscript[ε, 1]) + 1

where (Subscript[ε, i]) ∼ N(0, (σ^2)) with σ = 5. Suppose that ten independent observations will be taken at the values of x identified by the list {1, 2, 3, ..., 10}. We set the parameters to the appropriate values and generate a set of data {((Subscript[x, i]),(Subscript[y, i])): i = 1, 2, 3, ... 10}, that we name "sample1":

A=1;B=3;sigma=5;indep=Range[10];

sample1=sample[A,B,sigma,indep]

The data is visualized as a scatter plot:

\(\(scatter = 
      ListPlot[sample1, 
		
        PlotRange -> {{Min[indep] - 1, Max[indep]}, Automatic}]; \)\)

The students are reminded that, in practice, the true parameters of the model are unknown and estimates of the parameters A and B will be made from the observed data. To find the best fitting line using the least squares criterion, the following command is executed:

line1[x_]=Fit[sample1,{1,x},x]

That is, the equation of the least squares line based on this set of data is 
y = 0.966688 + 3.20227 x.

A plot of the data, the least squares line, and the regression line together illustrates the connection between them. The least squares line is shown in blue and the regression line is shown as a dashed red line.

\(plot2 = 
      Plot[{line1[x], A + B\ x}, 
		{x, Min[indep] - 1, Max[indep] + 1}, 

        		PlotStyle 
          -> {{RGBColor[0, 0, 1]}, 

            				{Thickness[0.005], Dashing[{0.02}], RGBColor[1, 0, 0]}}, 
        
		DisplayFunction -> Identity]\)

Show[scatter,plot2]

This results in a clear visual representation that connects three important conceptual parts of the theory: the linear, deterministic part of the model; the observed data, which reflects the random component of the model; and the least squares line determined from the data that serves as an estimate of the linear component of the model.

We now begin to emphasize the sample dependency of the least squares estimates. Assuming the same linear model and the same fixed values of the independent variable, students obtain a second set of observed data ("sample2"). The least squares line based on this data is determined, and a plot of this line (in green), the previous least squares line( in blue) generated by the first set of data, and the regression line (dashed and in red) representing the fixed component of the model are displayed together.

sample2=sample[A,B,sigma,indep]

line2[x_]=Fit[sample2,{1,x},x]

\(plot3 = 
      Plot[line2[x], 
		{x, Min[indep] - 1, Max[indep] + 1}, 
		
        PlotStyle -> {RGBColor[0, 1, 0]}, 
		
        DisplayFunction -> Identity]\)

\(Show[plot2, plot3, 
	DisplayFunction -> $DisplayFunction]\)

The two samples, both generated with the same fixed values of the independent variable are different—a result of the random component of the model. The least squares estimates of A and B are sample dependent. That is, a and b are statistics whose values vary from sample to sample. The change in the least squares line from one set of data to the next is evidence of this.

	The following short program, which combines elements of the previous steps,  generates a set of data from the model and displays the resulting scatter plot, least squares line (in blue) and the regression line (dashed in red) in one plot.

\(least[A_, B_, sigma_, indep_] := 
	
      \((sample1 = sample[A, B, sigma, indep]; 
		
        scatter = 
          ListPlot[sample1, 
				
            PlotRange 
              -> {{Min[indep] - 1, 
							Max[indep]}, 
                Automatic}, 
				PlotStyle -> PointSize[0.02], 

            				DisplayFunction -> Identity]; 
		
        line1[x_] := Fit[sample1, {1, x}, x]; 
		
        plot1 = Plot[{line1[x], A + B\ x}, 

            				{x, Min[indep] - 1, Max[indep] + 1}, 
				
            PlotStyle 
              -> {{RGBColor[0, 0, 1]}, 

                						{Thickness[0.005], Dashing[{0.02}], 

                  							RGBColor[1, 0, 0]}}, 
				
            DisplayFunction -> Identity]; 
		
        Show[scatter, plot1, 
			
          DisplayFunction -> $DisplayFunction])\)\)

Here is an example of the output produced by this program. Students can repeatedly execute the program which further illustrates the sample variation described earlier.

least[A,B,sigma,indep]

A visual examination of the least squares estimates

At this point in the session, students have observed that the least squares line varies from sample to sample about the regression line. As a result, they understand that the least squares estimates, a and b, of the parameters A and B in the linear model are statistics and vary from sample to sample. How well do these estimators perform? Is there an association between them?

In this portion of the laboratory session, we begin to answer these questions by generating 1000 samples of observations from the model. For each sample, we extract the corresponding estimates (a and b) for the intercept and slope (A and B) of the regression line , and generate 1000 points of the form (a, b). These 1000 points are then plotted together with the horizontal and vertical axes set to intersect at the point (A, B).

\(sims[A_, B_, sigma_, indep_List] := 
	
      \((s = Table[{indep[[i]], 

              					
              Random[NormalDistribution[
                  A + B\ indep[[i]], 
                  sigma]]}, 
				{j, 1000}, {i, Length[indep]}]; 
		
        t = \((Fit[#1, {1, x}, x]&)\)/@s; 
		
        c = \((CoefficientList[#1, x]&)\)/@t; 
		
        ListPlot[c, AxesOrigin -> {A, B}, 
			
          AxesLabel -> {"a", "b"}])\)\)

sims[A,B,sigma,indep]

What elements of the plot stand out?  We see a dramatic display of the sample variability of the least squares estimates and at the same time we observe a clear pattern to this variation.

The cloud of points appears to be centered at the point (A, B). (Students can reexecute the previous command to verify that this is no fluke!). What does this suggest? Students are led to the hypothesis that a and b are unbiased estimators of  A and B respectively. The sample means of these estimates are easily computed by Mathematica as further support for this hypothesis.

Map[Mean,Transpose[c]]

This leads students, naturally, to a greater conceptual understanding of unbiasedness.The mathematical verification that a and b are, respectively, unbiased estimators of A and B is now a more meaningful and interesting exercise .

Students also observe dramatic visual evidence that  the estimators a and b are negatively correlated. Are a and b always negatively correlated? It is a good intuitive exercise to consider whether or not this is true. There is now strong  motivation for students to derive the covariance of a and b:

Cov(a,  b) = \(\(\(-σ\^2\) 
            x\&_\)\/\(∑\+\(i = 1\)\%n\((x\_i - x\&_)\)\^2\)\)

Now, the answer to the question is clear. Examples such as the following illustrate this in the  case where \(TraditionalForm\`x\&_\)< 0 and in the case where \(TraditionalForm\`x\&_\)= 0. Of course, the centering of the estimates at the point (A, B) is still evident.

sims[A, B, sigma, -indep]

sims[A, B, 
      sigma, {-5, -4, -3, -2, -1, 1, 2, 3, 4, 5}]

Conclusion

The incorporation of computer laboratory assignments can add an interesting and dynamically effective pedagogical element to the calculus-based mathematical statistics course.  Mathematica simulations and graphical representations of them provide students with an experiential connection to random phenomena that explains, suggests, and motivates the statistical theory that is studied at the undergraduate level. The present example is designed to include all three elements. Such laboratory activities can result in increased insight and understanding that enhance the educational experience of students.

References

1. Larsen, R.J. and Marx, M.L. (1986). An Introduction to Mathematical Statistics and its Applications, 2nd edition. Prentice-Hall, Englewood Cliffs, New Jersey.



2. Rice, John A. (1995). Mathematical Statistics and Data Analysis, 2nd edition. Wadsworth, Belmont, California.

About the author

Evan Fisher is an Associate Professor in the 
Department of Mathematics at Lafayette College. He uses Mathematica frequently 
in his teaching of probability and statistics. He received his Ph.D. in 
mathematics from the University of Illinois at Urbana-Champaign. His area of 
research is in probability.