#set document(title: "6.3 Machine Learning in Regression Analysis", author: "OpenStax / XYZ Homework")
#set page(width: 8.5in, height: auto, margin: 1in)
#import "@preview/cetz:0.5.2"
#set text(font: ("STIX Two Text", "Libertinus Serif", "New Computer Modern"), size: 10.5pt, lang: "en")
#show math.equation: set text(font: ("STIX Two Math", "New Computer Modern Math"))
#set par(justify: true, leading: 0.62em, spacing: 0.9em)
#set enum(spacing: 1.1em)   // room between list items so tall inline fractions don't collide
#set list(spacing: 1.1em)
#set table(stroke: 0.5pt + rgb("#c7ccd3"))

#let BLUE = rgb("#183B6F")     // brand navy — section bars + example/solution labels (white on navy 11.09:1)
#let ORANGE = rgb("#A94509")   // brand primary-700 — AA-safe deep orange for TEXT (5.93:1 on white; raw brand #F37021 is 2.94:1 and must never carry text)
#let RED = rgb("#DC2626")      // brand error-600
#let GREEN = rgb("#059669")    // brand success-600 (decoration only; small green text uses green-text #007942)

#show heading.where(level: 1): it => block(width: 100%, above: 0pt, below: 16pt,
  fill: gradient.linear(BLUE, rgb("#2C5AA0")), inset: (x: 14pt, y: 12pt), radius: 3pt,
  text(fill: white, weight: "bold", size: 19pt, it.body))
#show heading.where(level: 2): it => block(width: 100%, above: 18pt, below: 10pt,
  fill: BLUE, inset: (x: 10pt, y: 6pt), radius: 2pt,
  text(fill: white, weight: "bold", size: 12pt, it.body))
#show heading.where(level: 3): it => text(fill: ORANGE, weight: "bold", size: 12.5pt, it.body)
#show heading.where(level: 4): it => text(fill: BLUE, weight: "bold", size: 10.5pt, it.body)

#let examplebox(label, title, body) = block(width: 100%, breakable: true,
  fill: rgb("#EFF1F5"), stroke: 0.5pt + rgb("#CFDDF0"), radius: 4pt, inset: 10pt,
  above: 12pt, below: 12pt)[
  #block(below: 6pt)[#box(fill: BLUE, inset: (x: 6pt, y: 2pt), radius: 2pt,
    text(fill: white, weight: "bold", size: 8.5pt, label)) #h(0.4em) #strong[#title]]
  #body]

// rail = decorative left rule (raw brand token); labelcolor = AA-safe label text shade
#let notebox(label, rail, labelcolor, tint, body) = block(width: 100%, breakable: true, fill: tint,
    stroke: (left: 3pt + rail), inset: (left: 10pt, rest: 8pt), radius: (right: 4pt),
    above: 11pt, below: 11pt)[
    #text(fill: labelcolor, weight: "bold", size: 7.5pt, tracking: 0.5pt)[#upper(label)]
    #linebreak() #body]

#let solutionbox(body) = block(above: 4pt, below: 8pt)[
  #text(fill: BLUE, weight: "bold", size: 8.5pt)[Solution] #linebreak() #body]

#let figph(msg) = block(width: 100%, height: 60pt, fill: rgb("#f6f7f9"),
  stroke: (paint: rgb("#c7ccd3"), dash: "dashed"), radius: 4pt, inset: 10pt)[
  #align(center + horizon, text(fill: rgb("#889"), style: "italic", size: 9pt, msg))]

// Standardize inlined figure sizes: measure the natural CeTZ canvas, then scale to a
// consistent envelope (aspect-aware; see build_typst.py FIG_* constants). Unlike the
// print preamble, dimensions are FLOORED: in an editor a user can trim a figure to a
// degenerate 1-D shape (a bare line), and w/h or tw/w would then divide by zero.
#let _STD_W = 3.5
#let _WIDE_W = 5.6
#let _MAX_H = 3.4
#let _ASPECT_WIDE = 2.2
#let _UPSCALE_MAX = 1.15
#let stdfig(body) = context {
  let m = measure(body)
  let w = calc.max(m.width / 1in, 0.01)
  let h = calc.max(m.height / 1in, 0.01)
  let tw = if w / h > _ASPECT_WIDE { _WIDE_W } else { _STD_W }
  let s = calc.min(tw / w, _MAX_H / h, _UPSCALE_MAX)
  align(center, box(scale(x: s * 100%, y: s * 100%, reflow: true, body)))
}

