LogisticRegression.Rmd

---
title: "Logistic Regression"
output: 
  html_document:
    code_folding: hide
header-includes: 
  - \usepackage{amsmath}
---

<script type="text/javascript">
 function showhide(id) {
    var e = document.getElementById(id);
    e.style.display = (e.style.display == 'block') ? 'none' : 'block';
 }
  
 function openTab(evt, tabName) {
    var i, tabcontent, tablinks;
    tabcontent = document.getElementsByClassName("tabcontent");
    for (i = 0; i < tabcontent.length; i++) {
        tabcontent[i].style.display = "none";
    }
    tablinks = document.getElementsByClassName("tablinks");
    for (i = 0; i < tablinks.length; i++) {
        tablinks[i].className = tablinks[i].className.replace(" active", "");
    }
    document.getElementById(tabName).style.display = "block";
    evt.currentTarget.className += " active";
 }
</script>

```{r, include=FALSE}
library(plotly)
library(pander)
library(reshape2)
library(ResourceSelection)
library(mosaic)
library(tidyverse)
library(DT)
```


----

Regression for a qualitative binary response variable $(Y_i = 0$ or $1)$. The explanatory variables can be either quantitative or qualitative.

----

## Simple Logistic Regression Model {.tabset .tabset-pills .tabset-fade}

<div style="float:left;width:125px;" align=center>
<img src="./Images/BinomYQuantX.png" width=58px;>
</div>

Regression for a qualitative binary response variable $(Y_i = 0$ or $1)$ using a single (typically quantitative) explanatory variable.

### Overview

<div style="padding-left:125px;">

The probability that $Y_i = 1$ given the observed value of $x_i$ is called $\pi_i$ and is modeled by the equation

<div style="float:right;font-size:.8em;background-color:lightgray;padding:5px;border-radius:4px;"><a style="color:darkgray;" href="javascript:showhide('simplelogisticlatexrcode')">Math Code</a></div>


<div id="simplelogisticlatexrcode" style="display:none;">
```{}
$$
  P(Y_i = 1|\, x_i) = \frac{e^{\beta_0 + \beta_1 x_i}}{1+e^{\beta_0 + \beta_1 x_i}} = \pi_i
$$
```
</div>

<center>
<span class="tooltipr">
$P($
  <span class="tooltiprtext">The "P" stands for "Probability that..."</span>
</span><span class="tooltipr">
$Y_i$
  <span class="tooltiprtext">The response variable. The "i" denotes that this is the y-value for individual "i", where "i" is 1, 2, 3,... and so on up to $n$, the sample size.</span>
</span><span class="tooltipr">
$= 1$
  <span class="tooltiprtext">Equals 1... This states that we are assuming that the probability that the response variable $Y_i$ is a 1 for the current individual.</span>
</span><span class="tooltipr">
$| x_i)$
  <span class="tooltiprtext">Given $x_i$... in other words, the "|" says "given" and $x_i$ means the x-value of the current individual.</span>
</span><span class="tooltipr">
$=$
  <span class="tooltiprtext">Equals sign.</span>
</span><span class="tooltipr">
$\displaystyle\frac{e^{\beta_0 + \beta_1 x_i}}{1 + e^{\beta_0 + \beta_1 x_i}}$
  <span class="tooltiprtext">The logistic regression equation where $e=2.71828...$ is the "natural constant" number and $\beta_0$ is the y-intercept and $\beta_1$ is the slope.</span>
</span><span class="tooltipr">
$= \pi_i$
  <span class="tooltiprtext">The $\pi_i$ stands for the probability of individual $i$ having a y-value equal to 1 given their $x_i$ value. It is the short hand notation for $P(Y_i = 1 |x_i)$. (It is NOT the number 3.14...)</span>
</span>
</center>

<br/>







The coefficents $\beta_0$ and $\beta_1$ are difficult to interpret directly. Typicall $e^{\beta_0}$ and $e^{\beta_1}$ are interpreted instead. The value of $e^{\beta_0}$ or $e^{\beta_1}$ denotes the relative change in the odds that $Y_i=1$. The odds that $Y_i=1$ are $\frac{\pi_i}{1-\pi_i}$.


----

**Examples:** [challenger](./Analyses/Logistic Regression/Examples/challengerLogisticReg.html) | [mouse](./Analyses/Logistic Regression/Examples/mouseLogisticReg.html) 

----

</div>

### R Instructions

<div style="padding-left:125px;">

**Console** Help Command: `?glm()`

#### Perform a Logistic Regression

<a href="javascript:showhide('logistic1')">
<div class="hoverchunk">
<span class="tooltipr">
YourGlmName
  <span class="tooltiprtext">This is some name you come up with that will become the R object that stores the results of your logistic regression `glm()` command.</span>
</span><span class="tooltipr">
&nbsp;<-&nbsp;
  <span class="tooltiprtext">This is the "left arrow" assignment operator that stores the results of your `glm()` code into `YourGlmName`.</span>