#show figure: set block(breakable: false)
#set figure(gap: 8pt)
#show figure.caption: set text(size: 8.5pt, fill: rgb("#555"))


  
== 6.3#h(0.6em)Machine Learning in Regression Analysis


  
    
=== Learning Outcomes


    By the end of this section, you should be able to:


    - Use bootstrapping to analyze variation in linear regression.
- Outline assumptions relevant to a multiple linear regression.
- Perform multiple linear regressions and analyze significance of coefficients.


  
  Regression is a term that applies to many different techniques in data analysis and machine learning. The main idea is to model a relationship between one or more input (or independent) variables #math.equation(block: false, alt: "X")[$X$] and an output (or dependent) variable #emph[y] by creating an explicit functional equation #math.equation(block: false, alt: "y equals f open parenthesis X close parenthesis")[$y = f ( X )$]. Once the equation is found, new inputs can be plugged into it and the results would serve as predictions. The types of regressions are based on the form of the equations that are generated: linear, polynomial, logistic, exponential, etc. Most of these types of regression analyses may be done on single inputs or using multiple input variables. Some of these topics have already been discussed in this textbook, and so in this chapter, we delve into a more advanced analysis. First we explore #emph[bootstrapping] as a way to find more information about the reliability and variability of the parameters of a linear regression. Then we discuss multiple linear and logistic regressions, including how to perform these tasks in Python.


  
    
=== Linear Regression and Bootstrapping


    The fundamentals of linear regression were presented in Inferential Statistics and Regression Analysis. While linear regression is a very powerful tool, it may suffer from inaccuracies due to noisy data. How confident can we be in the parameters of a linear model? In this section we use a resampling technique called #strong[bootstrapping], which is aimed at estimating variation to build a linear regression model with parameters within a specified confidence interval. Recall that we defined bootstrapping in Bootstrapping Methods as a technique in which repeated samples are taken hundreds or thousands of times and then the sample mean is calculated for each sample.


    
      
==== Bootstrapping Fundamentals


      Recall, for a set of bivariate data—that is, a set of points   #math.equation(block: false, alt: "X equals { open parenthesis x sub i , y sub i close parenthesis }")[$X = \{ ( x_(i) , y_(i) ) \}$] in the plane—the linear regression or line of best fit is a line with equation #math.equation(block: false, alt: "y equals a plus b x")[$y = a + b x$], in which the values of #math.equation(block: false, alt: "a")[$a$]  and #math.equation(block: false, alt: "b")[$b$] are chosen to minimize the squared vertical distance from each #math.equation(block: false, alt: "y sub i")[$y_(i)$] to the predicted value, #math.equation(block: false, alt: "y ^ sub i equals a plus b x sub i")[$hat(y)_(i) = a + b x_(i)$]. The values of #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$] clearly depend on the points of the dataset #math.equation(block: false, alt: "X")[$X$], and so these parameters are just our best guess estimates for the “true” values of #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$]. Perhaps a different sampling of #math.equation(block: false, alt: "X")[$X$] would yield slightly different values of #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$]. We may not have the luxury of resampling from the population, and so we do the next best thing: Resample from the set #math.equation(block: false, alt: "X")[$X$] itself! This process is called bootstrapping because you are attempting to “pull yourself up by your own bootstraps,” as the old saying goes. Here is an outline of the process.


      + Choose a large number #math.equation(block: false, alt: "N")[$N$] of samples taken from the dataset #math.equation(block: false, alt: "X")[$X$], with replacement (so repetitions of data are allowed).
+ Find the linear regression for this set of #math.equation(block: false, alt: "N")[$N$] samples. Record the #math.equation(block: false, alt: "a")[$a$] (intercept) and #math.equation(block: false, alt: "b")[$b$] (slope) values into lists #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$], respectively.
+ Repeat steps 1 and 2 many times.
+ The sets #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$] represent distributions of possible values of the parameters #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b")[$b$] for the true linear fit, #math.equation(block: false, alt: "y equals a plus b x")[$y = a + b x$].
+ Let #math.equation(block: false, alt: "a minus equals M E A N open parenthesis A close parenthesis")[$limits(a)^(−) = M E A N ( A )$] and #math.equation(block: false, alt: "b minus equals M E A N open parenthesis B close parenthesis")[$limits(b)^(−) = M E A N ( B )$]. These are the parameters we will chose for our (bootstrapped) linear regression, #math.equation(block: false, alt: "y equals a minus plus b minus x")[$y = limits(a)^(−) + limits(b)^(−) x$].
+ Find the standard deviation of #math.equation(block: false, alt: "A")[$A$] and #math.equation(block: false, alt: "B")[$B$] as well in order to produce a confidence interval for each parameter.


      These steps will next be illustrated with the help of Python.


    
    
      
==== Linear Regression with Bootstrapping in Python


      In NCAA Division I basketball, a statistic called BARTHAG provides a measure of many different factors that affect a team’s chances of winning against an average team in the same division. This statistic (along with many others) can be found in the dataset #link("https://openstax.org/r/NCAA")[#strong[NCAA-2021-stats.csv]]. Let #math.equation(block: false, alt: "X")[$X$] be the dataset consisting of pairs (BARTHAG, W), where W stands for the number of wins in 2021. There are 346 data points in this set. Here is the code that produces the scatterplot and regression line.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      import pandas as pd  \#\# for dataset management
      import matplotlib.pyplot as plt \#\# for data visualization
      from sklearn.linear\_model import LinearRegression
      
      \# Read data
      data = pd.read\_csv('NCAA-2021-stats.csv')
      data = data\[\['BARTHAG', 'W'\]\].dropna()
      
      \# Find the linear regression model
      model = LinearRegression()
      x = data\['BARTHAG'\].values.reshape(-1,1)
      y = data\['W'\].values.reshape(-1,1)
      model.fit(x,y)
      y\_pred = model.predict(x)
      
      \# Visualizing the linear regression on the scatterplot
      plt.scatter(x,y)
      plt.plot(x, y\_pred, color='red')
      plt.show()
      
      a = model.intercept\_
      b = model.coef\_
      print("Linear model: y = %1.3f + %2.3f x" % (a, b))
    
        The resulting output will look like this:


        #figure(figph[A scatterplot with 346 blue data points. The X axis ranges from 0 to 1 at increments of 0.2. The Y axis ranges from 0 to 25 at increments of 5. A red regression line increases from left to right from about 5 to 20.], alt: "A scatterplot with 346 blue data points. The X axis ranges from 0 to 1 at increments of 0.2. The Y axis ranges from 0 to 25 at increments of 5. A red regression line increases from left to right from about 5 to 20.", caption: none)


        
    Linear model: y = 4.820 + 14.682 x
]


      Then, the set is resampled 50 times, using 100 data points each time. A linear regression is run each time, producing a set of 50 values of #math.equation(block: false, alt: "a")[$a$] and 50 values of #math.equation(block: false, alt: "b")[$b$]. All 50 regression lines are shown on the same graph.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
from sklearn.utils import resample

      data.plot.scatter(x='BARTHAG', y='W')
      
      \# Build arrays for slopes and intercepts
      slopes = \[\]
      intercepts = \[\]
      
      \# Resampling 100 points 50 times
      for i in range(50):
        data1 = resample(data, n\_samples=100, replace=True)
      
        x = data1\['BARTHAG'\].values.reshape(-1,1)
        y = data1\['W'\].values.reshape(-1,1)
        model.fit(x,y)
        y\_pred = model.predict(x)
      
        slopes.append(model.coef\_\[0\]\[0\])
        intercepts.append(model.intercept\_\[0\])
      
        plt.plot(x, y\_pred, color='red')
      
      plt.show()
    
        The resulting output will look like this (results for bootstrapping models may slightly vary):


        #figure(figph[A scatterplot with 346 blue data points. The X axis ranges from 0 to 1 at increments of 0.2. The Y axis ranges from 0 to 25 at increments of 5. Fifty red regression lines increasing from left to right from about 5 to 20 represent the data being resampled 50 times, using 100 data points each time.], alt: "A scatterplot with 346 blue data points. The X axis ranges from 0 to 1 at increments of 0.2. The Y axis ranges from 0 to 25 at increments of 5. Fifty red regression lines increasing from left to right from about 5 to 20 represent the data being resampled 50 times, using 100 data points each time.", caption: none)
]


      Note how all the lines are very close to one another, but there is some variation. The mean and standard deviation of intercepts and slopes may be found using the #link("https://openstax.org/r/docspython1")[#strong[statistics]] library in Python.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