</span><span class="tooltipr">
glm(
  <span class="tooltiprtext">glm( is an R function that stands for "General Linear Model". It works in a similar way that the `lm(` function works except that it requires a `family=` option to be specified at the end of the command.</span>
</span><span class="tooltipr">
Y&nbsp;
  <span class="tooltiprtext">Y is your binary response variable. It must consist of only 0's and 1's. Since TRUE's = 1's and FALSE's = 0's in R, Y could be a logical statement like (Price > 100) or (Animal == "Cat") if your Y-variable wasn't currently coded as 0's and 1's.</span>
</span><span class="tooltipr">
~&nbsp;
  <span class="tooltiprtext">The tilde symbol ~ is used to tell R that Y should be treated as a function of the explanatory variable X.</span>
</span><span class="tooltipr">
X, 
  <span class="tooltiprtext">X is the explanatory variable (typically quantitative) that will be used to explain the probability that the response variable Y is a 1.</span>
</span><span class="tooltipr">
&nbsp;data = NameOfYourDataset,  
  <span class="tooltiprtext">NameOfYourDataset is the name of the dataset that contains Y and X. In other words, one column of your dataset would be called Y and another column would be called X.</span>
</span><span class="tooltipr">
&nbsp;family=binomial)
  <span class="tooltiprtext">The family=binomial command tells the `glm(` function to perform a logistic regression. It turns out that `glm` can perform many different types of regressions, but we only study it as a tool to perform a logistic regression in this course.</span>
</span><br/><span class="tooltipr">
summary(YourGlmName)
  <span class="tooltiprtext">The `summary` command allows you to print the results of your logistic regression that were previously saved in `YourGlmName`.</span>
</span>
</div>
</a>

<div id="logistic1" style="display:none;">

Example output from a regression. Hover each piece to learn more.

<table class="rconsole">
<tr>
<td>
<span class="tooltiprout">
  Call:<br/>
  glm(formula = am ~ disp, family = binomial, data = mtcars)
    <span class="tooltiprouttext">This is simply a statement of your original glm(...) "call" that you made when performing your regression. It allows you to verify that you ran what you thought you ran in the glm(...).</span>
</span>
</td>
</tr>
</table>
      

<br/>

<table class="rconsole">
<tr>
<td colspan="2">
  <span class="tooltiprout">
    Deviance Residuals:
      <span class="tooltiprouttext">Deviance residuals are a measure of how far the fitted probability for $\pi_i$ has differed from the actual outcome of $Y_i$ in terms of the log of the fitted probability space. (This is a fairly complicated idea.)
  </span>
</td>
</tr>
<tr>
<td align="right">
  <span class="tooltiprout">
    Min<br/>
    &nbsp; -1.5651
      <span class="tooltiprouttext">"min" gives the value of the residual that is furthest below the regression line. Ideally, the magnitude of this value would be about equal to the magnitude of the largest positive residual (the max) because the hope is that the residuals are normally distributed around the line.</span>
  </span>
</td><td align="right">   
  <span class="tooltiprout">
    1Q<br/>
    &nbsp; -0.6648
      <span class="tooltiprouttext">"1Q" gives the first quartile of the residuals, which will always be negative, and ideally would be about equal in magnitude to the third quartile.</span>
  </span>
</td><td align="right"> 
  <span class="tooltiprout">
    Median<br/>
    &nbsp; -0.2460
      <span class="tooltiprouttext">"Median" gives the median of the residuals, which would ideally would be about equal to zero. Note that because the regression line is the least squares line, the mean of the residuals will ALWAYS be zero, so it is never included in the output summary. This particular median value of -0.2460 is a little smaller than zero than we would hope for and suggests a right skew in the data because the mean (0) is greater than the median (-0.2460) witnessing the residuals are right skewed. This can also be seen in the maximum being much larger in magnitude than the minimum.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    3Q<br/>
    &nbsp; 0.7276
      <span class="tooltiprouttext">"3Q" gives the third quartile of the residuals, which would ideally would be about equal in magnitude to the first quartile. In this case, it is pretty close, which helps us see that the first quartile of residuals on either side of the line is behaving fairly normally.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    Max</br>
    &nbsp; 2.2691
      <span class="tooltiprouttext">"Max" gives the maximum positive residuals, which would ideally would be about equal in magnitude to the minimum residual. In this case, it is much larger than the minimum, which helps us see that the residuals are likely right skewed.</span>
  </span>
</td>
</tr>
</table>

<br/>


<table class="rconsole">
<tr>
<td colspan="2">
  <span class="tooltiprout">
    Coefficients:
      <span class="tooltiprouttext">Notice that in your glm(...) you used only $Y$ and $X$. You did type out any coefficients, i.e., the $\beta_0$ or $\beta_1$ of the regression model. These coefficients are estimated by the glm(...) function and displayed in this part of the output along with standard errors, t-values, and p-values.</span>
  </span>
</td>
</tr><tr>
<td align="left"></td>
<td align="right">
  <span class="tooltiprout">
    &nbsp; Estimate
      <span class="tooltiprouttext">To learn more about the "Estimates" of the "Coefficients" see the "Explanation" tab, "Estimating the Model Parameters" section for details.</span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; Std. Error
      <span class="tooltiprouttext">To learn more about the "Standard Errors" of the "Coefficients" see the "Explanation" tab, "Inference for the Model Parameters" section.</span>
  </span>    
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; z value 
      <span class="tooltiprouttext">The test statistic is a regular old z-score. It is most reliable when the sample size is "large." It is a measurement of the number of standard errors the estimate is from 0.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; Pr(>|z|)
      <span class="tooltiprouttext">This is the p-value, the probability of observing a test statistic more extreme than Z. </span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    (Intercept)
      <span class="tooltiprouttext">This always says "Intercept" for any glm(...) you run in R. That is because R always assumes there is a y-intercept for your regression function.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 2.630849
      <span class="tooltiprouttext">This is the estimate of the y-intercept, $\beta_0$. It is called $b_0$. It is the value of the log of the odds that $Y_i=1$ when $x_i$ is zero. Remember to use $e^{b_0}$ to interpret this values actual effect on the odds.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 1.050170
      <span class="tooltiprouttext">This is the standard error of $b_0$. It tells you how much $b_0$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    2.505
      <span class="tooltiprouttext">The test statistic for testing the hypothesis that $\beta_0 = 0$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.01224
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_0 = 0$. It measures the probability of observing a z-score as extreme as the one observed. To compute it yourself in R, use `pnorm(-abs(your z-value))*2`.</span>
  </span>
</td><td align="left">
  <span class="tooltiprout">
    *
      <span class="tooltiprouttext">This is called a "star". One star means significant at the 0.1 level of $\alpha$.</span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    disp
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; -0.014604
      <span class="tooltiprouttext">This is the estimate of the slope, $\beta_1$. It is called $b_1$. It is the change in the log of the odds that $Y_i = 1$ as X is increased by 1 unit. Remember to use $e^{b_1}$ to compute the actual effect on the odds.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.005168
      <span class="tooltiprouttext">This is the standard error of $b_1$. It tells you how much $b_1$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    -2.826
      <span class="tooltiprouttext">This is the test statistic for testing the hypothesis that $\beta_1 = 0$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.00471
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_1 = 0$. To compute it yourself in R, use `pnorm(-abs(your z-value))*2`</span>
  </span>
</td><td align="left">
<span class="tooltiprout">
    \*\*
      <span class="tooltiprouttext">This is called a "star". Three stars means significant at the 0.001 level of $\alpha$.</span>
  </span>
</td>
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span>
  \-\-\-
  </span>
</td>
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
    Signif. codes:  0 ‘\*\*\*’ 0.001 ‘\*\*’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      <span class="tooltiprouttext">These "codes" explain what significance level the p-value is smaller than based on how many "stars" * the p-value is labeled with in the Coefficients table above.</span>
  </span>
</td>
</tr>
</table>

<br/>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
    (Dispersion parameter for binomial family taken to be 1)
      <span class="tooltiprouttext">This is a simplifying assumption of the logistic regression. Overdispersion is a common problem with logistic regression data, but is typically ignored. Unless you become an expert in statistics, this is not something you need to worry about.</span>
  </span>
</td>
</tr>
</table>

<br/>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
    &nbsp; Null Deviance:
      <span class="tooltiprouttext">The deviance of the null model. This is the model that excludes any information from the x-variable, i.e., $\beta_1=0$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;43.230
      <span class="tooltiprouttext"></span>
  </span>    
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;on 31 degrees of freedom 
      <span class="tooltiprouttext">The residual degrees of freedom. The higher this number, the more reliable the p-values will be from the logistic regression.</span>
  </span> 
</td> 
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
   Residual deviance:
      <span class="tooltiprouttext">The sum of log of the squared residuals. Essentially the resulting statistic of a goodness of fit test measuring how well the data works with the logistic regression model. Using pchisq(residual deviance, df residual deviance, lower.tail=FALSE) gives the p-value for this goodness of fit test. However, the residual deviance only follows a chi-squared distribution with df residual deviance when there are many repeated x-values, and all x-values have at least a few replicates.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;29.732
      <span class="tooltiprouttext">This can be calculated by `sum(log(myglm$res^2))`.</span>
  </span>    
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;on 30 degrees of freedom 
      <span class="tooltiprouttext">This is $n-p$ where $n$ is the sample size and $p$ is the number of parameters in the regression model. In this case, there is a sample size of 32 and two parameters, $\beta_0$ and $\beta_1$, so 32-2 = 30.</span>
  </span> 
</td> 
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
   AIC:
      <span class="tooltiprouttext">As stated in the R Help file for ?glm, "A version of Akaike's An Information Criterion..." The AIC is useful for comparing different models for the same Y-variable. The glm model with the lowest AIC (which can go negative) is the best model.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;33.732
      <span class="tooltiprouttext">The AIC for this particular model is 33.732. So if a different model (using the same Y-variable as this model) can get a lower AIC, it is a better model.</span>
  </span>    
</td> 
</tr>
</table>

<br/>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
   Number of Fisher Scoring iterations:
      <span class="tooltiprouttext">If you have taken a class in Numerical Analysis, this tells you how many iterations of the maximization algorithm were required before converging to the "Estimates" of the parameters $\beta_0$ and $\beta_1$ found in the summary.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;5
      <span class="tooltiprouttext">This implementation of glm required 5 Fisher Scoring iterations to converge. Fewer iterations hints that the model is a better fit than when many iterations are required.</span>
  </span>    
</td> 
</tr>
</table>


</div>

<br/>

#### Diagnose the Goodness-of-Fit

There are two ways to check the **goodness of fit** of a logistic regression model.

<div style="padding-left:25px;">
**Option 1**: Hosmer-Lemeshow Goodness-of-Fit Test (Most Common)

To check the **goodness of fit** of a logistic regression model where there are **few or no replicated $x$-values** use the Hosmer-Lemeshow Test.

<a href="javascript:showhide('goodnessoffit2')">
<div class="hoverchunk">
<span class="tooltipr">
library(ResourceSelection)
  <span class="tooltiprtext">This loads the ResourceSelection R package so that you can access the hoslem.test() function. You may need to run the code: install.packages("ResourceSelection") first.</span>
</span><br/><span class="tooltipr">
hoslem.test(
  <span class="tooltiprtext">This R function performs the Hosmer-Lemeshow Goodness of Fit Test. See the "Explanation" file to learn about this test.</span>
</span><span class="tooltipr">
YourGlmName
  <span class="tooltiprtext">`YourGlmName` is the name of your glm(...) code that you created previously.</span>
</span><span class="tooltipr">
\$y,&nbsp;
  <span class="tooltiprtext">ALWAYS type a "y" here. This gives you the actual binary (0,1) y-values of your logistic regression. The goodness of fit test will compare these actual values to your predicted probabilities for each value in order to see if the model is a "good fit."</span>
</span><span class="tooltipr">
YourGlmName
  <span class="tooltiprtext">`YourGlmName` is the name you used to save the results of your glm(...) code.</span>
</span><span class="tooltipr">
\$fitted,&nbsp;
  <span class="tooltiprtext">ALWAYS type "fitted" here. This gives you the fitted probabilities $\pi_i$ of your logistic regression.</span>
</span><span class="tooltipr">
g=10)
  <span class="tooltiprtext">The "g=10" is the default option for the value of g. The g is the number of groups to run the goodness of fit test on. Just leave it at 10 unless you are told to do otherwise. Ask your teacher for more information if you are interested.</span>
</span>
</div>
</a>

<div id="goodnessoffit2" style="display:none;">
```{r, echo=FALSE}
myglm <- glm(am ~ disp, data = mtcars, family = binomial)
library(ResourceSelection)
hoslem.test(myglm$y, myglm$fitted)
```

Note that the null hypothesis of the goodness-of-fit test is that "the logistic regression is a good fit." So we actually don't want to "reject the null" in this case. So a large p-value here means our logistic regression fits the data satisfactorily. A small p-value implies a poor fit and the results of the logistic regression should not be fully trusted.

</div>

<br/>

**Option 2**: Deviance Goodness-of-fit Test (Less Common)

In some cases, there are **many replicated $x$-values** for **all** x-values, i.e.,  each value of x is repeated more than 50 times. Though this is rare, it is good to use the *deviance goodness-of-fit test* whenever this happens.

<a href="javascript:showhide('goodnessoffit1')">
<div class="hoverchunk">
<span class="tooltipr">
pchisq(
  <span class="tooltiprtext">The `pchisq` command allows you to compute p-values from the chi-squared distribution.</span>
</span><span class="tooltipr">
residual deviance,&nbsp;
  <span class="tooltiprtext">The residual deviance is shown at the bottom of the output of your `summary(YourGlmName)` and should be typed in here as a number like 25.3.</span>
</span><span class="tooltipr">
df for residual deviance,&nbsp;
  <span class="tooltiprtext">The df for the residual deviance is also shown at the bottom of the output of your `summary(YourGlmName)`.</span>
</span><span class="tooltipr">
lower.tail=FALSE)
<span class="tooltiprtext">This command ensures you find the probability of the chi-squared distribution being as extreme or more extreme than the observed value of residual deviance.</span>
</span>
</div>
</a>


<div id="goodnessoffit1" style="display:none;">
```{r, echo=FALSE}
myglm <- glm(am ~ disp, data = mtcars, family = binomial)
summary(myglm)
pchisq(29.732, 30, lower.tail = FALSE)
```

The null hypothesis of the goodness-of-fit test is that the logistic regression is a good fit of the data. So a large p-value (like 0.479) is good because it allows us to trust the results of our logistic regression. When the p-value becomes very small, we must "reject the null" and conclude a poor fit, which implies that we should not trust the results of the logistic regression.

</div>

</div>
<br/>

#### Plot the Regression

<div class="tab">
  <button class="tablinks" onclick="openTab(event, 'LogisticBaseplot')">Base R</button>
  <button class="tablinks" onclick="openTab(event, 'Logisticggplot')">ggplot2</button>
</div>

<div id="LogisticBaseplot" class="tabcontent">
  <p>
<a href="javascript:showhide('LogisticBaseplot2')">
<div class="hoverchunk"><span class="tooltipr">
plot(
  <span class="tooltiprtext">The plot function allows us to draw a scatterplot where the y-axis is only 0's or 1's and the x-axis is quantitative.</span>
</span><span class="tooltipr">
Y
  <span class="tooltiprtext">Y should be some logical statement like Fail > 0 or sex == "B" or height < 68. In other words, Y needs to be a collection of 0's and 1's.</span>
</span><span class="tooltipr">
&nbsp;~&nbsp;
  <span class="tooltiprtext">The tilde is read "Y on X."</span>
</span><span class="tooltipr">
X,
  <span class="tooltiprtext">X is some quantitative variable.</span>
</span><span class="tooltipr">
&nbsp;data=
  <span class="tooltiprtext">Tell the plot which data set to use. X and Y are columns of that data set.</span>
</span><span class="tooltipr">
YourDataSet
  <span class="tooltiprtext">The name of your data set.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for plot(...) function.</span>
</span><br/><span class="tooltipr">
curve(
  <span class="tooltiprtext">A function in R that draws a "curve" on a plot. Using add=TRUE puts the curve onto the current plot.</span>
</span><span class="tooltipr">
exp(
  <span class="tooltiprtext">This allows you to compute e to the power of something. So exp(1) = e, exp(2) = e^2 and so on.</span>
</span><span class="tooltipr">
$b_0$ +
  <span class="tooltiprtext">This is the "Intercept" estimate from your logistic regression summary output.</span>
</span><span class="tooltipr">
&nbsp;$b_1$
  <span class="tooltiprtext">This is the slope estimate from your logistic regression summary output.</span>
</span><span class="tooltipr">
\*x
  <span class="tooltiprtext">The curve function demands you ALWAYS use a lower-case "x" in this function.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis.</span>
</span><span class="tooltipr">
/
  <span class="tooltiprtext">The division symbol.</span>
</span><span class="tooltipr">
(1 + 
  <span class="tooltiprtext">Required for the formula.</span>
</span><span class="tooltipr">
&nbsp;exp($b_0$ + $b_1$*x)
  <span class="tooltiprtext">This needs to match exactly the first version of this statement.</span>
</span><span class="tooltipr">
),
  <span class="tooltiprtext">Closing parenthesis for the denominator.</span>
</span><span class="tooltipr">
&nbsp;add = TRUE)
  <span class="tooltiprtext">This makes the curve be added to the current plot. If it is left out, the curve will be drawn in a new plot.</span>
</span>
</div>
</a>

<div id="LogisticBaseplot2" style="display:none;">
```{r, echo=FALSE}
myglm <- glm(am ~ disp, data = mtcars, family = binomial)
plot(am ~ disp, data = mtcars, pch=16)
curve(exp(2.631-0.015*x)/(1+exp(2.631-0.015*x)), add=TRUE)
```
</div>

  </p>
</div>

<div id="Logisticggplot" class="tabcontent">
  <p>
<a href="javascript:showhide('Logisticggplot2')">
<div class="hoverchunk"><span class="tooltipr">
ggplot(
  <span class="tooltiprtext">The ggplot(...) function is used to create a basic ggplot frame.</span>
</span><span class="tooltipr">
data=YourDataSetName,
  <span class="tooltiprtext">Use this to specify the name of your data set.</span>
</span><span class="tooltipr">
&nbsp;aes(
  <span class="tooltiprtext">The aes(...) function stands for "aesthetics" and tells the ggplot which variables to match up with the x-axis and y-axis of the graph as well as other visual things like the type of plotting characters, size of plotting characters, color, fill, and so on.</span>
</span><span class="tooltipr">
x=X,
  <span class="tooltiprtext">Use x=nameOfYourXvariable to declare the x-axis of your graph.</span>
</span><span class="tooltipr">
y=Y
  <span class="tooltiprtext">Use y=nameOfYourYvariable to declare the y-axis of your graph. If your y-variable is a logical expression, like height>60 then you must use y=as.numeric(height>60).</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for aes(...) function.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for ggplot(...) function.</span>
</span><span class="tooltipr">
&nbsp;+
  <span class="tooltiprtext">The plus sign adds a new layer to the ggplot.</span>
</span><br/><span class="tooltipr">
&nbsp; geom_point(
  <span class="tooltiprtext">Tells the plot to add the physical geometry of "points" to the plot.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for the geom_point() function.</span>
</span><span class="tooltipr">
&nbsp;+
  <span class="tooltiprtext">Add another layer to the plot.</span>
</span><br/><span class="tooltipr">
&nbsp; geom_smooth(
  <span class="tooltiprtext">Add a smoothing line to the graph.</span>
</span><span class="tooltipr">
method="glm",
  <span class="tooltiprtext">This adds a general linear model to the graph.</span>
</span><span class="tooltipr">
method.args = list(family="binomial"),
  <span class="tooltiprtext">This tells the method="glm" to choose specifically the "binomial" model, otherwise known as the "logistic regression" model.</span>
</span><span class="tooltipr">
&nbsp;se=FALSE
  <span class="tooltiprtext">Turn off the displaying of the confidence band around the logistic regression. You can turn this on if you know what it means.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for the geom_smooth() function.</span>
</span><span class="tooltipr">
&nbsp;+
  <span class="tooltiprtext">Add another layer to the ggplot.</span>
</span><br/><span class="tooltipr">
&nbsp; theme_bw()
  <span class="tooltiprtext">Give the graph a basic black and white theme. Other themes are possible, see ?theme_ in your Console.</span>
</span>
</div>
</a>

<div id="Logisticggplot2" style="display:none;">
```{r, echo=FALSE}
ggplot(mtcars, aes(disp, am)) +
  geom_point() +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  theme_bw()
```
</div>

  </p>
</div>

<br/>

#### Predict Probabilities 

To predict the probability that $Y_i=1$ for a given $x$-value, use the code

<a href="javascript:showhide('LogisticPredict')">
<div class="hoverchunk">
<span class="tooltipr">
predict(
  <span class="tooltiprtext">The predict() function allows us to use the regression model that was obtained from glm() to predict the probability that $Y_i = 1$ for a given $X_i$.</span>
</span><span class="tooltipr">
YourGlmName,&nbsp;
  <span class="tooltiprtext">`YourGlmName` is the name of the object you created when you performed your logistic regression using glm().</span>
</span><span class="tooltipr">
newdata =&nbsp;
<span class="tooltiprtext">The `newdata = ` command allows you to specify the x-values for which you want to obtain predicted probabilities that $Y_i=1$.</span>
</span><span class="tooltipr">
data.frame(xVariableName = someNumericValue),&nbsp;
<span class="tooltiprtext">The xVariableName is the same one you used in your glm(y ~ x, ...) statement for "x". Input any desired value for the "someNumericValue" spot. Then, the predict code uses the logistic regression model equation to calculate a predicted probability that $Y_i = 1$ for the given $x_i$ value that you specify.</span>
</span><span class="tooltipr">
type = "response")
  <span class="tooltiprtext">The type = "response" options specifies that you want predicted probabilities. There are other options available. See ?predict.glm for details.</span>
</span>
</div>
</a>

<div id="LogisticPredict" style="display:none;">
```{r}
myglm <- glm(am ~ disp, data = mtcars, family = binomial)
predict(myglm, newdata = data.frame(disp = 200), type = "response")
```
</div>





----

</div>

### Explanation

<div style="padding-left:125px;">


Simple Logistic Regression is used when 

* the response variable is binary $(Y_i=0$ or $1)$, and
* there is a single explanatory variable $X$ that is typically quantitative but could be qualitative (if $X$ is binary or ordinal).


#### The Model 

Since $Y_i$ is binary (can only be 0 or 1) the model focuses on describing the probability that $Y_i=1$ for a given scenario. The probability that $Y_i = 1$ given the observed value of $x_i$ is called $\pi_i$ and is modeled by the equation

$$
  P(Y_i = 1|\, x_i) = \frac{e^{\beta_0 + \beta_1 x_i}}{1+e^{\beta_0 + \beta_1 x_i}} = \pi_i
$$

The assumption is that for certain values of $X$ the probability that $Y_i=1$ is higher than for other values of $X$. 


#### Interpretation

This model for $\pi_i$ comes from modeling the *log of the odds* that $Y_i=1$ using a linear regression, i.e., 
$$
  \log\underbrace{\left(\frac{\pi_i}{1-\pi_i}\right)}_{\text{Odds for}\ Y_i=1} = \underbrace{\beta_0 + \beta_1 x_i}_{\text{linear regression}}
$$ 
Beginning to solve this equation for $\pi_i$ leads to the intermediate, but important result that
$$
    \underbrace{\frac{\pi_i}{1-\pi_i}}_{\text{Odds for}\ Y_i=1} = e^{\overbrace{\beta_0 + \beta_1 x_i}^{\text{linear regression}}} = e^{\beta_0}e^{\beta_1 x_i}
$$
Thus, while the coefficients $\beta_0$ and $\beta_1$ are difficult to interpret directly, $e^{\beta_0}$ and $e^{\beta_1}$ have a valuable interpretation. The value of $e^{\beta_0}$ is interpreted as the odds for $Y_i=1$ when $x_i = 0$. It may not be possible for a given model to have $x_i=0$, in which case $e^{\beta_0}$ has no interpretation. The value of $e^{\beta_1}$ denotes the proportional change in the odds that $Y_i=1$ for every one unit increase in $x_i$. 

Notice that solving the last equation for $\pi_i$ results in the logistic regression model presented at the beginning of this page.


#### Hypothesis Testing

Similar to linear regression, the hypothesis that
$$
  H_0: \beta_1 = 0 \\
  H_a: \beta_1 \neq 0
$$
can be tested with a logistic regression. If $\beta_1 = 0$, then there is no relationship between $x_i$ and the log of the odds that $Y_i = 1$. In other words, $x_i$ is not useful in predicting the probability that $Y_i = 1$. If $\beta_1 \neq 0$, then there is information in $x_i$ that can be utilized to predict the probability that $Y_i = 1$, i.e., the logistic regression is meaningful.


#### Checking Model Assumptions {#diagnostics}

The model assumptions are not as clear in logistic regression as they are in linear regression. For our purposes we will focus only on considering the goodness of fit of the logistic regression model. If the model appears to fit the data well, then it will be assumed to be appropriate.

<div style="padding-left:15px;">

##### Deviance Goodness of Fit Test

If there are replicated values of each $x_i$, then the deviance goodness of fit test tests the hypotheses
$$
  H_0: \pi_i = \frac{e^{\beta_0 + \beta_1 x_i}}{1+e^{\beta_0 + \beta_1 x_i}}
$$
$$
  H_a: \pi_i \neq \frac{e^{\beta_0 + \beta_1 x_i}}{1+e^{\beta_0 + \beta_1 x_i}}
$$


##### Hosmer-Lemeshow Goodness of Fit Test

If there are very few or no replicated values of each $x_i$, then the Hosmer-Lemeshow goodness of fit test can be used to test these same hypotheses. In each case, the null assumes that logistic regression is a good fit for the data while the alternative is that logistic regression is not a good fit. 

</div>

#### Prediction

One of the great uses of Logistic Regression is that it provides an estimate of the probability that $Y_i=1$ for a given value of $x_i$. This probability is often referred to as the *risk* that $Y_i=1$ for a certain individual. For example, if $Y_i=1$ implies a person has a disease, then $\pi_i=P(Y_i=1)$ represents the risk of individual $i$ having the disease based on their value of $x_i$, perhaps a measure of their cholesterol or some other predictor of the disease. 

</div>

----



###



----


## Multiple Logistic Regression Model {.tabset .tabset-pills .tabset-fade}


<div style="float:left;width:125px;" align=center>
<img src="./Images/BinomYMultX.png" width=98px;>
</div>

Logistic regression for multiple explanatory variables that can either be quantitative or qualitative or a mixture of the two.


### Overview

<div style="padding-left:125px;">

Select a model to see interpretation details, an example, and R Code help.

<div class="tab">
  <button class="tablinks" onclick="openTab(event, 'LearnMoresimpleLogisticModel')">Simple</button>
  <button class="tablinks" onclick="openTab(event, 'LearnMoreQuadraticModel')">Quadratic</button>
  <button class="tablinks" onclick="openTab(event, 'LearnMoreCubicModel')">Cubic</button>
  <button class="tablinks" onclick="openTab(event, 'LearnMoreTwoLinesModel')">Two-Lines</button>
  <button class="tablinks" onclick="openTab(event, 'LearnMorethreeDModel')">3D</button>
  <button class="tablinks" onclick="openTab(event, 'LearnMoreHDModel')">HD</button>
</div>

<div id="LearnMoresimpleLogisticModel" class="tabcontent" style="display:block;">
  <p>


<table>
  <tr>
  <td>

```{r, fig.width=1.5, fig.height=1.5, echo=FALSE}
par(mai=c(0.1,0.1,0.3,0.1))
plot(Wind > 10 ~ Temp, data=airquality, pch=21, cex=1.5, bg=rgb(.1,.1,.1,.1), col="skyblue", yaxt='n', xaxt='n', axes=FALSE)
glm.slr <- glm(Wind > 10 ~ Temp, data=airquality, family=binomial)
curve(exp(glm.slr$coef[1] + glm.slr$coef[2]*x)/(1+exp(glm.slr$coef[1] + glm.slr$coef[2]*x)), col="orange", lwd=2, add=TRUE)
```

  </td><td style="text-align: center;padding-left:15px;">

$$
  P(Y_i = 1|\, x_i) = \frac{e^{\beta_0 + \beta_1 x_i}}{1+e^{\beta_0 + \beta_1 x_i}} = \pi_i
$$

  </td>
  </tr>
</table>

<br/>

The Simple Logistic Regression model uses a single x-variable once: $X_i$.

| Parameter | Effect |
|-----------|--------|
| $\beta_0$ | Y-intercept of the Model. Only interpreted by computing $e^{\beta_0}$, which gives the "baseline odds" of the model. Technically only meaningful for when $x_i=0$ is a reasonable value. |
| $\beta_1$ | Slope term of the model. Interpreted by computing $e^{\beta_1}$, which gives the multiplicative change in the odds for each 1 unit increase in x. Say, $e^{\beta_1} = 1.2$, then the odds are 1.2 times what they were before whenever x is increased by 1 unit. That's a twenty percent increase in odds. |




</p>
</div>






<div id="LearnMoreQuadraticModel" class="tabcontent">
  <p>
  
<table>
  <tr>
  <td>

```{r, fig.width=1.5, fig.height=1.5, echo=FALSE}
par(mai=c(0.1,0.1,0.3,0.1))
plot(ccf > 70 ~ month, data=Utilities, pch=21, bg=rgb(.1,.1,.1,.1), col="skyblue", yaxt='n', xaxt='n', main="", cex.main=1, cex=1.5, axes=FALSE)
glm.quad <- glm(ccf > 70 ~ month + I(month^2), data=Utilities, family=binomial)
b <- coef(glm.quad)
curve(exp(b[1] + b[2]*x + b[3]*x^2)/(1+exp(b[1] + b[2]*x + b[3]*x^2)), col="orange", lwd=2, add=TRUE)
```

  </td><td style="text-align: center;padding-left:15px;">

$$
 P(Y_i = 1| X_i) = \frac{e^{\beta_0 + \beta_1 X_i + \beta_2 X_i^2}}{1 + e^{\beta_0 + \beta_1 X_{i } + \beta_2 X_{i}^2}}
$$

  </td>
  </tr>
</table>

<br/>


The Quadratic Logistic Regression model uses the same $X$-variable twice, once with a $\beta_1 X_i$ *slope term* and once with a $\beta_2 X_i^2$ *quadratic term*. 

| Parameter | Effect |
|-----------|-------------------------------------------------------------------------|
| $\beta_0$ | Y-intercept of the Model. Interpreted by computing $e^{\beta_0}$, which gives the baseline odds of a success. Only interpreteable when $X_i=0$ is meaningful. |
| $\beta_1$ | Controls the x-position of the center point of the "vertex" of the quadratic logistic model by $\frac{-\beta_1}{2\cdot\beta_2}$. Not directly interpretable, even with $e^{\beta_1}$, as the effect of $\beta_1$ is not independent of $\beta_2$. |
| $\beta_2$ | Controls the concavity and "steepness" of the Model: negative values result in a logistic model with a maximum point on the curve, positive values result in a minimum point on the curve; large values imply "steeper" curves and low values imply "flatter" curves. Also involved in the position of the vertex, see $\beta_1$'s explanation. Not directly interpretable. |





**An Example**

Using the `airquality` data set, we run the following "quadratic" logistic regression. Pay careful attention to how the mathematical model for $P(Y_i=1 | X_i) = \ldots$ is translated to R-Code inside of `glm(...)` by using the log(Odds) model instead of the $P(Y_i = 1 | X_i)$ model.

$$
  \underbrace{\log\left(\overbrace{\frac{\pi_i}{1-\pi_i}}^{\text{Odds}\ Y_i = 1}\right)}_\text{Temp>80} \underbrace{=}_{\sim} \overbrace{\beta_0}^{\text{y-int}} + \overbrace{\beta_1}^{\stackrel{\text{slope}}{\text{term}}} \underbrace{X_{i}}_\text{Month} \underbrace{+}_{+} \overbrace{\beta_2}^{\stackrel{\text{quadratic}}{\text{term}}}  \underbrace{X_{i}^2}_\text{I(Month^2)}
$$
Then, by using the statement in `glm(..., family=binomial)` the model that is put into the glm is translated back to the logistic model for you. Note that whether we look at the log(Odds) model or the $P(Y_i=1)$ model, the values of each $\beta$ are the same.



<a href="javascript:showhide('quadraticregressionexamplecode')">
<div class="hoverchunk">
<span class="tooltipr">
glm.quad <- 
  <span class="tooltiprtext">A name we made up for our "quadratic" logistic regression.</span>
</span><span class="tooltipr">
glm(
  <span class="tooltiprtext">R function glm used to perform generalized linear regressions, among which the logistic regression model is a specific example. But glm can do much more. </span>
</span><span class="tooltipr">
Temp > 75
  <span class="tooltiprtext">Y-variable, should be 0's and 1's. In this case the `Temp > 75` statement is translated into TRUE (1) and FALSE (0) values. It may be preferred to use a `mutate(y = ifelse(Temp > 75, 1, 0))` statement prior to running the glm to create a column of 0's and 1's.</span>
</span><span class="tooltipr">
&nbsp;~&nbsp;
<span class="tooltiprtext">The tilde `~` is what glm(...) uses to state the regression equation $Y_i = ...$ for the log of the odds linear equation. Notice that the `~` is not followed by $\beta_0 + \beta_1$ like $Y_i = ...$. Instead, $X_{i}$ (Month in this case) is the first term following `~`. This is because the $\beta$'s are going to be estimated by the glm(...). These "Estimates" can be found using summary(glmObject) and looking at the **Estimates** column in the output.</span>
</span><span class="tooltipr">
Month
  <span class="tooltiprtext">$X_{i}$, should be quantitative.</span>
</span><span class="tooltipr">
&nbsp;+&nbsp;
  <span class="tooltiprtext">The plus `+` is used between each term in the model. Note that only the x-variables are included in the glm(...) from the $\text{log of the odds} = ...$ model. No beta's are included.</span>
</span><span class="tooltipr">
I(Month^2)
  <span class="tooltiprtext">$X_{i}^2$, where the function I(...) protects the squaring of Month from how glm(...) would otherwise interpret that statement. The I(...) function must be used anytime you raise an x-variable to a power in the glm(...) statement.</span>
</span><span class="tooltipr">
, data=airquality
  <span class="tooltiprtext">This is the data set we are using for the regression.</span>
</span><span class="tooltipr">
, family=binomial
  <span class="tooltiprtext">This declares that a logistic regression will be performed instead of a linear regression.</span>
</span><span class="tooltipr">
)  
  <span class="tooltiprtext">Closing parenthsis for the glm(...) function.</span>
</span><span class="tooltipr">
&nbsp;&nbsp;&nbsp;&nbsp;  
  <span class="tooltiprtext">Press Enter to run the code.</span>
</span><span class="tooltipr" style="float:right;">
&nbsp;...&nbsp; 
  <span class="tooltiprtext">Click to View Output.</span>
</span>
</div>
</a>
<div id="quadraticregressionexamplecode" style="display:none;">
Pay special attention to how the glm(...) code uses $Y_i \sim X_{i} + X_{i}^2$ and drops all $\beta$'s from the model statement. This is because the estimates for the $\beta$'s are given by the output of the glm(...) funtion in the "Estimates" column of summary(....).
</div>


```{r}
glm.quad <- glm(Temp > 75 ~ Month + I(Month^2), data=airquality, family=binomial)
emphasize.strong.cols(1)
pander(summary(glm.quad)$coefficients, )
```


The **estimates** shown in the summary output table above approximate the $\beta$'s in the log of the odds logistic regression model: 

* $\beta_0$ is estimated by the (Intercept) value of -47.99, 
* $\beta_1$ is estimated by the `Month` value of 13.87, and
* $\beta_2$ is estimated by the `I(Month^2)` value of -0.9452.

Because the estimate of the $\beta_2$ term is negative (-0.9452), this parabola will "open down" (concave). This tells us that probabilitiy of the high temperature exceeding 75 degrees F will increase to a point, then decrease again. The vertex of this parabola will be at $-b_1/(2b_2) = -(13.87)/(2\cdot (-0.9452)) = 7.337072$ months, which tells us that the highest average temperature will occur around mid July (7.34 months to be exact). The y-intercept is -47.99, which would be an odds of nearly zero, $e^{-47.99}\approx 1.44e-21$, that the temperature would exceed 75 degrees F if it were possible for the month to be "month zero." Since this is not possible, the y-intercept is not meaningful for this model.

Note that interpreting either $\beta_1$ or $\beta_2$ by themselves is quite difficult because they both work with together with $X_{i}$. 

$$
\log\left(\overbrace{\frac{\pi_i}{1-\pi_i}}^\text{odds of Temp >75}\right) \approx \overbrace{-47.99}^\text{y-int} + \overbrace{13.87}^{\stackrel{\text{slope}}{\text{term}}} X_{i} + \overbrace{-0.9452}^{\stackrel{\text{quadratic}}{\text{term}}} X_{i}^2
$$

The regression function is drawn as follows. Be sure to look at the "Code" to understand how this graph was created using the ideas in the equation above.

<table>
<tr>
<td>

**Using Base R**

```{r}
plot(Temp>75 ~ Month, data=airquality, col="skyblue", pch=21, bg="gray83", main="Quadratic Model using airquality data set", cex.main=1)

#get the "Estimates" automatically:
b <- coef(glm.quad)
# Then b will have 3 numbers stored inside:
# b[1] is the estimate of beta_0: -47.99
# b[2] is the estimate of beta_1: 13.87
# b[3] is the estimate of beta_2: -0.9452
curve(exp(b[1] + b[2]*x + b[3]*x^2)/(1+exp(b[1] + b[2]*x + b[3]*x^2)), col="skyblue", lwd=2, add=TRUE)
```
  
</td>
<td>

**Using ggplot2**

```{r}
#get the "Estimates" automatically:
b <- coef(glm.quad)
# Then b will have 3 estimates:
# b[1] is the estimate of beta_0: -47.99
# b[2] is the estimate of beta_1: 13.87
# b[3] is the estimate of beta_2: -0.9452

ggplot(airquality, aes(y=ifelse(Temp>75,1,0), x=Month)) +
  geom_point(pch=21, bg="gray83", color="skyblue", alpha=0.1, size=6) +
  #geom_smooth(method="lm", se=F, formula = y ~ poly(x, 2)) + #easy way, but the more involved manual way using stat_function (see below) is more dynamic.
  stat_function(fun = function(x) exp(b[1] + b[2]*x + b[3]*x^2)/(1+exp(b[1] + b[2]*x + b[3]*x^2)), color="skyblue") +
  labs(title="Quadratic Model using airquality data set") 
  

```

</td>
</tr>
</table>


</p>
</div>






<div id="LearnMoreCubicModel" class="tabcontent" style="display:none;">
  <p>


<table>
  <tr>
  <td>

```{r, fig.width=1.5, fig.height=1.5, echo=FALSE}
par(mai=c(0.1,0.1,0.3,0.1))
B2 <- filter(Births78, day_of_week %in% 3:6, day_of_year >50, day_of_year < 320)
plot(births ~ day_of_year, data=B2, pch=21, bg="gray83", col="skyblue", yaxt='n', xaxt='n', axes=FALSE)
lm.cubic <- lm(births ~ day_of_year + I(day_of_year^2) + I(day_of_year^3), data=B2)
b <- coef(lm.cubic)
curve(b[1] + b[2]*x + b[3]*x^2 + b[4]*x^3, col="orange", lwd=2, add=TRUE)
```

  </td><td style="text-align: center;padding-left:15px;">

$$
 Y_i = \overbrace{\underbrace{\beta_0 + \beta_1 X_i + \beta_2 X_i^2 + \beta_3 X_i^3}_{E\{Y_i\}}}^\text{Cubic Model} + \epsilon_i
$$

  </td>
  </tr>
</table>


<br/>

The Cubic model uses the same $X$-variable thrice, once with a $\beta_1 X_i$ term, once with a $\beta_2 X_i^2$ term, and once with a $\beta_3 X_i^3$ term. The $X_i^3$ term is called the "cubic" term.

| Parameter | Effect |
|-----------|-------------------------------------------------------------------------|
| $\beta_0$ | Y-intercept of the Model. |
| $\beta_1$ | No clear interpretation, but could be called the "base slope coefficient" and contributes to the position of the inflection points of the cubic function. |
| $\beta_2$ | No clear interpretation, but it also contributes to the location of the inflection points. |
| $\beta_3$ | This is the coefficient of the cubic term. No clear interpretation, but it determines the concavity of the model by its sign. |



**An Example**

Using the `CO2` data set, we run the following "cubic" regression.

$$
  \underbrace{Y_i}_\text{uptake} \underbrace{=}_{\sim} \overbrace{\beta_0}^{\text{y-int}} + \overbrace{\beta_1}^{\stackrel{\text{slope}}{\text{term}}} \underbrace{X_{i}}_\text{conc} + \overbrace{\beta_2}^{\stackrel{\text{quadratic}}{\text{term}}}  \underbrace{X_{i}^2}_\text{I(conc^2)} +  \overbrace{\beta_3}^{\stackrel{\text{cubic}}{\text{term}}}  \underbrace{X_{i}^3}_\text{I(conc^3)} + \epsilon_i
$$


<a href="javascript:showhide('cubicregressionexamplecode')">
<div class="hoverchunk">
<span class="tooltipr">
lm.cubic <- 
  <span class="tooltiprtext">A name we made up for our "cubic" regression.</span>
</span><span class="tooltipr">
lm(
  <span class="tooltiprtext">R function lm used to perform linear regressions in R. The lm stands for "linear model".</span>
</span><span class="tooltipr">
uptake
  <span class="tooltiprtext">Y-variable, should be quantitative.</span>
</span><span class="tooltipr">
&nbsp;~&nbsp;
<span class="tooltiprtext">The tilde `~` is what lm(...) uses to state the regression equation $Y_i = ...$. Notice that the `~` is not followed by $\beta_0 + \beta_1$ like $Y_i = ...$. Instead, $X_i$ is the first term following `~`. This is because the $\beta$'s are going to be estimated by the lm(...). These estimates can be found using summary(lmObject).</span>
</span><span class="tooltipr">
conc
  <span class="tooltiprtext">$X_{i}$, should be quantitative.</span>
</span><span class="tooltipr">
&nbsp;+&nbsp;
  <span class="tooltiprtext">The plus `+` is used between each term in the model. Note that only the x-variables are included in the lm(...) from the $Y_i = ...$ model. No beta's are included.</span>
</span><span class="tooltipr">
I(conc^2)
  <span class="tooltiprtext">$X_{i}^2$, where the function I(...) protects the squaring of conc from how lm(...) would otherwise interpret that statement. The I(...) function must be used anytime you raise an x-variable to a power in the lm(...) statement.</span>
</span><span class="tooltipr">
&nbsp;+&nbsp;
  <span class="tooltiprtext">The plus `+` is used between each term in the model. Note that only the x-variables are included in the lm(...) from the $Y_i = ...$ model. No beta's are included.</span>
</span><span class="tooltipr">
I(conc^3)
  <span class="tooltiprtext">$X_{i}^3$, where the function I(...) protects the cubing of conc from how lm(...) would otherwise interpret that statement. The I(...) function must be used anytime you raise an x-variable to a power in the lm(...) statement.</span>
</span><span class="tooltipr">
, data=CO2
  <span class="tooltiprtext">This is the data set we are using for the regression.</span>
</span><span class="tooltipr">
)  
  <span class="tooltiprtext">Closing parenthsis for the lm(...) function.</span>
</span><span class="tooltipr">
&nbsp;&nbsp;&nbsp;&nbsp;  
  <span class="tooltiprtext">Press Enter to run the code.</span>
</span><span class="tooltipr" style="float:right;">
&nbsp;...&nbsp; 
  <span class="tooltiprtext">Click to View Output.</span>
</span>
</div>
</a>
<div id="cubicregressionexamplecode" style="display:none;">
Pay special attention to how the lm(...) code uses $Y_i \sim X_{i} + X_{i}^2$ and drops all $\beta$'s and $\epsilon$ from the model statement. This is because the estimates for the $\beta$'s and $\epsilon$ are given by the output of the lm(...) funtion in the "Estimates" column of summary(....) and in `lmObject$residuals`.
</div>


```{r}
lm.cubic <- lm(uptake ~ conc + I(conc^2) + I(conc^3), data=CO2)
pander(summary(lm.cubic)$coefficients)
```



The **estimates** shown above approximate the $\beta$'s in the regression model: $\beta_0$ is estimated by the (Intercept) value of -1.483, $\beta_1$ is estimated by the `conc` value of 0.1814, $\beta_2$ is estimated by the `I(conc^2)` value of -0.0003063, and $\beta_3$ is estimated by the `I(conc^3)` value of 1.601e-07, which translates to 0.0000001601.

Because the estimate of the $\beta_3$ term is positive, this cubic model will "open up". In other words, as the function moves from left to right, it will go off to positive infinity (up). If the term would have been negative, then the function would head to negative infinity (down) instead.

$$
\hat{Y}_i = \overbrace{-1.483}^\text{y-int} + \overbrace{0.1814}^{\stackrel{\text{slope}}{\text{term}}} X_{i} + \overbrace{-0.0003063}^{\stackrel{\text{quadratic}}{\text{term}}} X_{i}^2 + \overbrace{1.601e-07}^{\stackrel{\text{cubic}}{\text{term}}} X_{i}^3
$$

The regression function is drawn as follows. Be sure to look at the "Code" to understand how this graph was created using the ideas in the equation above.

<table>
<tr>
<td>

**Using Base R**

```{r}
plot(uptake ~ conc, data=CO2, col="skyblue", pch=21, bg="gray83", main="Cubic Model using CO2 data set", cex.main=1)

#get the "Estimates" automatically:
b <- coef(lm.cubic)
# Then b will have 4 estimates:
# b[1] is the estimate of beta_0: -1.483
# b[2] is the estimate of beta_1: 0.1814
# b[3] is the estimate of beta_2: -0.0003063
# b[4] is the estimate of beta_3: 1.601e-07
curve(b[1] + b[2]*x + b[3]*x^2 + b[4]*x^3, col="skyblue", lwd=2, add=TRUE)
```
  
</td>
<td>

**Using ggplot2**

```{r}
#get the "Estimates" automatically:
b <- coef(lm.cubic)
# Then b will have 4 estimates:
# b[1] is the estimate of beta_0: -1.483
# b[2] is the estimate of beta_1: 0.1814
# b[3] is the estimate of beta_2: -0.0003063
# b[4] is the estimate of beta_3: 1.601e-07

ggplot(CO2, aes(y=uptake, x=conc)) +
  geom_point(pch=21, bg="gray83", color="skyblue") +
  #geom_smooth(method="lm", se=F, formula = y ~ poly(x, 3)) + #easy way, but the more involved manual way using stat_function (see below) is more dynamic.
  stat_function(fun = function(x) b[1] + b[2]*x + b[3]*x^2 + b[4]*x^3, color="skyblue") +
  labs(title="Cubic Model using CO2 data set") 
  

```

</td>
</tr>
</table>

It should be stated, that the cubic function is not the best fit for this data. However, it is a lot better than just a simple line, or a quadratic model, as shown below.

```{r}
plot(uptake ~ conc, data=CO2, col="skyblue", pch=21, bg="gray83", main="Cubic Model using CO2 data set", cex.main=1)

#get the "Estimates" automatically:
b <- coef(lm.cubic)
# Then b will have 4 estimates:
# b[1] is the estimate of beta_0: -1.483
# b[2] is the estimate of beta_1: 0.1814
# b[3] is the estimate of beta_2: -0.0003063
# b[4] is the estimate of beta_3: 1.601e-07
curve(b[1] + b[2]*x + b[3]*x^2 + b[4]*x^3, col="skyblue", lwd=2, add=TRUE)
b <- coef(lm(uptake ~ conc + I(conc^2), data=CO2))
curve(b[1] + b[2]*x + b[3]*x^2, col="firebrick", lwd=2, add=TRUE)
b <- coef(lm(uptake ~ conc, data=CO2))
curve(b[1] + b[2]*x, col="orange", lwd=2, add=TRUE)

```





  </p>
</div>







<div id="LearnMoreTwoLinesModel" class="tabcontent">
  <p>
  
<table>
  <tr>
  <td>

```{r, fig.width=1.5, fig.height=1.5, echo=FALSE}
par(mai=c(0.1,0.1,0.3,0.1))
plot(mpg ~ qsec, data=mtcars, col=c("skyblue","orange")[as.factor(am)], pch=21, bg="gray83", xaxt='n', yaxt='n', main="", axes=FALSE)
lm.2lines <- lm(mpg ~ qsec + am + qsec:am, data=mtcars)
b <- coef(lm.2lines)
curve(b[1] + b[2]*x, col="skyblue", lwd=2, add=TRUE)
curve(b[1] + b[3] + (b[2] + b[4])*x, col="orange", lwd=2, add=TRUE)
```

  </td><td style="text-align: center;padding-left:15px;">  
  
$$
 Y_i = \overbrace{\underbrace{\beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{1i} X_{2i}}_{E\{Y_i\}}}^\text{Two-lines Model} + \epsilon_i 
$$

$$
 X_{2i} = \left\{\begin{array}{ll} 1, & \text{Group B} \\ 0, & \text{Group A} \end{array}\right.
$$

  </td>
  </tr>
</table>


The so called "two-lines" model uses a quantitative $X_{1i}$ variable and a 0,1 indicator variable $X_{2i}$. It is a basic example of how a "dummy variable" or "indicator variable" can be used to turn qualitative variables into quantitative terms. In this case, the indicator variable $X_{2i}$, which is either 0 or 1, produces two separate lines: one line for Group A, and one line for Group B.

| Parameter | Effect |
|-----------|-------------------------------------------------------------------------|
| $\beta_0$ | Y-intercept of the Model. |
| $\beta_1$ | Controls the slope of the "base-line" of the model, the "Group 0" line. |
| $\beta_2$ | Controls the **change in y-intercept** for the second line in the model as compared to the y-intercept of the "base-line" line. |
| $\beta_3$ | Called the "interaction" term. Controls the **change in the slope** for the second line in the model as compared to the slope of the "base-line" line. |




**An Example**

Using the `mtcars` data set, we run the following "two-lines" regression. Note that `am` has only 0 or 1 values: `View(mtcars)`.

$$
  \underbrace{Y_i}_\text{mpg} \underbrace{=}_{\sim} \overbrace{\beta_0}^{\stackrel{\text{y-int}}{\text{baseline}}} + \overbrace{\beta_1}^{\stackrel{\text{slope}}{\text{baseline}}} \underbrace{X_{1i}}_\text{qsec} + \overbrace{\beta_2}^{\stackrel{\text{change in}}{\text{y-int}}}  \underbrace{X_{2i}}_\text{am} + \overbrace{\beta_3}^{\stackrel{\text{change in}}{\text{slope}}} \underbrace{X_{1i}X_{2i}}_\text{qsec:am} + \epsilon_i
$$


<a href="javascript:showhide('twolinesregressionexamplecode')">
<div class="hoverchunk">
<span class="tooltipr">
lm.2lines <- 
  <span class="tooltiprtext">A name we made up for our "two-lines" regression.</span>
</span><span class="tooltipr">
lm(
  <span class="tooltiprtext">R function lm used to perform linear regressions in R. The lm stands for "linear model".</span>
</span><span class="tooltipr">
mpg
  <span class="tooltiprtext">Y-variable, should be quantitative.</span>
</span><span class="tooltipr">
&nbsp;~&nbsp;
<span class="tooltiprtext">The tilde `~` is what lm(...) uses to state the regression equation $Y_i = ...$. Notice that the `~` is not followed by $\beta_0 + \beta_1$ like $Y_i = ...$. Instead, $X_{1i}$ is the first term following `~`. This is because $\beta$'s are going to be estimated by the lm(...). These estimates can be found using summary(lmObject).</span>
</span><span class="tooltipr">
qsec
  <span class="tooltiprtext">$X_{1i}$, should be quantitative.</span>
</span><span class="tooltipr">
&nbsp;+&nbsp;
  <span class="tooltiprtext">The plus `+` is used between each term in the model. Note that only the x-variables are included in the lm(...) from the $Y_i = ...$ model. No beta's are included.</span>
</span><span class="tooltipr">
am
  <span class="tooltiprtext">$X_{2i}$, an indicator or 0,1 variable. This term allows the y-intercept of the two lines to differ.</span>
</span><span class="tooltipr">
&nbsp;+&nbsp;
  <span class="tooltiprtext">The plus `+` is used between each term in the model. Note that only the x-variables are included in the lm(...) from the $Y_i = ...$ model. No beta's are included.</span>
</span><span class="tooltipr">
qsec:am
  <span class="tooltiprtext">$X_{1i}X_{2i}$ the interaction term. This allows the slopes of the two lines to differ.</span>
</span><span class="tooltipr">
, data=mtcars
  <span class="tooltiprtext">This is the data set we are using for the regression.</span>
</span><span class="tooltipr">
)  
  <span class="tooltiprtext">Closing parenthsis for the lm(...) function.</span>
</span><span class="tooltipr">
&nbsp;&nbsp;&nbsp;&nbsp;  
  <span class="tooltiprtext">Press Enter to run the code.</span>
</span><span class="tooltipr" style="float:right;">
&nbsp;...&nbsp; 
  <span class="tooltiprtext">Click to View Output.</span>
</span>
</div>
</a>
<div id="twolinesregressionexamplecode" style="display:none;">
Pay special attention to how the lm(...) code uses $Y_i \sim X_{1i} + X_{2i} + X_{1i}X_{2i}$ and drops all $\beta$'s and $\epsilon$ from the model statement. This is because the estimates for the $\beta$'s and $\epsilon$ are given by the output of the lm(...) funtion in the "Estimates" column of summary(....) and in `lm.2lines$residuals`.
</div>


```{r}
lm.2lines <- lm(mpg ~ qsec + am + qsec:am, data=mtcars)
pander(summary(lm.2lines)$coefficients)
```

The **estimates** shown above approximate the $\beta$'s in the regression model: $\beta_0$ is estimated by the (Intercept), $\beta_1$ is estimated by the `qsec` value of 1.439, $\beta_2$ is estimated by the `am` value of -14.51, and $\beta_3$ is estimated by the `qsec:am` value of 1.321.

This gives two separate equations of lines.

**Automatic Transmission (am==0, $X_{2i} = 0$) Line**

$$
\hat{Y}_i = \overbrace{-9.01}^{\stackrel{\text{y-int}}{\text{baseline}}} + \overbrace{1.439}^{\stackrel{\text{slope}}{\text{baseline}}} X_{1i}
$$

**Manual Transmission (am==1 , $X_{2i} = 1$) Line**

$$
\hat{Y}_i = \underbrace{(\overbrace{-9.01}^{\stackrel{\text{y-int}}{\text{baseline}}} + \overbrace{-14.51}^{\stackrel{\text{change in}}{\text{y-int}}})}_{\stackrel{\text{y-intercept}}{-23.52}} + \underbrace{(\overbrace{1.439}^{\stackrel{\text{slope}}{\text{baseline}}} +\overbrace{1.321}^{\stackrel{\text{change in}}{\text{slope}}})}_{\stackrel{\text{slope}}{2.76}} X_{1i}
$$

These lines are drawn as follows. Be sure to look at the "Code" to understand how this graph was created using the ideas in the two equations above.

<table>
<tr>
<td>

**Using Base R**

```{r}
plot(mpg ~ qsec, data=mtcars, col=c("skyblue","orange")[as.factor(am)], pch=21, bg="gray83", main="Two-lines Model using mtcars data set", cex.main=1)

legend("topleft", legend=c("Baseline (am==0)", "Changed-line (am==1)"), bty="n", lty=1, col=c("skyblue","orange"), cex=0.8)

#get the "Estimates" automatically:
b <- coef(lm.2lines)
# Then b will have 4 estimates:
# b[1] is the estimate of beta_0: -9.0099
# b[2] is the estimate of beta_1:  1.4385
# b[3] is the estimate of beta_2: -14.5107
# b[4] is the estimate of beta_3: 1.3214
curve(b[1] + b[2]*x, col="skyblue", lwd=2, add=TRUE)  #baseline (in blue)
curve((b[1] + b[3]) + (b[2] + b[4])*x, col="orange", lwd=2, add=TRUE) #changed line (in orange)
```
  
</td>
<td>

**Using ggplot2**

```{r}
#get the "Estimates" automatically:
b <- coef(lm.2lines)
# Then b will have 4 estimates:
# b[1] is the estimate of beta_0: -9.0099
# b[2] is the estimate of beta_1:  1.4385
# b[3] is the estimate of beta_2: -14.5107
# b[4] is the estimate of beta_3: 1.3214

ggplot(mtcars, aes(y=mpg, x=qsec, color=factor(am))) +
  geom_point(pch=21, bg="gray83") +
  #geom_smooth(method="lm", se=F) + #easy way, but only draws the full interaction model. The manual way using stat_function (see below) is more involved, but more dynamic.
  stat_function(fun = function(x) b[1] + b[2]*x, color="skyblue") + #am==0 line
  stat_function(fun = function(x) (b[1]+b[3]) + (b[2]+b[4])*x,color="orange") + #am==1 line 
  scale_color_manual(name="Transmission (am)", values=c("skyblue","orange")) +
  labs(title="Two-lines Model using mtcars data set") 
  

```

</td>
</tr>
</table>


</p>
</div>






<div id="LearnMorethreeDModel" class="tabcontent">
  <p>


<table>
  <tr>
  <td>

```{r volcano, fig.width=1.5, fig.height=1.5, echo=FALSE}
par(mai=c(0.1,0.1,0.3,0.1))
  z <- volcano
  x = seq(0, 1, length.out = nrow(z))
  y = seq(0, 1, length.out = ncol(z))
  xlim = range(x, finite = TRUE) 
  ylim = range(y, finite = TRUE) 
  zlim = range(z, finite = TRUE)
  nlevels = 20
  levels <- pretty(zlim, nlevels)
  plot.new()
  mycolorpalette <- colorRampPalette(c("skyblue2", "orange"))
  .filled.contour(x, y, z, levels, col=mycolorpalette(length(levels)-1))
```

  </td><td style="text-align: center;padding-left:15px;">
  
$$
 Y_i = \overbrace{\underbrace{\beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{1i}X_{2i}}_{E\{Y_i\}}}^\text{3D Model} + \epsilon_i
$$

  </td>
  </tr>
</table>



The so called "3D" regression model uses two different quantitative x-variables, an $X_{1i}$ and an $X_{2i}$. Unlike the two-lines model where $X_{2i}$ could only be a 0 or a 1, this $X_{2i}$ variable is quantitative, and can take on any quantitative value.

| Parameter | Effect |
|-----------|--------------------------------------------------------------------------|
| $\beta_0$ | Y-intercept of the Model |
| $\beta_1$ | Slope of the line in the $X_1$ direction. |
| $\beta_2$ | Slope of the line in the $X_2$ direction. |
| $\beta_3$ | Interaction term that allows the model, which is a plane in three-dimensional space, to "bend". If this term is zero, then the regression surface is just a flat plane. |




**An Example**

Here is what a 3D regression looks like when there is no interaction term. The two x-variables of `Month` and `Temp` are being used to predict the y-variable of `Ozone`.

$$
  \underbrace{Y_i}_\text{Ozone} \underbrace{=}_{\sim} \overbrace{\beta_0}^{\stackrel{\text{y-int}}{\text{baseline}}} + \overbrace{\beta_1}^{\stackrel{\text{slope}}{\text{baseline}}} \underbrace{X_{1i}}_\text{Temp} + \overbrace{\beta_2}^{\stackrel{\text{change in}}{\text{y-int}}}  \underbrace{X_{2i}}_\text{Month} + \epsilon_i
$$

```{r}
air_lm <- lm(Ozone ~ Temp + Month, data= airquality)
pander(air_lm$coefficients)
```

Notice how the slope, $\beta_1$, in the "Temp" direction is estimated to be 2.659 and the slope in the "Month" direction, $\beta_2$, is estimated to be -3.522. Also, the y-intercept, $\beta_0$, is estimated to be -139.6.


```{r, message=FALSE, warning=FALSE}
## Hint: library(car) has a scatterplot 3d function which is simple to use
#  but the code should only be run in your console, not knit.

## library(car)
## scatter3d(Y ~ X1 + X2, data=yourdata)



## To embed the 3d-scatterplot inside of your html document is harder.
#library(plotly)
#library(reshape2)

#Perform the multiple regression
air_lm <- lm(Ozone ~ Temp + Month, data= airquality)

#Graph Resolution (more important for more complex shapes)
graph_reso <- 0.5

#Setup Axis
axis_x <- seq(min(airquality$Temp), max(airquality$Temp), by = graph_reso)
axis_y <- seq(min(airquality$Month), max(airquality$Month), by = graph_reso)

#Sample points
air_surface <- expand.grid(Temp = axis_x, Month = axis_y, KEEP.OUT.ATTRS=F)
air_surface$Z <- predict.lm(air_lm, newdata = air_surface)
air_surface <- acast(air_surface, Month ~ Temp, value.var = "Z") #y ~ x

#Create scatterplot
plot_ly(airquality, 
        x = ~Temp, 
        y = ~Month, 
        z = ~Ozone,
        text = rownames(airquality), 
        type = "scatter3d", 
        mode = "markers") %>%
  add_trace(z = air_surface,
            x = axis_x,
            y = axis_y,
            type = "surface")

```


Here is a second view of this same regression with what is called a contour plot, contour map, or density plot.

```{r, eval=FALSE}
mycolorpalette <- colorRampPalette(c("skyblue2", "orange"))
filled.contour(x=axis_x, y=axis_y, z=matrix(air_surface$Z, length(axis_x), length(axis_y)), col=mycolorpalette(26))
```


**Including the Interaction Term**

Here is what a 3D regression looks like when the interaction term is present. The two x-variables of `Month` and `Temp` are being used to predict the y-variable of `Ozone`.

$$
  \underbrace{Y_i}_\text{Ozone} \underbrace{=}_{\sim} \overbrace{\beta_0}^{\stackrel{\text{y-int}}{\text{baseline}}} + \overbrace{\beta_1}^{\stackrel{\text{slope}}{\text{baseline}}} \underbrace{X_{1i}}_\text{Temp} + \overbrace{\beta_2}^{\stackrel{\text{change in}}{\text{y-int}}}  \underbrace{X_{2i}}_\text{Month} + \overbrace{\beta_3}^{\stackrel{\text{change in}}{\text{slope}}} \underbrace{X_{1i}X_{2i}}_\text{Temp:Month} + \epsilon_i
$$

```{r}
air_lm <- lm(Ozone ~ Temp + Month + Temp:Month, data= airquality)
pander(air_lm$coefficients)
```

Notice how all coefficient estimates have changed. The y-intercept, $\beta_0$ is now estimated to be $-3.915$. The slope term, $\beta_1$, in the Temp-direction is estimated as $0.77$, while the slope term, $\beta_2$, in the Month-direction is estimated to be $-23.01$. This change in estimated coefficiets is due to the presence of the interaction term's coefficient, $\beta_3$, which is estimated to be $0.2678$. As you should notice in the graphic, the interaction model allows the "slopes" in each direction to change, creating a "curved" surface for the regression surface instead of a flat surface.


```{r, message=FALSE, warning=FALSE}
#Perform the multiple regression
air_lm <- lm(Ozone ~ Temp + Month + Temp:Month, data= airquality)

#Graph Resolution (more important for more complex shapes)
graph_reso <- 0.5

#Setup Axis
axis_x <- seq(min(airquality$Temp), max(airquality$Temp), by = graph_reso)
axis_y <- seq(min(airquality$Month), max(airquality$Month), by = graph_reso)

#Sample points
air_surface <- expand.grid(Temp = axis_x, Month = axis_y, KEEP.OUT.ATTRS=F)
air_surface <- air_surface %>% mutate(Z=predict.lm(air_lm, newdata = air_surface))
air_surface <- acast(air_surface, Month ~ Temp, value.var = "Z") #y ~ x

#Create scatterplot
plot_ly(airquality, 
        x = ~Temp, 
        y = ~Month, 
        z = ~Ozone,
        text = rownames(airquality), 
        type = "scatter3d", 
        mode = "markers") %>%
  add_trace(z = air_surface,
            x = axis_x,
            y = axis_y,
            type = "surface")
```

And here is that same plot as a contour plot.

```{r, eval=FALSE}
air_surface <- expand.grid(Temp = axis_x, Month = axis_y, KEEP.OUT.ATTRS=F)
air_surface$Z <- predict.lm(air_lm, newdata = air_surface)
mycolorpalette <- colorRampPalette(c("skyblue2", "orange"))
filled.contour(x=axis_x, y=axis_y, z=matrix(air_surface$Z, length(axis_x), length(axis_y)), col=mycolorpalette(27))
```

</p>
</div>





<div id="LearnMoreHDModel" class="tabcontent" style="display:none;">
  <p>


<table>
  <tr>
  <td>

```{r, fig.width=1.5, fig.height=1.5, echo=FALSE}
par(mai=c(0.1,0.1,0.3,0.1))
mycolorpalette <- colorRampPalette(c("skyblue2", "orange"))
palette(rev(mycolorpalette(10)))
plot(Ozone ~ Solar.R, data=airquality, pch=21, col="gray83", bg=as.factor(pretty(Temp, 10)), yaxt='n', xaxt='n', axes=FALSE, cex=Month/18 *5)
lm.hd <- lm(Ozone ~ Solar.R + Month + Temp + Temp:Month:Solar.R, data=airquality)
b <- coef(lm.hd)
curve(b[1] + b[2]*x + b[3]*1 + b[4]*55 + b[5]*55*1*x, add=TRUE, col=palette()[1], lwd=1)
curve(b[1] + b[2]*x + b[3]*5 + b[4]*60 + b[5]*60*5*x, add=TRUE, col=palette()[2], lwd=2)
curve(b[1] + b[2]*x + b[3]*7 + b[4]*65 + b[5]*65*7*x, add=TRUE, col=palette()[3], lwd=3)
curve(b[1] + b[2]*x + b[3]*7 + b[4]*75 + b[5]*75*7*x, add=TRUE, col=palette()[5], lwd=3)
curve(b[1] + b[2]*x + b[3]*12 + b[4]*80 + b[5]*80*12*x, add=TRUE, col=palette()[8], lwd=4)
curve(b[1] + b[2]*x + b[3]*5 + b[4]*85 + b[5]*85*5*x, add=TRUE, col=palette()[9], lwd=2)
curve(b[1] + b[2]*x + b[3]*2 + b[4]*95 + b[5]*100*2*x, add=TRUE, col=palette()[10], lwd=3)
```

  </td><td style="text-align: center;padding-left:15px;">
  
$$
  Y_i = \overbrace{\underbrace{\beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \ldots + \beta_{p-1}X_{p-1,i}}_{E\{Y_i\}}}^\text{"High Dimensional Models"} + \epsilon_i
$$

  </td>
  </tr>
</table>



The so called "HD", or "High Dimensional", regression model uses three or more different quantitative x-variables, an $X_{1i}$, an $X_{2i}$, and at least an $X_{3i}$, but could use many, many other variables as well. Unlike the 3D model where the final regression could be shown as either a contour plot or a 3D-graphic, the high dimensional model exists in 4 or more dimensions. Thus, it is impossible to graph this model in its full form. Further, it isn't really even possible to "mentally connect" with this type of model is it exists beyond what our 3D minds can really comprehend.

| Parameter | Effect |
|-----------|--------------------------------------------------------------------------|
| $\beta_0$ | Y-intercept of the Model |
| $\beta_1$ | Slope of the line in the $X_1$ direction. |
| $\beta_2$ | Slope of the line in the $X_2$ direction. |
| $...$     | Slopes in other directions depending on how many other variables are included in the model. |
| $\beta_{p-1}$ | Final term in the model where there are $p$ total $\beta$'s. The reason for the $p-1$ on the last term is because we started with $\beta_0$ for the first term, leaving $\beta_{p-1}$ as the last term. |



**An Example**

Suppose we used three x-variables of `Wind`, `Temp`, and `Solar.R` to predict the y-variable of `Ozone`.

$$
  \underbrace{Y_i}_\text{Ozone} \underbrace{=}_{\sim} \overbrace{\beta_0}^{\stackrel{\text{y-int}}{\text{baseline}}} + \overbrace{\beta_1}^{\stackrel{\text{slope in}}{\text{Wind Direction}}} \underbrace{X_{1i}}_\text{Wind} + \overbrace{\beta_2}^{\stackrel{\text{slope in}}{\text{Temp Direction}}}  \underbrace{X_{2i}}_\text{Temp} + \overbrace{\beta_3}^{\stackrel{\text{slope in}}{\text{Solar.R Direction}}}  \underbrace{X_{3i}}_\text{Solar.R} + \epsilon_i
$$

```{r}
air_lm <- lm(Ozone ~ Wind + Temp + Solar.R, data= airquality)
pander(air_lm$coefficients)
```

Notice how the slope, $\beta_1$, in the "Wind" direction is estimated to be -3.334. The slope in the "Temp" direction, $\beta_2$, is estimated to be 1.652. The slope in the "Solar.R" direction, $\beta_3$, is estimated to be 0.05982. Also, the y-intercept, $\beta_0$, is estimated to be -64.34.

Visualizing this model is not really possible in its full form. However, we can draw the regression from three different angles or vantage points. This is a limited view of the full regression model, but at least provides some visual understanding. To do this, we draw $Y$ against each $X$-variable in separate scatterplots, one for each $X$-variable used in our model.

```{r message=FALSE, warning=FALSE}
b <- coef(air_lm)

par(mfrow=c(1,3))

  plot(Ozone ~ Wind, data=airquality)
  curve(b[1] + b[2]*x + b[3]*79 + b[4]*205, add=TRUE, col="skyblue")
  # The x-variable of this plot is "Wind"
  # The values of Temp=79 and Solar.R=205 are fixed at some interesting value,
  # in this case, their respective medians.

  plot(Ozone ~ Temp, data=airquality)
  curve(b[1] + b[2]*9.7 + b[3]*x + b[4]*205, add=TRUE, col="orange")
  # The x-variable of this plot is "Temp"
  # The values of Wind=9.7 and Solar.R=205 are fixed at some interesting value,
  # in this case, their respective medians.
  
  plot(Ozone ~ Solar.R, data=airquality)
  curve(b[1] + b[2]*9.7 + b[3]*79 + b[4]*x, add=TRUE, col="firebrick")
  # The x-variable of this plot is "Solar.R"
  # The values of Wind = 9.7 and Temp=79 are fixed at some interesting value,
  # in this case, their respective medians.

```



</p>
</div>


----



The coefficient $\beta_j$ is interpreted as the change in the expected value of $Y$ for a unit increase in $X_{j}$, holding all other variables constant, for $j=1,\ldots,p-1$. However, this interpretation breaks down when higher order terms (like $X^2$) or interaction terms (like $X1:X2$) are included in the model.

See the **Explanation** tab for details about possible hypotheses here.

----

</div>

















<div style="padding-left:125px;">

The probability that $Y_i = 1$ given the observed data $(x_{i1},\ldots,x_{ip})$ is called $\pi_i$ and is modeled by the equation

$$
  P(Y_i = 1|\, x_{i1},\ldots,x_{ip}) = \frac{e^{\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}}}{1+e^{\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip} }} = \pi_i
$$

The coefficents $\beta_0,\beta_1,\ldots,\beta_p$ are difficult to interpret directly. Typically $e^{\beta_k}$ for $k=0,1,\ldots,p$ is interpreted instead. The value of $e^{\beta_k}$ denotes the relative change in the odds that $Y_i=1$. The odds that $Y_i=1$ are $\frac{\pi_i}{1-\pi_i}$.

----


**Examples:** [GSS](./Analyses/Logistic Regression/Examples/GSSMultipleLogisticReg.html) 

----

</div>

### R Instructions

<div style="padding-left:125px;">


**Console** Help Command: `?glm()`

#### Perform the Logistic Regression

To perform a logistic regression in R use the commands

<a href="javascript:showhide('logistic2')">
<div class="hoverchunk">
<span class="tooltipr">
YourGlmName
  <span class="tooltiprtext">This is some name you come up with that will become the R object that stores the results of your logistic regression `glm()` command.</span>
</span><span class="tooltipr">
&nbsp;<-&nbsp;
  <span class="tooltiprtext">This is the "left arrow" assignment operator that stores the results of your `glm()` code into `YourGlmName`.</span>
</span><span class="tooltipr">
glm(
  <span class="tooltiprtext">glm( is an R function that stands for "General Linear Model". It works in a similar way that the `lm(` function works except that it requires a `family=` option to be specified at the end of the command.</span>
</span><span class="tooltipr">
Y&nbsp;
  <span class="tooltiprtext">Y is your binary response variable. It must consist of only 0's and 1's. Since TRUE's = 1's and FALSE's = 0's in R, Y could be a logical statement like (Price > 100) or (Animal == "Cat") if your Y-variable wasn't currently coded as 0's and 1's.</span>
</span><span class="tooltipr">
~&nbsp;
  <span class="tooltiprtext">The tilde symbol ~ is used to tell R that Y should be treated as a function of the explanatory variable X.</span>
</span><span class="tooltipr">
X1
  <span class="tooltiprtext">X1 is the first explanatory variable (typically quantitative) that will be used to explain the probability that the response variable Y is a 1.</span>
</span><span class="tooltipr">
\*
  <span class="tooltiprtext">The times symbol allows a shortcut for writing X1 + X2 + X1:X2 = X1*X2.</span>
</span><span class="tooltipr">
X2&nbsp;
  <span class="tooltiprtext">X2 is second the explanatory variable either quantitative or qualitative that will be used to explain the probability that the response variable Y is a 1.</span>
</span><span class="tooltipr">
..., 
  <span class="tooltiprtext">In theory, you could have many other explanatory variables, interaction terms, or even squared, cubed, or other transformations of terms added to this model.</span>
</span><span class="tooltipr">
&nbsp;data = NameOfYourDataset,  
  <span class="tooltiprtext">NameOfYourDataset is the name of the dataset that contains Y and X. In other words, one column of your dataset would be called Y and another column would be called X.</span>
</span><span class="tooltipr">
&nbsp;family=binomial)
  <span class="tooltiprtext">The family=binomial command tells the `glm(` function to perform a logistic regression. It turns out that `glm` can perform many different types of regressions, but we only study it as a tool to perform a logistic regression in this course.</span>
</span><br/><span class="tooltipr">
summary(YourGlmName)
  <span class="tooltiprtext">The `summary` command allows you to print the results of your logistic regression that were previously saved in `YourGlmName`.</span>
</span>
</div>
</a>

<div id="logistic2" style="display:none;">

Example output from a regression. Hover each piece to learn more.

<table class="rconsole">
<tr>
<td>
<span class="tooltiprout">
  Call:<br/>
 glm(formula = weight > 100 ~ Time*Diet, family = binomial, 
    data = ChickWeight)
    <span class="tooltiprouttext">This is simply a statement of your original glm(...) "call" that you made when performing your regression. It allows you to verify that you ran what you thought you ran in the glm(...).</span>
</span>
</td>
</tr>
</table>
      

<br/>

<table class="rconsole">
<tr>
<td colspan="2">
  <span class="tooltiprout">
    Deviance Residuals:
      <span class="tooltiprouttext">Deviance residuals are a measure of how far the fitted probability for $\pi_i$ has differed from the actual outcome of $Y_i$ in terms of the log of the fitted probability space. (This is a fairly complicated idea.)
  </span>
</td>
</tr>
<tr>
<td align="right">
  <span class="tooltiprout">
    Min<br/>
    &nbsp; -2.9446
      <span class="tooltiprouttext">"min" gives the value of the residual that is furthest below the regression line. Ideally, the magnitude of this value would be about equal to the magnitude of the largest positive residual (the max) because the hope is that the residuals are normally distributed around the line.</span>
  </span>
</td><td align="right">   
  <span class="tooltiprout">
    1Q<br/>
    &nbsp; -0.2364
      <span class="tooltiprouttext">"1Q" gives the first quartile of the residuals, which will always be negative, and ideally would be about equal in magnitude to the third quartile.</span>
  </span>
</td><td align="right"> 
  <span class="tooltiprout">
    Median<br/>
    &nbsp; 0.0000
      <span class="tooltiprouttext">"Median" gives the median of the residuals, which would ideally would be about equal to zero. Note that because the regression line is the least squares line, the mean of the residuals will ALWAYS be zero, so it is never included in the output summary.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    3Q<br/>
    &nbsp; 0.2776
      <span class="tooltiprouttext">"3Q" gives the third quartile of the residuals, which would ideally would be about equal in magnitude to the first quartile. In this case, it is pretty close, which helps us see that the first quartile of residuals on either side of the line is behaving fairly normally.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    Max</br>
    &nbsp; 1.9968
      <span class="tooltiprouttext">"Max" gives the maximum positive residuals, which would ideally would be about equal in magnitude to the minimum residual. In this case, it is much larger than the minimum, which helps us see that the residuals are likely right skewed.</span>
  </span>
</td>
</tr>
</table>

<br/>


<table class="rconsole">
<tr>
<td colspan="2">
  <span class="tooltiprout">
    Coefficients:
      <span class="tooltiprouttext">Notice that in your glm(...) you used only $Y$ and $X$. You did type out any coefficients, i.e., the $\beta_0$ or $\beta_1$ of the regression model. These coefficients are estimated by the glm(...) function and displayed in this part of the output along with standard errors, t-values, and p-values.</span>
  </span>
</td>
</tr><tr>
<td align="left"></td>
<td align="right">
  <span class="tooltiprout">
    &nbsp; Estimate
      <span class="tooltiprouttext">To learn more about the "Estimates" of the "Coefficients" see the "Explanation" tab, "Estimating the Model Parameters" section for details.</span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; Std. Error
      <span class="tooltiprouttext">To learn more about the "Standard Errors" of the "Coefficients" see the "Explanation" tab, "Inference for the Model Parameters" section.</span>
  </span>    
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; z value 
      <span class="tooltiprouttext">The Z-score testing the hypothesis that $\beta_j=0$</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; Pr(>|z|)
      <span class="tooltiprouttext">The p-value for the corresponding Z-score.</span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    (Intercept)
      <span class="tooltiprouttext">This always says "Intercept" for any glm(...) you run in R. That is because R always assumes there is a y-intercept for your regression function.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; -4.97603
      <span class="tooltiprouttext">This is the estimate of the y-intercept, $\beta_0$. It is called $b_0$. It is the average y-value when X is zero.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.65316
      <span class="tooltiprouttext">This is the standard error of $b_0$. It tells you how much $b_0$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    -7.618
      <span class="tooltiprouttext">Z-score for $b_0$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    2.57e-14
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_0 = 0$. It measures the probability of observing a Z-value as extreme as the one observed. To compute it yourself in R, use `pnorm(-abs(your z-value))*2`.</span>
  </span>
</td><td align="left">
<span class="tooltiprout">
    \*\*\*
      <span class="tooltiprouttext">This is called a "star". Three stars means significant at the 0.01 level of $\alpha$.</span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Time
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.39111 
      <span class="tooltiprouttext">This is the estimate of the slope, $\beta_1$. It is called $b_1$. It is the change in the average y-value as X is increased by 1 unit.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.04979
      <span class="tooltiprouttext">This is the standard error of $b_1$. It tells you how much $b_1$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    7.854
      <span class="tooltiprouttext">Z-score for $b_1$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    4.02e-15
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_1 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td><td align="left">
<span class="tooltiprout">
    \*\*\*
      <span class="tooltiprouttext">This is called a "star". Three stars means significant at the 0.01 level of $\alpha$.</span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Diet2
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.58283
      <span class="tooltiprouttext">This is the estimate of the change in the intercept, $\beta_2$. It is called $b_2$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 1.04061
      <span class="tooltiprouttext">This is the standard error of $b_2$. It tells you how much $b_2$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.560
      <span class="tooltiprouttext">Z-score for $b_2$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.5754
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_2 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Diet3
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; -8.44476
      <span class="tooltiprouttext">This is the estimate of the change in the intercept, $\beta_3$. It is called $b_3$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 4.32324
      <span class="tooltiprouttext">This is the standard error of $b_3$. It tells you how much $b_3$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    -1.953
      <span class="tooltiprouttext">Z-score for $b_3$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.0508
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_3 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td><td align="left">
<span class="tooltiprout">
    .
      <span class="tooltiprouttext">This is called a period. It means significant at the 0.1 level of $\alpha$.</span>
  </span>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Diet4
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; -67.41995
      <span class="tooltiprouttext">This is the estimate of the change in the intercept, $\beta_4$. It is called $b_4$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 3768.10023
      <span class="tooltiprouttext">This is the standard error of $b_4$. It tells you how much $b_4$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    -0.018
      <span class="tooltiprouttext">Z-score for $b_4$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.9857
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_2 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Time:Diet2
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.02391
      <span class="tooltiprouttext">This is the estimate of the change in the slope, $\beta_5$. It is called $b_5$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.08679
      <span class="tooltiprouttext">This is the standard error of $b_5$. It tells you how much $b_5$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.275
      <span class="tooltiprouttext">Z-score for $b_5$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.7830
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_2 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Time:Diet3
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 1.20955
      <span class="tooltiprouttext">This is the estimate of the change in the slope, $\beta_6$. It is called $b_6$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 0.51249
      <span class="tooltiprouttext">This is the standard error of $b_6$. It tells you how much $b_6$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    2.360
      <span class="tooltiprouttext">Z-score for $b_6$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.0183
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_2 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td><td align="left">
<span class="tooltiprout">
    \*
      <span class="tooltiprouttext">This is called a "star". One star means significant at the 0.05 level of $\alpha$.</span>
  </span>
</td>
</td>
</tr><tr>
<td align="left">
  <span class="tooltiprout">
    Time:Diet4
      <span class="tooltiprouttext">This is always the name of your X-variable in your glm(Y ~ X, ...).</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 8.83168
      <span class="tooltiprouttext">This is the estimate of the change in the slope, $\beta_7$. It is called $b_7$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp; 471.01259
      <span class="tooltiprouttext">This is the standard error of $b_7$. It tells you how much $b_7$ varies from sample to sample. The closer to zero, the better.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.019
      <span class="tooltiprouttext">Z-score for $b_7$.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    0.9850
      <span class="tooltiprouttext">This is the p-value of the test of the hypothesis that $\beta_2 = 0$. To compute it yourself in R, use `pt(-abs(your t-value), df of your regression)*2`</span>
  </span>
</td>
</td>
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span>
  \-\-\-
  </span>
</td>
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
    Signif. codes:  0 ‘\*\*\*’ 0.001 ‘\*\*’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      <span class="tooltiprouttext">These "codes" explain what significance level the p-value is smaller than based on how many "stars" * the p-value is labeled with in the Coefficients table above.</span>
  </span>
</td>
</tr>
</table>

<br/>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
    (Dispersion parameter for binomial family taken to be 1)
      <span class="tooltiprouttext">A simplifying assumption of the logistic regression. This can be changed if you know what you are doing. See ?summary.glm and look at the dispersion option for details. Unless you are pursuing an advanced degree in statistics, it is not recommended that you explore this option.</span>
  </span>
</td>
</tr>
</table>

<br/>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
    &nbsp; Null Deviance:
      <span class="tooltiprouttext">The deviance of the null model.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;800.44
      <span class="tooltiprouttext">In this case, the null deviance is 800.44</span>
  </span>    
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;on 577 degrees of freedom 
      <span class="tooltiprouttext">The residual degrees of freedom for the null model.</span>
  </span> 
</td> 
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
   Residual deviance:
      <span class="tooltiprouttext">The test statistic from the chi-squared goodness of fit test of the logistic regression.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;256.30
      <span class="tooltiprouttext">This is computed by the sum(log(myglm$res^2)). As long as this value is similar in size to the degrees of freedom, then the model is a good fit.</span>
  </span>    
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;on 570 degrees of freedom 
      <span class="tooltiprouttext">This is $n-p$ where $n$ is the sample size and $p$ is the number of parameters in the regression model. In this case, there is a sample size of 578 and five parameters, $\beta_0$, $\beta_1$, $\beta_2$, $\beta_3$, $\beta_4$, $\beta_5$,$\beta_6$, and $\beta_7$ so 578-8 = 570.</span>
  </span> 
</td> 
</tr>
</table>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
   AIC:
      <span class="tooltiprouttext">A form of the Akiake Information Criterion. Smaller values of the AIC indicate a better fitting model.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;272.3
      <span class="tooltiprouttext">If a model for the same Y-variable can be found with an AIC value lower than 272.3, then it is a better fit of the data.</span>
  </span>    
</td> 
</tr>
</table>

<br/>

<table class="rconsole">
<tr>
<td>
  <span class="tooltiprout">
   Number of Fisher Scoring iterations:
      <span class="tooltiprouttext">If you have taken a class in Numerical Analysis, this tells you how many iterations of the maximization algorithm were required before converging to the "Estimates" of the parameters $\beta_0$ and $\beta_1$ found in the summary.</span>
  </span>
</td><td align="right">
  <span class="tooltiprout">
    &nbsp;20
      <span class="tooltiprouttext">This implementation of glm required 20 Fisher Scoring iterations to converge. Fewer iterations hints that the model is a better fit than when many iterations are required.</span>
  </span>    
</td> 
</tr>
</table>


</div>
<br/>




#### Diagnose the Goodness-of-Fit

There are two ways to check the **goodness of fit** of a logistic regression model.

<div style="padding-left:25px;">
**Option 1**: Hosmer-Lemeshow Goodness-of-Fit Test

To check the **goodness of fit** of a logistic regression model where there are **few or no replicated $x$-values** use the Hosmer-Lemeshow Test.

<a href="javascript:showhide('goodnessoffit2')">
<div class="hoverchunk">
<span class="tooltipr">
library(ResourceSelection)
  <span class="tooltiprtext">This loads the ResourceSelection R package so that you can access the hoslem.test() function. You may need to run the code: install.packages("ResourceSelection") first.</span>
</span><br/><span class="tooltipr">
hoslem.test(
  <span class="tooltiprtext">This R function performs the Hosmer-Lemeshow Goodness of Fit Test. See the "Explanation" file to learn about this test.</span>
</span><span class="tooltipr">
YourGlmName
  <span class="tooltiprtext">`YourGlmName` is the name of your glm(...) code that you created previously.</span>
</span><span class="tooltipr">
\$y,&nbsp;
  <span class="tooltiprtext">ALWAYS type a "y" here. This gives you the actual binary (0,1) y-values of your logistic regression. The goodness of fit test will compare these actual values to your predicted probabilities for each value in order to see if the model is a "good fit."</span>
</span><span class="tooltipr">
YourGlmName
  <span class="tooltiprtext">`YourGlmName` is the name you used to save the results of your glm(...) code.</span>
</span><span class="tooltipr">
\$fitted,&nbsp;
  <span class="tooltiprtext">ALWAYS type "fitted" here. This gives you the fitted probabilities $\pi_i$ of your logistic regression.</span>
</span><span class="tooltipr">
g=10)
  <span class="tooltiprtext">The "g=10" is the default option for the value of g. The g is the number of groups to run the goodness of fit test on. Just leave it at 10 unless you are told to do otherwise. Ask your teacher for more information if you are interested.</span>
</span>
</div>
</a>

<div id="goodnessoffit2" style="display:none;">
```{r, echo=FALSE}
myglm <- glm(weight>100 ~ Time*Diet, data = ChickWeight, family = binomial)
library(ResourceSelection)
hoslem.test(myglm$y, myglm$fitted)
```

Note that the null hypothesis of the goodness-of-fit test is that "the logistic regression is a good fit." So we actually don't want to "reject the null" in this case. So a large p-value here means our logistic regression fits the data satisfactorily. A small p-value implies a poor fit and the results of the logistic regression should not be fully trusted.

</div>

<br/>

**Option 2**: Deviance Goodness-of-fit Test

In some cases, there are **many replicated $x$-values** for **all** x-values. Though this is rare, it is good to use the *deviance goodness-of-fit test* whenever this happens.

<a href="javascript:showhide('goodnessoffit1')">
<div class="hoverchunk">
<span class="tooltipr">
pchisq(
  <span class="tooltiprtext">The `pchisq` command allows you to compute p-values from the chi-squared distribution.</span>
</span><span class="tooltipr">
residual deviance,&nbsp;
  <span class="tooltiprtext">The residual deviance is shown at the bottom of the output of your `summary(YourGlmName)` and should be typed in here as a number like 25.3.</span>
</span><span class="tooltipr">
df for residual deviance,&nbsp;
  <span class="tooltiprtext">The df for the residual deviance is also shown at the bottom of the output of your `summary(YourGlmName)`.</span>
</span><span class="tooltipr">
lower.tail=FALSE)
<span class="tooltiprtext">This command ensures you find the probability of the chi-squared distribution being as extreme or more extreme than the observed value of residual deviance.</span>
</span>
</div>
</a>


<div id="multipleGOF2" style="display:none;">
```{r, warning=FALSE}
myglm <- glm(weight>100 ~ Time*Diet, data = ChickWeight, family = binomial)
summary(myglm)

pchisq(256.30, 570, lower.tail = TRUE)
```
</div>


The null hypothesis of the goodness-of-fit test is that the logistic regression is a good fit of the data. So a large p-value (like 0.479) is good because it allows us to trust the results of our logistic regression. When the p-value becomes very small, we must "reject the null" and conclude a poor fit, which implies that we should not trust the results of the logistic regression.

</div>


<br/>

#### Predict Probabilities

To predict the probability that $Y_i=1$ for a given $x$-value, use the code

<a href="javascript:showhide('MultipleLogisticPredict')">
<div class="hoverchunk">
<span class="tooltipr">
predict(
  <span class="tooltiprtext">The predict() function allows us to use the regression model that was obtained from glm() to predict the probability that $Y_i = 1$ for a given $X_i$.</span>
</span><span class="tooltipr">
YourGlmName,&nbsp;
  <span class="tooltiprtext">`YourGlmName` is the name of the object you created when you performed your logistic regression using glm().</span>
</span><span class="tooltipr">
newdata =&nbsp;
<span class="tooltiprtext">The `newdata = ` command allows you to specify the x-values for which you want to obtain predicted probabilities that $Y_i=1$.</span>
</span><span class="tooltipr">
NewDataFrame,&nbsp;
<span class="tooltiprtext">Typically, NewDataFrame is created in real time using the data.frame( X1 = c(Value 1, Value 2, ...), X2 = c(Value 1, Value 2, ...), ...) command. You should see the GSS example file for an example of how to use this function.</span>
</span><span class="tooltipr">
type = "response")
  <span class="tooltiprtext">The type = "response" options specifies that you want predicted probabilities. There are other options available. See ?predict.glm for details.</span>
</span>
</div>
</a>

<div id="MultipleLogisticPredict" style="display:none;">
```{r, echo=FALSE, warning=FALSE}
myglm <- glm(weight>100 ~ Time*Diet, data = ChickWeight, family = binomial)
predict(myglm, newdata = data.frame(Time = 13, Diet = "2"), type = "response")
```
</div>

<br/>

#### Plot the Regression

<div class="tab">
  <button class="tablinks" onclick="openTab(event, 'MultipleLogisticBaseplot')">Base R</button>
  <button class="tablinks" onclick="openTab(event, 'MultipleLogisticggplot')">ggplot2</button>
</div>

<div id="MultipleLogisticBaseplot" class="tabcontent">
  <p>
<a href="javascript:showhide('MultipleLogisticBaseplot2')">
<div class="hoverchunk"><span class="tooltipr">
b <- coef(myglm)
  <span class="tooltiprtext">This stores the estimated coefficients from the regression into the "b" vector.</span>
</span><br/><span class="tooltipr">
palette(
  <span class="tooltiprtext">The palette() function allows you to specify the colors R chooses for the plot.</span>
</span><span class="tooltipr">
c("SomeColor","DifferentColor",...)
  <span class="tooltiprtext">Specify as many colors for the palette as you wish. R will choose a color for each group that you specify later on in your plot(...) code.</span>
</span><span class="tooltipr">
))
  <span class="tooltiprtext">Closing parentheses.</span>
</span><br/><span class="tooltipr">
plot(
  <span class="tooltiprtext">Create the binary scatterplot for the logistic regression.</span>
</span><span class="tooltipr">
Y
  <span class="tooltiprtext">Note that Y must be binary (0,1) values. If it is not, use a logical statement to make it binary, like height>60 or sex=="B".</span>
</span><span class="tooltipr">
&nbsp;~&nbsp;
  <span class="tooltiprtext">The formula operator in R.</span>
</span><span class="tooltipr">
X1,
  <span class="tooltiprtext">The first x-variable in your glm code.</span>
</span><span class="tooltipr">
&nbsp;data = YourDataSet,
  <span class="tooltiprtext">Specify the name of your data set.</span>
</span><span class="tooltipr">
&nbsp;pch = 16)
  <span class="tooltiprtext">Select the type of plotting characters to use for the plot.</span>
</span><br/><span class="tooltipr">
curve(
  <span class="tooltiprtext">Use the curve function to add the logistic regression curve to the plot.</span>
</span><span class="tooltipr">
exp(
  <span class="tooltiprtext">The exp(...) function computes e^(stuff) in R and stands for the "exponential function."</span>
</span><span class="tooltipr">
(b[1]+
  <span class="tooltiprtext">The first coefficient found in b is the y-intercept.</span>
</span><span class="tooltipr">
b[2]
  <span class="tooltiprtext">The second coefficient found in b is the slope term.</span>
</span><span class="tooltipr">
\*x)
  <span class="tooltiprtext">The curve(...) function requires that you call the x-variable "x" although you can change this behavior using xname="SomeOtherName" if you want.</span>
</span><span class="tooltipr">
/
  <span class="tooltiprtext">The logistic model is e^(stuff) / (1 + e^(stuff)).</span>
</span><span class="tooltipr">
(1+exp(b[1]+b[2]*x))
  <span class="tooltiprtext">The denominator of the logistic regression. Be careful to group the entire denominator (1+exp(stuff)).</span>
</span><span class="tooltipr">
&nbsp;col = palette()[1],
  <span class="tooltiprtext">Pull the first color from the color palette for the first curve.</span>
</span><span class="tooltipr">
&nbsp;add = TRUE)
  <span class="tooltiprtext">Add the curve to the current plot.</span>
</span><br/><span class="tooltipr">
curve(exp((b[1]+b[3])+(b[2]+b[4])*x)
  <span class="tooltiprtext">Note how the numerator of this second logistic curve uses and adjusted y-intercept (b[1]+b[3]) and an adjusted slope (b[2]+b[4])*x.</span>
</span><span class="tooltipr">
/
  <span class="tooltiprtext">Begin the denominator.</span>
</span><span class="tooltipr">
(1+exp((b[1]+b[3])+(b[2]+b[4])*x))
  <span class="tooltiprtext">The denominator consists of (1+exp(stuff)).</span>
</span><span class="tooltipr">
,
  <span class="tooltiprtext">Begin optional parameters.</span>
</span><span class="tooltipr">
&nbsp;col = palette()[2],
  <span class="tooltiprtext">Set the color of this second curve to the second color in the palette.</span>
</span><span class="tooltipr">
&nbsp;add = TRUE)
  <span class="tooltiprtext">Add the curve to the current plot.</span>
</span><br/><span class="tooltipr">
legend(
  <span class="tooltiprtext">The legend function adds a legend to the current plot.</span>
</span><span class="tooltipr">
"topright",
  <span class="tooltiprtext">Typically the legend is placed in the top-right corner of the graph. But it could be placed "top", "topleft", "left", "bottomleft", "bottom", "center", "right", or "bottomright".</span>
</span><span class="tooltipr">
&nbsp;legend =
  <span class="tooltiprtext">Why you have to write legend again, no one knows, but you do.</span>
</span><span class="tooltipr">
c("Lable 1", "Label 2"),
  <span class="tooltiprtext">These are the values that will appear in the legend, one entry on each line of the legend.</span>
</span><span class="tooltipr">
&nbsp;col = palette(),
  <span class="tooltiprtext">This specifies the colors of the symbols in the legend. First color goes with "Label 1", and second goes with "Label 2" and so on.</span>
</span><span class="tooltipr">
&nbsp;lty = 1,
  <span class="tooltiprtext">Specifies that solid lines should be used in the legend. If you wanted dots instead, use pch=16. If you wanted dashed lines, use lty=2. Many other options exist.</span>
</span><span class="tooltipr">
&nbsp;bty =
  <span class="tooltiprtext">Specifies the "box type" (bty) that should be drawn around the legend.</span>
</span><span class="tooltipr">
'n')
  <span class="tooltiprtext">By placing a "no" option ('n') here, no box is drawn around the legend.</span>
</span>
</div>
</a>

<div id="MultipleLogisticBaseplot2" style="display:none;">
```{r, warning=FALSE}
myglm <- glm(weight>100 ~ Time*Diet, data = ChickWeight, family = binomial)
palette(c("firebrick", "green", "skyblue", "darkorchid4"))
b <- myglm$coefficients
plot(weight>100 ~ Time, data = ChickWeight, pch=16)
curve(exp(b[1]+b[2]*x)/(1+exp(b[1]+b[2]*x)), add=TRUE, col = palette()[1])
curve(exp((b[1]+b[3])+(b[2]+b[6])*x)/(1+exp((b[1]+b[3])+(b[2]+b[6])*x)), add=TRUE, col = palette()[2])
curve(exp((b[1]+b[4])+(b[2]+b[7])*x)/(1+exp((b[1]+b[4])+(b[2]+b[7])*x)), add=TRUE, col = palette()[3])
curve(exp((b[1]+b[5])+(b[2]+b[8])*x)/(1+exp((b[1]+b[5])+(b[2]+b[8])*x)), add=TRUE, col = palette()[4])
legend("topleft", legend = c("Diet 1", "Diet 2", "Diet 3", "Diet 4"), col = palette(), lty = 1, bty = 'n')
```
</div>

  </p>
</div>

<div id="MultipleLogisticggplot" class="tabcontent">
  <p>
<a href="javascript:showhide('MultipleLogisticggplot2')">
<div class="hoverchunk"><span class="tooltipr">
ggplot(data=YourDataName,
  <span class="tooltiprtext">Initialize a ggplot.</span>
</span><span class="tooltipr">
&nbsp;aes(
  <span class="tooltiprtext">Declare the aesthetics of the graph.</span>
</span><span class="tooltipr">
x=X1,
  <span class="tooltiprtext">Set the x-axis to be your X1 variable.</span>
</span><span class="tooltipr">
y=Y
  <span class="tooltiprtext">Set the y-axis to be your Y-variable.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for aes(...) function.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for ggplot(...) function.</span>
</span><span class="tooltipr">
&nbsp;+
  <span class="tooltiprtext">Add a layer to the ggplot.</span>
</span><br/><span class="tooltipr">
&nbsp; geom_point(
  <span class="tooltiprtext">Add dots to the ggplot.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for the geom_point() function.</span>
</span><span class="tooltipr">
&nbsp;+
  <span class="tooltiprtext">Add a layer to the ggplot.</span>
</span><br/><span class="tooltipr">
&nbsp; geom_smooth(method="glm", method.args=list(family="binomial"), se=FALSE,
  <span class="tooltiprtext">Add the logistic regression curve to the plot.</span>
</span><span class="tooltipr">
&nbsp; aes(color="X2)
  <span class="tooltiprtext">Create a different logistic regression curve for each group found in X2, and give each curve a different color.</span>
</span><span class="tooltipr">
)
  <span class="tooltiprtext">Closing parenthesis for the geom_smooth() function.</span>
</span><span class="tooltipr">
&nbsp;+
  <span class="tooltiprtext">Add a layer to the plot.</span>
</span><br/><span class="tooltipr">
&nbsp; theme_bw()
  <span class="tooltiprtext">Add a black and white them to the plot.</span>
</span>
</div>
</a>

<div id="MultipleLogisticggplot2" style="display:none;">
```{r, warning = FALSE}
ggplot(ChickWeight, aes(Time, as.numeric(weight>100))) +
  geom_point() +
  stat_smooth(aes(color = Diet), method = "glm", method.args = list(family = "binomial"), se = FALSE) +
  theme_bw() +
  labs(y = "weight>100")
```
</div>

  </p>
</div>

----

</div>

### Explanation

<div style="padding-left:125px;">


Multiple Logistic Regression is used when 

* the response variable is binary $(Y_i=0$ or $1)$, and
* there are multiple explanatory variables $X_1,\ldots,X_p$ that can be either quantitative or qualitative. 


#### The Model

Very little changes in multiple logistic regression from Simple Logistic Regression. The probability that $Y_i = 1$ given the observed data $(x_{i1},\ldots,x_{ip})$ is called $\pi_i$ and is modeled by the expanded equation

$$
  P(Y_i = 1|\, x_{i1},\ldots,x_{ip}) = \frac{e^{\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}}}{1+e^{\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip} }} = \pi_i
$$

The assumption is that for certain combinations of $X_1,\ldots,X_p$ the probability that $Y_i=1$ is higher than for other combinations. 


#### Interpretation

The model for $\pi_i$ comes from modeling the *log of the odds* that $Y_i=1$ using a linear regression, i.e., 
$$
  \log\underbrace{\left(\frac{\pi_i}{1-\pi_i}\right)}_{\text{Odds for}\ Y_i=1} = \underbrace{\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}}_{\text{linear regression}}
$$ 
Beginning to solve this equation for $\pi_i$ leads to the intermediate, but important result that
$$
    \underbrace{\frac{\pi_i}{1-\pi_i}}_{\text{Odds for}\ Y_i=1} = e^{\overbrace{\beta_0 + \beta_1 x_{i1} + \ldots + \beta_p x_{ip}}^{\text{liear regression}}} = e^{\beta_0}e^{\beta_1 x_{i1}}\cdots e^{\beta_p x_{ip}}
$$
As in Simple Linear Regression, the values of $e^{\beta_0}$, $e^{\beta_1}$, $\ldots$, $e^{\beta_p}$ are interpreted as the proportional change in odds for $Y_i=1$ when a given $x$-variable experiences a unit change, *all other variables being held constant*. 


#### Checking the Model Assumptions

Diagnostics are the same in multiple logistic regression as they are in simple logistic regression. 

#### Prediction

The idea behind prediction in multiple logistic regression is the same as in simple logistic regression. The only difference is that more than one explanatory variable is used to make the prediction of the *risk* that $Y_i=1$.

----

</div>

### 



----


<footer></footer>