import statistics as stats \# for mean and standard deviation

      a\_bar = stats.mean(intercepts)
      a\_std = stats.stdev(intercepts)
      b\_bar = stats.mean(slopes)
      b\_std = stats.stdev(slopes)
      
      print("Intercept: Mean = %1.3f, St.Dev. = %1.3f " % (a\_bar, a\_std) )
      print("Slope: Mean = %1.3f, St.Dev. = %1.3f " % (b\_bar, b\_std) )
    
        The resulting output will look like this:


        Intercept: Mean = 4.947, St.Dev. = 0.645 #linebreak() Slope: Mean = 14.489, St.Dev. = 1.267
]


      Recall that, before bootstrapping, we obtained a linear model with intercept 4.82 and slope 14.682. These are very close to the intercept and slope obtained from the bootstrapping method. The advantage to bootstrapping, though, is that we now have measures of spread (standard deviation), and so we can produce a confidence interval for both parameters. Recall that the 95% confidence interval is found at plus or minus #math.equation(block: false, alt: "open parenthesis 1.96 close parenthesis σ / the square root of N")[$( 1.96 ) σ / sqrt(N)$], where #math.equation(block: false, alt: "σ")[$σ$] is the standard deviation (See Bootstrapping Methods). This can be done in Python using the #strong[scipy.stats] package.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
import scipy.stats as st
      import numpy as np
      
      a\_min, a\_max = st.t.interval(0.95, df=len(intercepts)-1, loc=np.mean(intercepts), scale=st.sem(intercepts))
      b\_min, b\_max = st.t.interval(0.95, df=len(slopes)-1, loc=np.mean(slopes), scale=st.sem(slopes))
      
      print("95%% confidence interval for intercept: %1.2f \< a \< %1.2f " % (a\_min, a\_max) )
      print("95%% confidence interval for slope: %1.2f \< b \< %1.2f " % (b\_min, b\_max) )
    
        The resulting output will look like this:


        95% confidence interval for intercept: 4.76 \< a \< 5.13 #linebreak() 95% confidence interval for slope: 14.13 \< b \< 14.85
]


      Thus, we expect that the true value of the intercept is between 4.76 and 5.13 with 95% probability, and the true value of the slope is between 14.13 and 14.85 with 95% probability.


    
  
  
    
=== Multiple Regression Techniques


    In a simple linear regression, one variable #math.equation(block: false, alt: "open parenthesis x close parenthesis")[$( x )$] serves as the predictor for the output variable #math.equation(block: false, alt: "open parenthesis y close parenthesis")[$( y )$]. Many problems in data science, however, go well beyond a single predictor. There may be dozens of features in the data that work together to predict the output. In this context, we may use #strong[multiple regression] techniques, which use more than one input variable. In this section, we focus on two main concepts: multiple linear regression and multiple logistic regression.


    
      
==== Multiple Linear Regression


      A model for multiple linear regression involves some number of predictor variables, which we will label #math.equation(block: false, alt: "X sub 1 , X sub 2 , and so on X sub r")[$X_(1) , X_(2) , … X_(r)$]. Here, we are using capital letters to distinguish these feature variables from data values like #math.equation(block: false, alt: "x sub i")[$x_(i)$]. The multiple regression model takes the form


      #math.equation(block: true, alt: "y equals a plus b sub 1 X sub 1 plus b sub 2 X sub 2 plus and so on plus b sub r X sub r")[$y = a + b_(1) X_(1) + b_(2) X_(2) + … + b_(r) X_(r)$]
      In the equation, #math.equation(block: false, alt: "a")[$a$] is the intercept (value of #math.equation(block: false, alt: "y")[$y$] when all #math.equation(block: false, alt: "X sub k equals 0")[$X_(k) = 0$]), and the values of #math.equation(block: false, alt: "b sub k")[$b_(k)$] are the coefficients of the predictors, which function much like the slope #math.equation(block: false, alt: "b")[$b$] in the simple linear regression model. A positive #math.equation(block: false, alt: "b sub k")[$b_(k)$] indicates that the feature measured by #math.equation(block: false, alt: "X sub k")[$X_(k)$] has a positive correlation on values of #math.equation(block: false, alt: "y")[$y$], while a negative #math.equation(block: false, alt: "b sub k")[$b_(k)$] indicates a negative correlation. More sophisticated models also introduce an error term, #math.equation(block: false, alt: "E")[$E$], to account for unexplained variability in #math.equation(block: false, alt: "y")[$y$].


      The goal of multiple linear regression is to determine (or estimate) values of the parameters #math.equation(block: false, alt: "a")[$a$] and #math.equation(block: false, alt: "b sub k")[$b_(k)$] (for #math.equation(block: false, alt: "k equals 1 , 2 , and so on , r")[$k = 1 , 2 , … , r$]) that minimize the squared vertical distance between the actual and predicted outputs. There are well-known formulas for this task, and typically the regression is handled by a computer.


      There are certain assumptions that need to be checked in order to use multiple regression. The first two are about the nature of the input and output data itself. Generally, these data must be continuous.


      + The output #math.equation(block: false, alt: "y")[$y$] should be continuous, not discrete. However, regression could still be used if the output is a discrete variable with a very large, regular range of values (such as the number of people who live in a given country, which could be any positive integer).
+ The input variables, #math.equation(block: false, alt: "X sub k")[$X_(k)$], should be continuous, not discrete. However, the same comment as in (1) applies here: If the variable has a very large, regular range of possible values, then regression should work fine.


      The next two assumptions ensure that multiple regression would provide the best model. If these assumptions are not met, then a more flexible model would be required.


      + There should be a linear relationship between #math.equation(block: false, alt: "X sub k")[$X_(k)$] and for each #math.equation(block: false, alt: "k")[$k$], regarding the other variables as constant. This can be checked by running simple linear regressions on the pair #math.equation(block: false, alt: "open parenthesis X sub k , y close parenthesis")[$( X_(k) , y )$] for each #math.equation(block: false, alt: "k")[$k$] and verifying linearity.
+ The data should not have multicollinearity. Multicollinearity means that at least two of the input variables #math.equation(block: false, alt: "X sub j")[$X_(j)$], #math.equation(block: false, alt: "X sub k")[$X_(k)$] are highly correlated with one another, in which case one (or more) of the variables can be dropped.


      Finally, there are a number of technical assumptions that should be verified, including checking for outliers, making sure that the variance is more or less uniform at all parts of the domain (this is known as #emph[homoscedasticity]), and analyzing the residuals to determine if they are fairly independent and normally distributed. More advanced analysis might result in the determination that some variables are not significant in the model. These further considerations fall outside the scope of this text, however.


    
    
      
==== Multiple Linear Regression in Python


      In this section, we work through the Python code for multiple linear regression using the dataset #link("https://openstax.org/r/NCAA")[#strong[NCAA-2021-stats.csv]]. As in Machine Learning in Regression Analysis, we would like to predict the number of wins (W), but now we will include more than one predictor feature. In the first code block, we will use BARTHAG and WAB (which is known as #emph[wins above bubble], a measure of how many more wins this team has in comparison to a typical team that just makes it into the tournament on a similar schedule) to predict W.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      import pandas as pd  \#\# for dataset management
      import matplotlib.pyplot as plt \#\# for data visualization
      from sklearn.linear\_model import LinearRegression
      
      \# Read data
      data = pd.read\_csv('NCAA-2021-stats.csv')
      data = data.drop(\['SEED'\], axis=1).dropna() \#\# Remove the SEED column
      
      X = data\[\['BARTHAG', 'WAB'\]\]  \# Feature list
      y = data\[\['W'\]\]  \# Response variable
      
      \# Create the multiple linear regression model
      model = LinearRegression()
      model.fit(X, y)
      
      b = model.intercept\_\[0\]
      a1, a2 = model.coef\_\[0\]
      
      print("Intercept: ", round(b,3))
      print("Coefficients: ", \[round(a1,3), round(a2,3)\])
      print("Score: ", round(model.score(X,y),2))
    
        The resulting output will look like this:


        Intercept:  12.868 #linebreak() Coefficients:  \[4.794, 0.545\] #linebreak() Score:  0.58
]


      According to the output, the multiple regression equation for this model is:


      #math.equation(block: true, alt: "W equals 12.868 plus 4.794 times (BARTHAG) plus 0.545 times (WAB)")[$W = 12.868 + 4.794 × "(BARTHAG)" #h(0.2em) + 0.545 × "(WAB)"$]
      We can see that both BARTHAG and WAB have positive coefficients, and so as either variable increases, there is a corresponding increase in W. The R-squared score indicates that about 58% of the total variance in W is due to BARTHAG and WAB jointly. Let’s take a look at the data and see how well the points conform to the regression plane. First, here is the code to generate the scatterplot of points (BARTHAG, WAB, W) in three dimensions. (Make sure that you have already run the previous code beforehand.)


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# import libraries for graphing
      import matplotlib.pyplot as plt
      from mpl\_toolkits.mplot3d import Axes3D
      
      X1 = X\[\['BARTHAG'\]\]
      X2 = X\[\['WAB'\]\]
      
      \# Create a 3D figure
      fig = plt.figure(figsize=(10, 6))
      ax = fig.add\_subplot(111, projection='3d')
      
      \# Scatterplot of actual data points
      ax.set\_xlabel('BARTHAG')
      ax.set\_ylabel('WAB', rotation=45)
      ax.set\_zlabel('W')
      ax.set\_title('Actual Values')
      ax.set\_box\_aspect(aspect=None, zoom=0.85)
      ax.scatter(X1, X2, y, c='b', marker='o')
    
        The resulting output will look like this:


        #figure(figph[Three-dimensional scatterplot showing a cluster of blue dots representing data points. Axes are labeled “BARTHAG“ (0 to 1), “WAB” (-15 to 10), and “W” (0 to 25).], alt: "Three-dimensional scatterplot showing a cluster of blue dots representing data points. Axes are labeled “BARTHAG“ (0 to 1), “WAB” (-15 to 10), and “W” (0 to 25).", caption: none)
]


      Next, let’s take a look at the predicted values found by plugging BARTHAG and WAB into the multiple regression formula discussed previously in the chapter. We also include the regression plane. The Python code required to create the regression plane is somewhat involved, requiring another package, #link("https://openstax.org/r/numpy")[#strong[numpy]], to help create a mesh grid of #math.equation(block: false, alt: "open parenthesis x , y close parenthesis")[$( x , y )$] points, so don’t worry too much if you don’t understand that part of the code.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
import numpy as np

      \# Create a 3D figure
      fig = plt.figure(figsize=(10, 6))
      ax = fig.add\_subplot(111, projection='3d')
      
      \# Generate prediction values using the model
      y\_pred = model.predict(X)
      
      \# Create a grid of points for the regression plane
      X1min = min(X1.values)\[0\]
      X1max = max(X1.values)\[0\]
      X2min = min(X2.values)\[0\]
      X2max = max(X2.values)\[0\]
      X1\_space = np.linspace(X1min, X1max, 100)
      X2\_space = np.linspace(X2min, X2max, 100)
      x1, x2 = np.meshgrid(X1\_space, X2\_space)
      
      \# Use the regression equation to find z-values
      z = b + a1 \* x1 + a2 \* x2
      
      \# Plot the regression plane
      ax.set\_xlabel('BARTHAG')
      ax.set\_ylabel('WAB', rotation=45)
      ax.set\_zlabel('W')
      ax.set\_title('Multiple Linear Regression Plane and Predicted Values')
      ax.set\_box\_aspect(aspect=None, zoom=0.85)
      ax.scatter(X1, X2, y\_pred, c='r', marker='x')
      ax.plot\_surface(x1, x2, z, alpha=0.3, cmap='plasma')
    
        The resulting output will look like this:


        #figure(figph[Three-dimensional scatterplot showing multiple linear regression model fit. Red x's represent data points, and a tilted plane represents the regression model. Axes are labeled “BARTHAG” (0 to 1), “WAB” (-15 to 10), and “W” (2.5 to 22.5).], alt: "Three-dimensional scatterplot showing multiple linear regression model fit. Red x's represent data points, and a tilted plane represents the regression model. Axes are labeled “BARTHAG” (0 to 1), “WAB” (-15 to 10), and “W” (2.5 to 22.5).", caption: none)
]


      It may not be entirely clear from the graph, but each red “x” is a point on the regression plane. To see how well the regression formula does at modeling the actual data, we will remove the predicted values and put in the actual values together with line segments showing how close those values are to the plane.


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Create a 3D figure
      fig = plt.figure(figsize=(10, 6))
      ax = fig.add\_subplot(111, projection='3d')
      
      \# Add line segments from actual data points to predicted points
      for i in range(len(X1)):
          ax.plot(\[X1.values\[i,0\], X1.values\[i,0\]\], 
           \[X2.values\[i,0\], X2.values\[i,0\]\], 
            \[y.values\[i,0\], y\_pred\[i,0\]\], c='g', linewidth=1)
      
      \# Plot the regression plane
      ax.set\_xlabel('BARTHAG')
      ax.set\_ylabel('WAB', rotation=45)
      ax.set\_zlabel('W')
      ax.set\_title('Multiple Linear Regression')
      ax.set\_box\_aspect(aspect=None, zoom=0.85)
      ax.scatter(X1, X2, y, c='b', marker='o')
      ax.plot\_surface(x1, x2, z, alpha=0.3, cmap='plasma')
    
        The resulting output will look like this:


        #figure(figph[Three-dimensional scatterplot showing multiple linear regression model fit. Purple dots represent data points, and a tilted plane represents the regression model. Axes are labeled “BARTHAG” (0 to 1), “WAB” (-15 to 10), and “W” (0 to 25).], alt: "Three-dimensional scatterplot showing multiple linear regression model fit. Purple dots represent data points, and a tilted plane represents the regression model. Axes are labeled “BARTHAG” (0 to 1), “WAB” (-15 to 10), and “W” (0 to 25).", caption: none)
]


      Now, let’s see if we can improve the model by adding more features. Let’s add TOR (turnover percentage), 2P\_O (two-pointer shooting percentage), and 3P\_O (three-pointer shooting percentage).


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
X = data\[\['BARTHAG', 'WAB', 'TOR'\]\]  \# Feature list
      y = data\[\['W'\]\]  \# Response variable
      
      \# Create the multiple linear regression model
      model = LinearRegression()
      model.fit(X, y)
      
      b = model.intercept\_\[0\]
      a1, a2, a3 = model.coef\_\[0\]
      
      print("Intercept: ", round(b,3))
      print("Coefficients: ", \[round(a1,3), round(a2,3), round(a3,3)\])
      print("Score: ", round(model.score(X,y),2))
    
        The resulting output will look like this:


        Intercept:  16.116 #linebreak() Coefficients:  \[4.502, 0.526, -0.169\] #linebreak() Score:  0.59
]


      The formula now expands to:


      #math.equation(block: true, alt: "W equals 16.116 plus 4.502 times (BARTHAG) plus 0.526 times (WAB) minus 0.169 times (TOR)")[$W = 16.116 + 4.502 × "(BARTHAG)" #h(0.2em) + 0.526 × "(WAB)" #h(0.2em) − 0.169 × "(TOR)"$]
      The R-squared measure has increased slightly to 0.59, meaning that the new model explains a bit more of the variance. However, the small increase in prediction accuracy comes at a cost of adding complexity to the model. Simpler models tend to do better at avoiding overfitting in practice, so we will stick to the original two feature variables for this example.


    
    
      
==== Multiple Logistic Regression


      As the name suggests, multiple logistic regression is an extension of the single-variable logistic regression method to two or more input variables. The model estimates the probability of the dependent variable #math.equation(block: false, alt: "open parenthesis y close parenthesis")[$( y )$] taking on a particular category such as yes or no. Thus, the model is a binary classifier with multiple independent inputs. The standard multiple logistic regression model takes the form:


      #math.equation(block: true, alt: "logit open parenthesis y close parenthesis equals a plus b sub 1 X sub 1 plus b sub 2 X sub 2 plus and so on plus b sub r X sub r")[$"logit" #h(0.2em) ( y ) = a + b_(1) X_(1) + b_(2) X_(2) + … + b_(r) X_(r)$]
      Here, the logit function is the same as described in section 6.2.1. Multiple logistic regression works just like single logistic regression in the sense that the logit function is used to transform the data, a (multiple) linear regression is performed, and the results are transformed back by a logistic function. Note that the features may be continuous numerical data or discrete data, including Boolean (1 or 0 for True/False or Yes/No).


    
    
      
==== Multiple Logistic Regression in Python


      Revisit the dataset #link("https://openstax.org/r/CollegeCompletionData")[#strong[CollegeCompletionData.csv]] from The Confusion Matrix. It has two feature columns, “GPA” and “In-State.” Last time, we only used one feature to predict completion. It is just as easy to use multiple features!


      #notebox("Note", rgb("#8a94a6"), rgb("#556666"), rgb("#f7f8fa"))[
\# Import libraries
      import pandas as pd  \#\# for dataset management
      import matplotlib.pyplot as plt \#\# for data visualization
      import seaborn as sns \#\# for heatmap visualization
      from sklearn.linear\_model import LogisticRegression
      from sklearn.metrics import confusion\_matrix
      
      \# Read data
      data = pd.read\_csv('CollegeCompletionData.csv').dropna()
      X = data\[\['GPA','In\_State'\]\]
      y = data\['Completion'\]
      
      \# Build the logistic model
      model = LogisticRegression()
      model.fit(X,y)
      
      \# Generate model predictions
      y\_pred = model.predict(X)
      
      \# Generate the confusion matrix
      cf = confusion\_matrix(y, y\_pred)
      
      \# Plot the heatmap using seaborn and matplotlib
      sns.heatmap(cf, annot=True, fmt='d', cmap='Blues', cbar=True)
      plt.xlabel('Predicted')
      plt.ylabel('True')
      plt.yticks(rotation=0)
      plt.title('Confusion Matrix')
      plt.show()
    
        The resulting output will look like:


        #figure(figph[A confusion matrix visualizing the performance of a binary classification model. The matrix has two rows and two columns, representing the true and predicted classes, respectively. The diagonal cells (top-left and bottom-right) show the number of correctly classified instances (true positives and true negatives), while the off-diagonal cells show the number of misclassified instances (false positives and false negatives). Clockwise from top left the boxes read: 25, 5, 23, 9. A color scale on the right indicates the frequency of each value. It runs up the right side of the matrix starting as light blue at 5.0 and getting darker at the top at 25. The color of the boxes aligns with this scale.], alt: "A confusion matrix visualizing the performance of a binary classification model. The matrix has two rows and two columns, representing the true and predicted classes, respectively. The diagonal cells (top-left and bottom-right) show the number of correctly classified instances (true positives and true negatives), while the off-diagonal cells show the number of misclassified instances (false positives and false negatives). Clockwise from top left the boxes read: 25, 5, 23, 9. A color scale on the right indicates the frequency of each value. It runs up the right side of the matrix starting as light blue at 5.0 and getting darker at the top at 25. The color of the boxes aligns with this scale.", caption: none)
]


      Recall that the rows of the confusion matrix are labeled by actual values while the columns are labeled by predicted values. Thus, there are 25 true positives (TP) and 23 true negatives (TN), increasing the accuracy to #math.equation(block: false, alt: "48 / 62 equals 77 %")[$48 / 62 = 77 %$]. The number of false negatives (FN) dropped from 10 down to 5, but the number of false positives (FP) increased by 1, which is not a bad trade-off.