\n\n

Regression models are the workhorse of data science. They are the most\nwell described, practical and theoretically understood models in statistics.\nA data scientist well versed in regression models will be able to\nsolve an incredible array of problems.

\n\n

Perhaps the key insight for regression models is that they produce\nhighly interpretable model fits. This is unlike machine learning\nalgorithms, which often sacrifice interpretability for\nimproved prediction performance or automation. These are, of course,\nvaluable attributes in their own rights. However, the benefit of\nsimplicity, parsimony and interpretability offered by regression models\n(and their close generalizations) should make them a first tool of choice for any practical\nproblem.

\n\n

\nMotivating examples

\n\n

\nFrancis Galton’s height data

\n\n

Francis Galton, the 19th century polymath, can be credited with discovering\nregression. In his landmark paper Regression Toward Mediocrity in\nHereditary Stature\nhe compared the heights of parents and their children. He was particularly\ninterested in the idea that the children of tall parents tended to be tall also,\nbut a little shorter than their parents. Children of short parents tended\nto be short, but not quite as short as their parents. He referred to this\nas “regression to mediocrity” (or regression to the mean). In quantifying\nregression to the mean, he invented what we would call regression.

\n\n

It is perhaps surprising that Galton’s specific work on height is still relevant\ntoday. In fact this European Journal of Human Genetics manuscript compares\nGalton’s prediction models versus those using modern high throughput genomic\ntechnology (spoiler alert, Galton wins).

\n\n

Some questions from Galton’s data come to mind. How would one fit a model that relates\nparent and child heights? How would one predict a child’s height based on their\nparents? How would we quantify regression to the mean? In this class, we’ll\nanswer all of these questions plus many more.

\n\n

\nSimply Statistics versus Kobe Bryant

\n\n

Simply Statistics is a blog by Jeff Leek,\nRoger Peng and Rafael Irizarry. It is one of the most widely read statistics\nblogs, written by three of the top statisticians in academics. Rafa wrote\na (somewhat tongue in cheek) post regarding ball hogging among NBA basketball\nplayers. (By the way, your author has played basketball with Rafael, who\n is quite good, but certainly doesn’t pass up shots; glass\n houses and whatnot.)

\n\n

Here’s some key sentences:

\n\n

“Data supports the claim that if Kobe stops ball hogging the Lakers will win more”
“Linear regression suggests that an increase of 1% in % of shots taken by Kobe results in a drop of 1.16 points (+/- 0.22) in score differential.”

\n\n

In this book we will cover how to create summary statements like this using\nregression model building. Note the nice interpretability of\nthe linear regression model. With this model Rafa numerically relates the\nimpact of more shots taken on score differential.

\n\n

\nSummary notes: questions for this book

Regression models are incredibly handy statistical tools. One can use\nthem to answer all sorts of questions.\nConsider three of the most common tasks for regression models:

\n\n

\nPrediction e.g.: to use the parent’s heights to predict children’s heights.
\nModeling e.g.: to try to find a parsimonious, easily described mean\n relationship between parental and child heights.
\nCovariation e.g.: to investigate the variation in child heights that appears\n unrelated to parental heights (residual variation) and to quantify what impact genotype information has beyond parental height in explaining child height.

\n\n

An important aspect, especially in questions 2 and 3 is assessing\nmodeling assumptions. For example, it is important to figure out how/whether\nand what assumptions are needed to generalize findings beyond the data in question.\nPresumably, if we find a relationship between parental and child heights,\nwe’d like to extend that knowledge beyond the data used to build the model.\nThis requires assumptions. In this book, we’ll cover the main assumptions\nnecessary.

\n\n

\nExploratory analysis of Galton’s Data

\n\n

Let’s look at the data first. This data was created by Francis Galton in 1885.\nGalton was a statistician who invented the term and concepts\nof regression and correlation, founded the journal Biometrika, and\nwas the cousin of Charles Darwin.

\n\n

You may need to run install.packages(\"UsingR\") if the UsingR library is\nnot installed. Let’s look at the marginal (parents disregarding children and\nchildren disregarding parents) distributions first.\nThe parental distribution is all heterosexual couples. The parental average was corrected\nfor gender via multiplying female heights by 1.08. Remember, Galton didn’t have\nregression to help figure out a better way to do this correction!

\n\n

Loading and plotting Galton’s data.

library(UsingR); data(galton); library(reshape); long <- melt(galton)\ng <- ggplot(long, aes(x = value, fill = variable))\ng <- g + geom_histogram(colour = \"black\", binwidth=1)\ng <- g + facet_grid(. ~ variable)\ng\n

\n\n\n\n\n

\"Plotting — Plotting the `galton` dataset

\n\n\n

\nFinding the middle via least squares

Consider only the children’s heights. How could one describe the “middle”?\nConsider one definition. Let \$Y_i\$ be the height of child \$i\$ for\n\$i = 1, \\ldots, n = 928\$, then define the middle as the value\nof \$\\mu\$ that minimizes

\n\n\n$$\n\\sum_{i=1}^n (Y_i - \\mu)^2.\n$$\n\n

This is physical center of mass of the histogram.\nYou might have guessed that the answer \$\\mu = \\bar Y\$. This is called\nthe least squares estimate for \$\\mu\$. It is the point that minimizes\nthe sum of the squared distances between the observed data and itself.

\n\n

Note, if there was no variation in the data, every value of \$Y_i\$ was\nthe same, then there would be no error around the mean. Otherwise, our estimate\nhas to balance the fact that our estimate of \$\\mu\$ isn’t going to\npredict every observation perfectly. Minimizing the average (or sum of the)\nsquared errors seems like a reasonable strategy, though of course there are others.\nWe could minimize the average absolute deviation between the data \$\\mu\$\n (this leads to the median as the estimate instead of the mean).\nHowever, minimizing the squared error has many nice properties, so we’ll stick\nwith that for this class.

\n\n

\nExperiment

Let’s\nuse RStudio’s manipulate to see what value of\n\$\\mu\$ minimizes the sum of the squared deviations. The code below\nallows you to create a slider to investigate estimates and their\nmean squared error.

\n\n

Using manipulate to find the least squares estimate.

library(manipulate)\nmyHist <- function(mu){\n    mse <- mean((galton$child - mu)^2)\n    g <- ggplot(galton, aes(x = child)) + geom_histogram(fill = \"salmon\", colour = \"\\\nblack\", binwidth=1)\n    g <- g + geom_vline(xintercept = mu, size = 3)\n    g <- g + ggtitle(paste(\"mu = \", mu, \", MSE = \", round(mse, 2), sep = \"\"))\n    g\n}\nmanipulate(myHist(mu), mu = slider(62, 74, step = 0.5))\n

\n\n\n\n

The least squares estimate is the empirical mean.

\n\n

g <- ggplot(galton, aes(x = child)) + geom_histogram(fill = \"salmon\", colour = \"blac\\\nk\", binwidth=1)\ng <- g + geom_vline(xintercept = mean(galton$child), size = 3)\ng\n

\n\n\n\n\n

\"The — The best mean is the vertical line.

\n\n\n

\nThe math (not required)

\n\n

Why is the sample average the least squares estimate for \$\\mu\$?\nIt’s surprisingly easy to show. Perhaps more surprising is how generally\nthese results can be extended.

\n\n\n$$\n\\begin{eqnarray*}\n\\sum_{i=1}^n (Y_i - \\mu)^2 & = & \\sum_{i=1}^n (Y_i - \\bar Y + \\bar Y - \\mu)^2 \\\\\n & = & \\sum_{i=1}^n (Y_i - \\bar Y)^2 + 2 \\sum_{i=1}^n (Y_i - \\bar Y) (\\bar Y - \\mu) + \\sum_{i=1}^n (\\bar Y - \\mu)^2 \\\\\n & = & \\sum_{i=1}^n (Y_i - \\bar Y)^2 + 2 (\\bar Y - \\mu) \\sum_{i=1}^n (Y_i - \\bar Y) + \\sum_{i=1}^n (\\bar Y - \\mu)^2 \\\\\n & = & \\sum_{i=1}^n (Y_i - \\bar Y)^2 + 2 (\\bar Y - \\mu) (\\sum_{i=1}^n Y_i - n \\bar Y) + \\sum_{i=1}^n (\\bar Y - \\mu)^2 \\\\\n & = & \\sum_{i=1}^n (Y_i - \\bar Y)^2 + \\sum_{i=1}^n (\\bar Y - \\mu)^2\\\\\n & \\geq & \\sum_{i=1}^n (Y_i - \\bar Y)^2 \\\n\\end{eqnarray*}\n$$\n\n

\nComparing children’s heights and their parent’s heights

\n\n

Looking at either the parents or children on their own isn’t interesting.\nWe’re interested in how the relate to each other. Let’s plot the parent\nand child heights.

\n\n

ggplot(galton, aes(x = parent, y = child)) + geom_point()\n

\n\n\n\n\n

\n\n\n

The overplotting is clearly hiding some data. Here you can get the code\n\nto make the size and color of the points be the frequency.

\n\n\n

\n\n\n

\nRegression through the origin

A line requires two parameters to be specified, the intercept and the slope.\nLet’s first focus on the slope. We want to find the slope of the line that\nbest fits the data. However, we have to pick a good intercept. Let’s subtract\nthe mean from both the parent and child heights so that their subsequent\nmeans are 0.\nNow let’s find the line that goes through the origin (has intercept 0) by\npicking the best slope.

\n\n

Suppose that \$X_i\$ are the parent heights with the mean subtracted.\nConsider picking the slope \$\\beta\$ that minimizes

\n\n\n$$\n\\sum_{i=1}^n (Y_i - X_i \\beta)^2.\n$$\n\n

Each \$X_i \\beta\$ is the vertical height of\na line through the origin at point \$X_i\$.\nThus, \$Y_i - X_i \\beta\$\nis the vertical distance between the line\nat each observed \$X_i\$ point (parental height) and the\n\$Y_i\$ (child height).

\n\n

Our goal is exactly to use the origin as a pivot point and pick the\nline that minimizes the sum of the squared vertical distances\nof the points to the line. Use RStudio’s manipulate function to experiment. \nSubtract the means so that the origin is the mean of the parent\nand children heights.

\n\n

Code for plotting the data.

library(dplyr)\ny <- galton$child - mean(galton$child)\nx <- galton$parent - mean(galton$parent)\nfreqData <- as.data.frame(table(x, y))\nnames(freqData) <- c(\"child\", \"parent\", \"freq\")\nfreqData$child <- as.numeric(as.character(freqData$child))\nfreqData$parent <- as.numeric(as.character(freqData$parent))\nmyPlot <- function(beta){\n    g <- ggplot(filter(freqData, freq > 0), aes(x = parent, y = child))\n    g <- g  + scale_size(range = c(2, 20), guide = \"none\" )\n    g <- g + geom_point(colour=\"grey50\", aes(size = freq+20), show.legend = FALSE)\n    g <- g + geom_point(aes(colour=freq, size = freq))\n    g <- g + scale_colour_gradient(low = \"lightblue\", high=\"white\")\n    g <- g + geom_abline(intercept = 0, slope = beta, size = 3)\n    mse <- mean( (y - beta * x) ^2 )\n    g <- g + ggtitle(paste(\"beta = \", beta, \"mse = \", round(mse, 3)))\n    g\n}\nmanipulate(myPlot(beta), beta = slider(0.6, 1.2, step = 0.02))\n

\n\n\n\n

\nThe solution

In the next few lectures we’ll talk about why this is the solution. But,\nrather than leave you hanging, here it is:

\n\n

> lm(I(child - mean(child))~ I(parent - mean(parent)) - 1, data = galton)\nCall:\nlm(formula = I(child - mean(child)) ~ I(parent - mean(parent)) -\n    1, data = galton)\n\nCoefficients:\nI(parent - mean(parent))  \n                   0.646  \n

\n\n\n\n

Let’s plot the best fitting line. In the subsequent chapter we will learn all\nabout creating, interpreting and performing inference on such mode fits. (Note\nthat I shifted the origin back to the means of the original data.) The\nresults suggest that for every 1 inch increase in the parents’ height,\nwe estimate a 0.646 inch increase in the child’s height.

\n\n\n

\"Data — Data with the best fitting line.

\n\n\n

\nExercises

Consider the dataset given by x=c(0.725,0.429,-0.372 ,0.863). What value\nof mu minimizes sum((x - mu)^2)? Watch a video solution.\n
Reconsider the previous question. Suppose that weights were given,\nw = c(2, 2, 1, 1) so that we wanted to minimize sum(w * (x - mu) ^ 2) for mu.\nWhat value would we obtain? Watch a video solution.\n
Take the Galton dataset and obtain the regression through the origin slope estimate\nwhere the centered parental height is the outcome and the child’s height is\nthe predictor. Watch a video solution.\n

\n\n\n

\nNotation

\n\n

\nSome basic definitions

\n\n

In this chapter, we’ll cover some basic definitions and notation used throughout the\nbook. We will try to minimize the amount of mathematics required so that we\ncan focus on the concepts.

\n\n

\nNotation for data

We write \$X_1, X_2, \\ldots, X_n\$ to describe \$n\$ data points.\nAs an example, consider the data set \$\\{1, 2, 5\\}\$ then\n\$X_1 = 1\$, \$X_2 = 2\$, \$X_3 = 5\$ and \$n = 3\$.

\n\n

Of course, there’s nothing in particular about the variable \$X\$.\nWe often use a different letter, such as \$Y_1, \\ldots , Y_n\$ to describe\na data set. We will typically use Greek letters for things we don’t know.\nSuch as, \$\\mu\$ being a population mean that we’d like to estimate.

\n\n

\nThe empirical mean

\n\n

The empirical mean is a measure of center of our data. Under sampling assumptions,\nit estimates a population mean of interest. Define the empirical mean as

\n\n\n$$\n\\bar X = \\frac{1}{n}\\sum_{i=1}^n X_i.\n$$\n\n

Notice if we subtract the mean from data points, we get data that has mean 0. That is, if we define

\n\n\n$$\n\\tilde X_i = X_i - \\bar X.\n$$\n\n

then the mean of the \$\\tilde X_i\$ is 0.\nThis process is called centering the random variables.\nRecall from the previous lecture that the empirical mean is\nthe least squares solution for minimizing

\n\n\n$$\n\\sum_{i=1}^n (X_i - \\mu)^2\n$$\n\n

\nThe empirical standard deviation and variance

\n\n

The variance and standard deviation are measures of how spread out our data is.\nUnder sampling assumptions, they estimate variability in the population.\nWe define the empirical variance as:

\n\n\n$$\nS^2 = \\frac{1}{n-1} \\sum_{i=1}^n (X_i - \\bar X)^2\n= \\frac{1}{n-1} \\left( \\sum_{i=1}^n X_i^2 - n \\bar X ^ 2 \\right)\n$$\n\n

The empirical standard deviation is defined as\n\$S = \\sqrt{S^2}.\$

\n\n

Notice that the standard deviation has the same units as the data.\nThe data defined by \$X_i / s\$ have empirical standard deviation 1.\nThis is called scaling the data.

\n\n

\nNormalization

\n\n

We can combine centering and scaling of data as follows to get normalized data.\nIn particular, the data defined by:

\n\n\n$$\nZ_i = \\frac{X_i - \\bar X}{s}\n$$\n\n

has empirical mean zero and empirical standard deviation 1.\nThe process of centering then scaling the data is called normalizing the data.\nNormalized data are centered at 0 and have units equal to standard deviations of the original data.\nExample, a value of 2 from normalized data means that data point was two standard deviations larger than the mean.

\n\n

Normalization is very useful for creating data that comparable across experiments\nby getting rid of any shifting or scaling effects.

\n\n

\nThe empirical covariance

This class is largely considering how variables covary. This is estimated\nby the empirical covariance.\nConsider now when we have pairs of data, \$(X_i, Y_i)\$.\nTheir empirical covariance is defined as:

\n\n\n$$\nCov(X, Y) =\n\\frac{1}{n-1}\\sum_{i=1}^n (X_i - \\bar X) (Y_i - \\bar Y)\n= \\frac{1}{n-1}\\left( \\sum_{i=1}^n X_i Y_i - n \\bar X \\bar Y\\right)\n$$\n\n

This measure is of limited utility, since its units are the product of\nthe units of the two variables. A more useful definition normalizes the two\nvariables first.

\n\n

The correlation is defined as:

\n\n\n$$\nCor(X, Y) = \\frac{Cov(X, Y)}{S_x S_y}\n$$\n\n

where \$S_x\$ and \$S_y\$ are the estimates of standard deviations\nfor the \$X\$ observations and \$Y\$ observations, respectively. The correlation is simply the covariance of\nthe separately normalized X and Y data. Because the \ndata have been normalized, the correlation is a unit free\nquantity and thus has more of a hope of being interpretable\nacross settings.

\n\n

\nSome facts about correlation

First, the order of the arguments is irrelevant \$Cor(X, Y) = Cor(Y, X)\$\nSecondly, it has to be between -1 and 1, \$-1 \\leq Cor(X, Y) \\leq 1\$.\nThirdly, the correlation is exactly -1 or 1 only when the observations fall\nperfectly on a negatively or positively sloped, line, respectively.\nFourthly, \$Cor(X, Y)\$ measures the strength of the linear relationship between\nthe two variables, with stronger relationships as \$Cor(X,Y)\$ heads towards -1 or 1.\nFinally, \$Cor(X, Y) = 0\$ implies no linear relationship.

\n\n

\nExercises

Take the Galton dataset and find the mean, standard deviation and correlation between\nthe parental and child heights. Watch a video solution.\n
Center the parent and child variables and verify that the centered variable means are 0. Watch a video solution.\n
Rescale the parent and child variables and verify that the scaled variable standard deviations are 1. Watch a video solution.\n
Normalize the parental and child heights. Verify that\nthe normalized variables have mean 0 and standard deviation 1 and take the\ncorrelation between them. Watch a video solution.\n

\n\n\n

\nOrdinary least squares

\n\n

Ordinary least squares (OLS) is the workhorse of statistics. It gives\na way of taking complicated outcomes and explaining behavior (such as trends)\nusing linearity. The simplest application of OLS is fitting a line.

\n\n

\nGeneral least squares for linear equations

Consider again the parent and child height data from Galton.

\n\n\n

Let’s try fitting the best line.\nLet \$Y_i\$ be the \$i^{th}\$ child’s height and\n\$X_i\$ be the \$i^{th}\$ (average over the pair of) parental heights.\nConsider finding the best line of the form

\n\n\n$$\n\\mbox{Child Height} = \\beta_0 + \\mbox{Parent Height} \\beta_1,\n$$\n\n

Let’s try using least squares by minimizing the following equation\nover \$\\beta_0\$ and \$\\beta_1\$:

\n\n\n$$\n\\sum_{i=1}^n \\{Y_i - (\\beta_0 + \\beta_1 X_i)\\}^2.\n$$\n\n

Minimizing this equation will minimize the sum of the squared distances\nbetween the fitted line at the parents’ heights (\$ \\beta_1 X_i\$)\nand the observed child heights (\$Y_i\$).

\n\n

The result actually has a closed form. Specifically, the least squares\nof the line:

\n\n\n$$\nY = \\beta_0 + \\beta_1 X,\n$$\n\n

through the data pairs \$(X_i, Y_i)\$ with \$Y_i\$ as the outcome obtains\nthe line \$Y = \\hat \\beta_0 + \\hat \\beta_1 X\$\nwhere:

\n\n\n$$\n\\hat \\beta_1 = Cor(Y, X) \\frac{Sd(Y)}{Sd(X)} ~~~\\mbox{and}~~~ \\hat \\beta_0 = \\bar Y - \\hat \\beta_1 \\bar X.\n$$\n\n

At this point, a couple of notes are in order.\nFirst, the slope, \$\\hat \\beta_1\$, has the units of \$Y / X\$. Secondly,\nthe intercept, \$\\hat \\beta_0\$, has the units of \$Y\$.

\n\n

The line passes through the point \$(\\bar X, \\bar Y)\$. If you center your\nXs and Ys first, then the line will pass through the origin. Moreover,\nthe slope is the same one you would get if you centered the data,\n\$(X_i - \\bar X, Y_i - \\bar Y)\$, and either fit a linear regression\nor regression through the origin.

\n\n

To elaborate, regression through the origin, assuming that \$\\beta_0 = 0\$,\nyields the following solution to the least squares criteria:

\n\n\n$$\n\\hat \\beta_1 = \\frac{\\sum_{i=1}^n X_i Y_i}{\\sum_{i=1}^n X_i^2},\n$$\n\n

This is exactly the correlation times the\nratio in the standard deviations if\nthe both the Xs and Ys have been centered first. (Try it out using R to\nverify this!)

\n\n

It is interesting to think about what happens when you reverse the role of \$X\$\nand \$Y\$. Specifically, the slope of the regression line with \$X\$\nas the outcome and \$Y\$ as the predictor is \$Cor(Y, X) Sd(X)/ Sd(Y)\$.

\n\n

If you normalized the data,\n\$\\{ \\frac{X_i - \\bar X}{Sd(X)}, \\frac{Y_i - \\bar Y}{Sd(Y)}\\}\$,\nthe slope is simply the correlation, \$Cor(Y, X)\$, regardless of which variable\nis treated as the outcome.

\n\n

\nRevisiting Galton’s data

\n\n

Let’s double check our calculations using R

\n\n

Fitting Galton’s data using linear regression.

> y <- galton$child\n> x <- galton$parent\n> beta1 <- cor(y, x) *  sd(y) / sd(x)\n> beta0 <- mean(y) - beta1 * mean(x)\n> rbind(c(beta0, beta1), coef(lm(y ~ x)))\n     (Intercept)      x\n[1,]       23.94 0.6463\n[2,]       23.94 0.6463\n

\n\n\n\n

We can see that the result of lm is identical to hard coding the fit ourselves.\nLet’s reverse the outcome/predictor relationship.

\n\n

> beta1 <- cor(y, x) *  sd(x) / sd(y)\n> beta0 <- mean(x) - beta1 * mean(y)\n> rbind(c(beta0, beta1), coef(lm(x ~ y)))\n     (Intercept)      y\n[1,]       46.14 0.3256\n[2,]       46.14 0.3256\n

\n\n\n\n

Now let’s show that\nregression through the origin yields an equivalent slope if you center the data first

\n\n

> yc <- y - mean(y)\n> xc <- x - mean(x)\n> beta1 <- sum(yc * xc) / sum(xc ^ 2)\nc(beta1, coef(lm(y ~ x))[2])\n            x\n0.6463 0.6463\n

\n\n\n\n

Now let’s show that normalizing variables results in the slope being the correlation.

\n\n

> yn <- (y - mean(y))/sd(y)\n> xn <- (x - mean(x))/sd(x)\n> c(cor(y, x), cor(yn, xn), coef(lm(yn ~ xn))[2])\n                  xn\n0.4588 0.4588 0.4588\n

\n\n\n\n

The image below plots the data again, the best fitting line and standard\nerror bars for the fit. We’ll work up to that point later. But, understanding\nthat our fitted line is estimated with error is an important concept.\nYou can find the code for the plot here.

\n\n\n

\nShowing the OLS result

If you would like to see a proof of why the ordinary least squares result works out to be the way that it is:\nwatch this video.

\n\n

\nExercises

Install and load the package UsingR and load the father.son data with data(father.son). Get the\nlinear regression fit where the son’s height is the outcome and the father’s height is the predictor.\nGive the intercept and the slope, plot the data and overlay the fitted regression line. Watch a video solution.\n
Refer to problem 1. Center the father and son variables and refit the model omitting the intercept. Verify that the slope estimate is the same as the linear regression fit\nfrom problem 1. Watch a video solution.\n
Refer to problem 1. Normalize the father and son data and see that the fitted slope is the correlation. Watch a video solution.\n
Go back to the linear regression line from Problem 1. If a father’s height was 63 inches, what would you predict\nthe son’s height to be? Watch a video solution.\n
Consider a data set where the standard deviation of the outcome variable is double that of the predictor. Also, the variables\nhave a correlation of 0.3. If you fit a linear regression model, what would be the estimate of the slope? Watch a video solution.\n
Consider the previous problem. The outcome variable has a mean of 1 and the predictor has a mean of 0.5. What would be\nthe intercept? Watch a video solution.\n
True or false, if the predictor variable has mean 0, the estimated intercept from linear regression\nwill be the mean of the outcome? Watch a video solution.\n
Consider problem 5 again. What would be the estimated slope if the predictor and outcome were reversed? Watch a video solution.\n

\n\n\n

\nRegression to the mean

\n\n

\nA historically famous idea, regression to the mean

Here is a fundamental question.\nWhy is it that the children of tall parents tend to be tall, but not as tall as their parents?\nWhy do children of short parents tend to be short, but not as short as their parents?\nConversely, why do parents of very short children, tend to be short, but not a short as their child? And the same with parents of very tall children?

\n\n

We can try this with anything that is measured with error.\nWhy do the best performing athletes this year tend to do a little worse the following?\nWhy do the best performers on hard exams always do a little worse on the next hard exam?

\n\n

These phenomena are all examples of so-called regression to the mean.\nRegression to the mean, was invented by Francis Galton in the paper “Regression towards mediocrity in hereditary stature” The Journal of the Anthropological Institute of Great Britain and Ireland , Vol. 15, (1886). The idea served as a foundation for the discovery of linear regression.

\n\n

Think of it this way, imagine if you simulated pairs of random normals.\nThe largest first ones would be the largest by chance, and the probability that\nthere are smaller for the second simulation is high.\nIn other words \$P(Y < x | X = x)\$ gets bigger as \$x\$\nheads to the very large values.\nSimilarly \$P(Y > x | X = x)\$ gets bigger as \$x\$ heads to very\nsmall values. Think of the regression line as the intrinsic part and the\nregression to the mean as the result of noise.\nUnless \$Cor(Y, X) = 1\$ the intrinsic part isn’t perfect and so we\nshould think about how much regression to the mean should occur. In other words,\nwhat should we multiply tall parent’s heights by to predict their children’s\nheight?

\n\n

\nRegression to the mean

Let’s investigate this with Galton’s father and son data.\nSuppose that we normalize \$X\$ (child’s height) and \$Y\$\n(father’s height) so that they both have mean 0 and variance 1.\nThen, recall, our regression line passes through \$(0, 0)\$ (the mean of the X and Y).\nThe slope of the regression line is \$Cor(Y,X)\$, regardless of which variable is the outcome (recall, both standard deviations are 1).\nNotice if \$X\$ is the outcome and you create a plot where \$X\$\nis the horizontal axis, the slope of the least squares line that you plot is\n\$1/Cor(Y, X)\$. Let’s plot the normalized\nfather and son heights to investigate.

\n\n

Code for the plot.

library(UsingR)\ndata(father.son)\ny <- (father.son$sheight - mean(father.son$sheight)) / sd(father.son$sheight)\nx <- (father.son$fheight - mean(father.son$fheight)) / sd(father.son$fheight)\nrho <- cor(x, y)\nlibrary(ggplot2)\ng = ggplot(data.frame(x, y), aes(x = x, y = y))\ng = g + geom_point(size = 5, alpha = .2, colour = \"black\")\ng = g + geom_point(size = 4, alpha = .2, colour = \"red\")\ng = g + geom_vline(xintercept = 0)\ng = g + geom_hline(yintercept = 0)\ng = g + geom_abline(position = \"identity\")\ng = g + geom_abline(intercept = 0, slope = rho, size = 2)\ng = g + geom_abline(intercept = 0, slope = 1 / rho, size = 2)\ng = g + xlab(\"Father's height, normalized\")\ng = g + ylab(\"Son's height, normalized\")\ng\n

\n\n\n\n\n

\n\n\n

Let’s investigate the plot and the regression fits.\nIf you had to predict a son’s normalized height, it would be\n\$Cor(Y, X) * X_i\$ where \$X_i\$ was the normalized\nfather’s height. Conversely,\nif you had to predict a father’s normalized height, it would be\n\$Cor(Y, X) * Y_i\$.

\n\n

Multiplication by this correlation shrinks toward 0 (regression toward the mean).\nIt is in this way that Galton used regression to account for regression toward\nthe mean.\nIf the correlation is 1 there is no regression to the mean,\n(if father’s height perfectly determines child’s height and vice versa).

\n\n

Note since Galton’s original seminal paper, the idea of regression to the mean\nhas been generalized and expanded upon. However, the core remains. In\npaired measurements, if there’s randomness then the extreme values of\none element of the pair will be likely less extreme in the other element.

\n\n

The number of applications of this phenomena is staggering. Some financial\nadvisors recommend putting your money in your worst performing fund because\nof regression to the mean. (If there’s a lot of noise, those are the most\nlikely to gain in value.) An example that I’ve run into is that students performing\nthe best on midterm exams often do much worse on the final.\nAthletes often follow a phenomenal season with merely\na good season. It’s a useful exercise\nto think whenever paired observations are being evaluated whether real intrinsic\nproperties are being discussed, or just regression to the mean.

\n\n

\nExercises

You have two noisy scales and a bunch of people that you’d like to weigh. You\nweigh each person on both scales. The correlation was 0.75. If you normalized\neach set of weights, what would you have to multiply the weight on one scale\nto get a good estimate of the weight on the other scale? Watch a video solution.\n
Consider the previous problem. Someone’s weight was 2 standard deviations above\nthe mean of the group on the first scale. How many standard deviations above the\nmean would you estimate them to be on the second? Watch a video solution.\n
You ask a collection of husbands and wives to guess how many jellybeans are in a jar.\nThe correlation is 0.2. The standard deviation for the husbands is 10 beans while the\nstandard deviation for wives is 8 beans. Assume that the data were centered so\nthat 0 is the mean for each. The centered guess for a husband was 30 beans (above\n the mean). What would be your best estimate of the wife’s guess? Watch a video solution.\n

\n\n\n

\nStatistical linear regression models

\n\n

\n\n

Up to this point, we’ve only considered estimation. Estimation is useful,\nbut we also need to know how to extend our estimates to a population.\nThis is the process of statistical inference. Our approach to statistical\ninference will be through a statistical model. At the bare minimum, we\nneed a few distributional assumptions on the errors. However, we’ll focus on\nfull model assumptions under Gaussianity.

\n\n

\nBasic regression model with additive Gaussian errors.

Consider developing a probabilistic model for linear regression. Our\nstarting point will assume a systematic component via a line and then\nindependent and identically distributed Gaussian errors. We can write\nthe model out as:

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X_i + \\epsilon_{i}\n$$\n\n

Here, the \$\\epsilon_{i}\$ are assumed to be independent and\nidentically distributed as\n\$N(0, \\sigma^2)\$. Under this model,

\n\n\n$$\nE[Y_i ~|~ X_i = x_i] = \\mu_i = \\beta_0 + \\beta_1 x_i\n$$\n\n

and

\n\n\n$$\nVar(Y_i ~|~ X_i = x_i) = \\sigma^2.\n$$\n\n

This model implies\n that the \$Y_i\$ are independent and normally\ndistributed with means \$\\beta_0 + \\beta_1 x_i\$ and variance\n\$\\sigma^2\$. We could write this more compactly as

\n\n\n$$\nY_i ~|~ X_i = x_i \\sim N(\\beta_0 + \\beta_1 x_i, \\sigma^2).\n$$\n\n

While this specification of the model is a perhaps better for advanced\npurposes, specifying the model as linear with additive error terms is\ngenerally more useful. With that specification, we can hypothesize and\ndiscuss the nature of the errors. In fact, we’ll even cover ways to estimate\nthem to investigate our model assumption.

\n\n

Remember that our least squares estimates of\n\$\\beta_0\$ and \$\\beta_1\$ are:

\n\n\n$$\n\\hat \\beta_1 = Cor(Y, X) \\frac{Sd(Y)}{Sd(X)} ~~~\\mbox{and}~~~ \\hat \\beta_0 = \\bar Y - \\hat \\beta_1 \\bar X.\n$$\n\n

It is convenient that under our Gaussian additive error model\nthat the maximum likelihood estimates of\n\$\\beta_0\$ and \$\\beta_1\$ are the least squares estimates.

\n\n

\nInterpreting regression coefficients, the intercept

\n\n

\n\n

Our model allows us to attach statistical interpretations to our parameters.\nLet’s start with the intercept; \$\\beta_0\$ represents\nthe expected value of the response when the predictor is 0. We can show this\nas:

\n\n\n$$\nE[Y | X = 0] = \\beta_0 + \\beta_1 \\times 0 = \\beta_0.\n$$\n\n

Note, the intercept isn’t always of interest. For example,\nwhen \$X=0\$ is impossible or far outside of the range of data.\nTake as a specific instance, when X is blood pressure, no one is interested\nin studying blood pressure’s impact on anything for values near 0.

\n\n

There is a way to make your intercept more interpretable.\nConsider that:

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X_i + \\epsilon_i\n= \\beta_0 + a \\beta_1 + \\beta_1 (X_i - a) + \\epsilon_i\n= \\tilde \\beta_0 + \\beta_1 (X_i - a) + \\epsilon_i.\n$$\n\n

Therefore, shifting your \$X\$ values by value \$a\$\nchanges the intercept, but not the slope.\nOften \$a\$ is set to \$\\bar X\$, so that the intercept is\ninterpreted as the expected response at the average \$X\$ value.

\n\n

\nInterpreting regression coefficients, the slope

Now that we understand how to interpret the intercept, let’s try interpreting\nthe slope. Our slope, \$\\beta_1\$,\nis the expected change in response for a 1 unit change in the predictor.\nWe can show that as follows:

\n\n\n$$\nE[Y ~|~ X = x+1] - E[Y ~|~ X = x] =\n\\beta_0 + \\beta_1 (x + 1) - (\\beta_0 + \\beta_1 x ) = \\beta_1\n$$\n\n

Notice that the interpretation of \$\\beta_1\$ is tied to the\nunits of the X variable. Let’s consider the impact of changing the units.

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X_i + \\epsilon_i\n= \\beta_0 + \\frac{\\beta_1}{a} (X_i a) + \\epsilon_i\n= \\beta_0 + \\tilde \\beta_1 (X_i a) + \\epsilon_i\n$$\n\n

Therefore, multiplication of \$X\$ by a factor \$a\$\nresults in dividing the coefficient by a factor of \$a\$.

\n\n

As an example, suppose that \$X\$ is height in meters (m) and \$Y\$\nis weight in kilograms (kg). Then \$\\beta_1\$ is kg/m.\nConverting \$X\$ to centimeters implies multiplying \$X\$ by 100 cm/m.\nTo get \$\\beta_1\$ in the right units if we had fit the model in meters,\nwe have to divide by 100 cm/m. Or, we can write out the notation as:

\n\n\n$$\nX m \\times \\frac{100cm}{m} = (100 X) cm\n~~\\mbox{and}~~\n\\beta_1 \\frac{kg}{m} \\times\\frac{1 m}{100cm} =\n\\left(\\frac{\\beta_1}{100}\\right)\\frac{kg}{cm}\n$$\n\n

\nUsing regression for prediction

\n\n

\n\n

Regression is quite useful for prediction.\nIf we would like to guess the outcome at a particular\nvalue of the predictor, say \$X\$, the regression model guesses:

\n\n\n$$\n\\hat \\beta_0 + \\hat \\beta_1 X\n$$\n\n

In other words, just find the Y value on the line with the corresponding X\nvalue. Regression, especially linear regression, often doesn’t produce\nthe best prediction algorithms. However, it produces parsimonious and\ninterpretable models along with the predictions. Also, as we’ll see later\nwe’ll be able to get easily described estimates of uncertainty associated\nwith our predictions.

\n\n

\nExample

\n\n

Let’s analyze the diamond data set from the UsingR package.\nThe data is diamond prices (in Singapore dollars) and diamond weight\nin carats. Carats are a standard measure of diamond mass, 0.2 grams.\nTo get the data use library(UsingR); data(diamond)

\n\n

First let’s plot the data. Here’s the code I used

\n\n

library(UsingR)\ndata(diamond)\nlibrary(ggplot2)\ng = ggplot(diamond, aes(x = carat, y = price))\ng = g + xlab(\"Mass (carats)\")\ng = g + ylab(\"Price (SIN $)\")\ng = g + geom_point(size = 7, colour = \"black\", alpha=0.5)\ng = g + geom_point(size = 5, colour = \"blue\", alpha=0.2)\ng = g + geom_smooth(method = \"lm\", colour = \"black\")\ng\n

\n\n\n\n

and here’s the plot.

\n\n\n

First, let’s fit the linear regression model. This is done\nwith the lm function in R (lm stands for linear model). The\nsyntax is lm(Y ~ X) where Y is the response and X is the\npredictor.

\n\n

> fit <- lm(price ~ carat, data = diamond)\n> coef(fit)\n(Intercept)       carat\n     -259.6      3721.0\n

\n\n\n\n

The function coef grabs the fitted coefficients and conveniently names them\nfor you. Therefore, we estimate an expected 3721.02 (SIN) dollar increase in price\nfor every carat increase in mass of diamond.\nThe intercept -259.63 is the expected price of a 0 carat diamond.

\n\n

We’re not interested in 0 carat diamonds (it’s hard to get a good\nprice for them ;-). Let’s fit the model with a more interpretable intercept\nby centering our X variable.

\n\n

> fit2 <- lm(price ~ I(carat - mean(carat)), data = diamond)\ncoef(fit2)\n           (Intercept) I(carat - mean(carat))\n                 500.1                 3721.0\n

\n\n\n\n

Thus the new intercept, 500.1, is the expected price for\nthe average sized diamond of the data (0.2042 carats). Notice\nthe estimated slope didn’t change at all.

\n\n

Now let’s try changing the scale.\nThis is useful since a one carat increase in a diamond is pretty big.\nWhat about changing units to 1/10th of a carat?\nWe can just do this by just dividing the coefficient by 10, no need to\nrefit the model.

\n\n

Thus, we expect a 372.102 (SIN) dollar change in price for every\n1/10th of a carat increase in mass of diamond.

\n\n

Let’s show via R that this is the same as rescaling our X variable\nand refitting. To go from 1 carat to 1/10 of a carat units, we need\nto multiply our data by 10.

\n\n

> fit3 <- lm(price ~ I(carat * 10), data = diamond)\n> coef(fit3)\n  (Intercept) I(carat * 10)\n       -259.6         372.1\n

\n\n\n\n

Now, let’s predict the price of a diamond. This should be as\neasy as just evaluating the fitted line at the price we want to

\n\n

> newx <- c(0.16, 0.27, 0.34)\n> coef(fit)[1] + coef(fit)[2] * newx\n[1]  335.7  745.1 1005.5\n

\n\n\n\n

Therefore, we predict the price to be 335.7, 745.1 and 1005.5 for\na 0.16, 0.26 and 0.34 carat diamond. Of course, our prediction models\nwill get more elaborate and R has a generic function, predict,\nto put our X values into the model for us. The data has to go into\nthe model as a data frame with the same named X variables.

\n\n

> predict(fit, newdata = data.frame(carat = newx))\n     1      2      3\n 335.7  745.1 1005.5\n

\n\n\n\n

Let’s visualize our prediction. In the following plot, the\npredicted values at the observed Xs are the red points on the fitted\nline. The new X values are the at vertical lines, which are connected to\nthe predicted values via the connected horizontal lines.

\n\n\n

\"Illustrating — Illustrating prediction with regression.

\n\n\n

\nExercises

Fit a linear regression model to the father.son dataset with\nthe father as the predictor and the son as the outcome. Give a p-value\nfor the slope coefficient and perform the relevant hypothesis test. Watch a video solution.\n
Refer to question 1. Interpret both parameters. Recenter for the\nintercept if necessary. Watch a video solution.\n
Refer to question 1. Predict the son’s height if the father’s height\nis 80 inches. Would you recommend this prediction? Why or why not? Watch a video solution.\n
Load the mtcars dataset. Fit a linear regression with\nmiles per gallon as the outcome and horsepower as the predictor. Interpret\nyour coefficients, recenter for the intercept if necessary. Watch a video solution.\n
Refer to question 4. Overlay the fit onto a scatterplot. Watch a video solution.\n
Refer to question 4. Test the hypothesis of no linear relationship between horsepower and miles per gallon. Watch a video solution.\n
Refer to question 4. Predict the miles per gallon for a horsepower of 111. Watch a video solution.\n

\n\n\n

\nResiduals

\n\n

\n\n

\nResidual variation

Residuals represent variation left unexplained by our model. We emphasize the\ndifference between residuals and errors. The errors are the unobservable true\ndeviations from the known coefficients, while residuals are the observable deviations\nfrom the estimated coefficients. In a sense, the residuals are estimates of\nthe errors.

\n\n

Consider again the diamond data set from UsingR. Recall that\nthe data is diamond prices (Singapore dollars) and diamond weight\nin carats (standard measure of diamond mass, 0.2g).\nTo get the data use library(UsingR); data(diamond). Recall that the data\nand our linear regression fit looked like the following:

\n\n\n

Recall our linear model was

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X_i + \\epsilon_i\n$$\n\n

where we are assuming that \$\\epsilon_i \\sim N(0, \\sigma^2)\$.\nOur observed outcome is \$Y_i\$ with associated\npredictor value, \$X_i\$. Let’s label the predicted outcome\nfor index \$i\$ as \$\\hat Y_i\$. Recall that we obtain\nour predictions by plugging our observed \$X_i\$ into the\nlinear regression equation:

\n\n\n$$\n\\hat Y_i = \\hat \\beta_0 + \\hat \\beta_1 X_i\n$$\n\n

The residual is defined as the difference the between the observed and predicted outcome

\n\n\n$$\ne_i = Y_i - \\hat Y_i.\n$$\n\n

The residuals are exactly the vertical distance between the observed data point\nand the associated point on the regression line. Positive residuals have associated\nY values above the fitted line and negative residuals have values below.

\n\n\n

\"Picture — Picture of the residuals for the diamond data. Residuals are the signed length\nof the red lines.

\n\n\n

Least squares minimizes the sum of the squared residuals, \$\\sum_{i=1}^n e_i^2\$.\nNote that the \$e_i\$ are observable, while the errors, \$\\epsilon_i\$ are\nnot. The residuals can be thought of as estimates of the errors.

\n\n

\nProperties of the residuals

Let’s consider some properties of the residuals.\nFirst, under our model, their expected value is 0, \$E[e_i] = 0\$.\nIf an intercept is included, \$\\sum_{i=1}^n e_i = 0\$. Note this tells\nus that the residuals are not independent. If we know \$n-1\$ of them,\nwe know the \$n^{th}\$. In fact, we will only have \$n-p\$ free\nresiduals, where \$p\$ is the number of coefficients in our regression model,\nso \$p=2\$ for linear regression with an intercept and slope.\nIf a regressor variable, \$X_i\$, is included in the model\nthen \$\\sum_{i=1}^n e_i X_i = 0\$.

\n\n

What do we use residuals for? Most importantly,\nresiduals are useful for investigating poor model fit.\nResidual plots highlight poor model fit.

\n\n

Another use for residuals is to create covariate adjusted variables. Specifically,\nresiduals can be thought of as the outcome (Y) with the\nlinear association of the predictor (X) removed. So, for example, if you wanted\nto create a weight variable with the linear effect of height removed, you would\nfit a linear regression with weight as the outcome and height as the predictor\nand take the residuals. (Note this only works if the relationship is linear.)

\n\n

Finally, we should note the different sorts of variation one encounters in regression.\nThere’s the total variability in our response, usually called total variation.\nOne then differentiates residual variation (variation after removing\nthe predictor) from systematic variation (variation explained by the regression model).\nThese two kinds of variation add up to the total variation, which we’ll see later.

\n\n

\nExample

\n\n

\n\n

The code below shows how to obtain the residuals.

\n\n

> data(diamond)\n> y <- diamond$price; x <- diamond$carat; n <- length(y)\n> fit <- lm(y ~ x)\n## The easiest way to get the residuals\n> e <- resid(fit)\n## Obtain the residuals manually, get the predicted Ys first\n> yhat <- predict(fit)\n## The residuals are y - yhat. Let's check by comparing this\n## with R's build in resid function\n> max(abs(e -(y - yhat)))\n[1] 9.486e-13\n## Let's do it again hard coding the calculation of Yhat\nmax(abs(e - (y - coef(fit)[1] - coef(fit)[2] * x)))\n[1] 9.486e-13\n

\n\n\n\n

\nResiduals versus X

\n\n

A useful plot is the residuals versus the X values. This allows us to\nzoom in on instances of poor model fit.\nWhenever we look at a residual plot, we are searching for systematic patterns\nof any sort. Here’s the plot for diamond data.

\n\n\n

Let’s go through some contrived examples to highlight. Here’s a plot of\nnonlinear data where we’ve fit a line.

\n\n\n

Here’s what happens when you focus in on the residuals.

\n\n\n

Another common feature where our model fails is when the variance around\nthe regression line is not constant. Remember our errors are assumed to be\nGaussian with a constant error. Here’s an example where heteroskedasticity\nis not at all apparent in the scatterplot.

\n\n\n

\"Scatterplot — Scatterplot demonstrating heteroskedasticity.

\n\n\n

Now look at the consequences of focusing in on the residuals.

\n\n\n

If we look at the residual plot for the diamond data, things don’t look so bad.

\n\n\n

\nEstimating residual variation

Watch this before beginning.

\n\n

We’ve talked at length about how to estimate \$\\beta_0\$ and \$\\beta_1\$.\nHowever, there’s another parameter in our model, \$\\sigma\$.\nRecall that our model is \$Y_i = \\beta_0 + \\beta_1 X_i + \\epsilon_i\$,\nwhere \$\\epsilon_i \\sim N(0, \\sigma^2)\$.

\n\n

It seems natural to use our residual variation to estimate population\nerror variation.\nIn fact, the maximum likelihood estimate of \$\\sigma^2\$ is\n\$\\frac{1}{n}\\sum_{i=1}^n e_i^2\$, the average squared residual. Since\nthe residuals have a zero mean (if an intercept is included), this is close\nto the the calculating the variance of the residuals. However, to obtain\nunbiasedness, most people use

\n\n\n$$\n\\hat \\sigma^2 = \\frac{1}{n-2}\\sum_{i=1}^n e_i^2.\n$$\n\n

The \$n-2\$ instead of \$n\$ is so that \$E[\\hat \\sigma^2] = \\sigma^2\$.\nThis is exactly analogous to dividing by \$n-1\$ in the ordinary variance\ncalculation. In fact, the ordinary variance (using var in R on a vector) is\nexactly the same as the residual variance estimate from a model that has an\nintercept and no slope. The \$n-2\$ instead of \$n-1\$ when we\ninclude a slope can be thought of as losing a degree of freedom from having\nto estimate an extra parameter (the slope).

\n\n

Most of this is typically opaque to the user, since we just grab the correct\nresidual variance output from lm. But, to solidify the concepts, let’s\ngo through the diamond example to make sure that we can hard code the estimates\non our own. (And from then on we’ll just use lm.)

\n\n

\nDiamond example

\n\n

Finding residual variance estimates.

> y <- diamond$price; x <- diamond$carat; n <- length(y)\n> fit <- lm(y ~ x)\n## the estimate from lm\n> summary(fit)$sigma\n[1] 31.84\n## directly calculating from the residuals\n> sqrt(sum(resid(fit)^2) / (n - 2))\n[1] 31.84\n

\n\n\n\n

\nSummarizing variation

\n\n

A way to think about regression is in the decomposition of variability of\nour response. The total variability in our response is the variability\naround an intercept. This is also the variance estimate from a model with\nonly an intercept:

\n\n\n$$\n\\mbox{Total variability}=\\sum_{i=1}^n (Y_i - \\bar Y)^2\n$$\n\n

The regression variability is the variability that is explained by adding the\npredictor. Mathematically, this is:

\n\n

\$\\mbox{Regression variability} = \\sum_{i=1}^n (\\hat Y_i - \\bar Y)^2\$.

\n\n

The residual variability is what’s leftover around the regression line

\n\n\n$$\n\\mbox{Residual variability} = \\sum_{i=1}^n (Y_i - \\hat Y_i)^2\n$$\n\n

It’s a nice fact that the error and regression variability add up to the\ntotal variability:

\n\n\n$$\n\\sum_{i=1}^n (Y_i - \\bar Y)^2\n= \\sum_{i=1}^n (Y_i - \\hat Y_i)^2 + \\sum_{i=1}^n (\\hat Y_i - \\bar Y)^2\n$$\n\n

Thus, we can think of regression as explaining away variability. The fact\nthat all of the quantities are positive and that they add up this way allows\nus to define the proportion of the total variability explained by the model.

\n\n

Consider our diamond example again. The plot below shows the variation\nexplained by a model with an intercept only (representing\n total variation) and then the mass is included as a linear predictor.\nNotice how much the variation decreases when including the diamond mass.

\n\n

Here’s the code:

\n\n

e = c(resid(lm(price ~ 1, data = diamond)),\n      resid(lm(price ~ carat, data = diamond)))\nfit = factor(c(rep(\"Itc\", nrow(diamond)),\n               rep(\"Itc, slope\", nrow(diamond))))\ng = ggplot(data.frame(e = e, fit = fit), aes(y = e, x = fit, fill = fit))\ng = g + geom_dotplot(binaxis = \"y\", size = 2, stackdir = \"center\", binwidth = 20)\ng = g + xlab(\"Fitting approach\")\ng = g + ylab(\"Residual price\")\ng\n

\n\n\n\n\n

\n\n\n

\nR squared

R squared is the percentage of the total variability that is explained\nby the linear relationship with the predictor

\n\n\n$$\nR^2 = \\frac{\\sum_{i=1}^n (\\hat Y_i - \\bar Y)^2}{\\sum_{i=1}^n (Y_i - \\bar Y)^2}\n$$\n\n

Here are some summary notes about R squared.

\n\n

\$R^2\$ is the percentage of variation explained by the regression model.
\n\n$$\n0 \\leq R^2 \\leq 1\n$$\n
\$R^2\$ is the sample correlation squared
\$R^2\$ can be a misleading summary of model fit.\n
- Deleting data can inflate it.
- (For later.) Adding terms to a regression model always increases \$R^2\$.
\n

\n\n

Anscombe’s residuals (named after their inventor)\nare a famous example of how R squared doesn’t tell the whole story\nabout model fit. In this example, four data sets have equivalent R squared\nvalues and beta values, but dramatically different model fits. The result is to suggest\nthat reducing data to a single number, be it R squared, a test statistic\nor a P-value, often masks important aspects of the data.\nThe code is quite simple: data(anscombe);example(anscombe).

\n\n\n

\nExercises

Fit a linear regression model to the father.son dataset with\nthe father as the predictor and the son as the outcome. Plot the\nfather’s height (horizontal axis) versus the residuals (vertical axis).\nWatch a video solution.\n
Refer to question 1. Directly estimate the residual variance and\ncompare this estimate to the output of lm. Watch a video solution.\n
Refer to question 1. Give the R squared for this model. Watch a video solution.\n
Load the mtcars dataset. Fit a linear regression with\nmiles per gallon as the outcome and horsepower as the predictor.\nPlot horsepower versus the residuals. Watch a video solution.\n
Refer to question 4. Directly estimate the residual variance and\ncompare this estimate to the output of lm. Watch a video solution.\n
Refer to question 4. Give the R squared for this model. Watch a video solution.\n

\n\n\n

\nRegression inference

\n\n

In this chapter, we’ll consider statistical inference for regression models.

\n\n

\nReminder of the model

Consider our regression model:

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X_i + \\epsilon_i\n$$\n\n

where \$\\epsilon \\sim N(0, \\sigma^2)\$. Let’s consider some\nways for doing inference for our regression parameters.\nFor this development, we assume that the true model is known.\nWe also assume that you’ve seen confidence intervals and hypothesis tests\nbefore. If not, consider taking the Statistical Inference course and reading the accompanying book\nbefore approaching this material.

\n\n

Remember our estimates:

\n\n\n$$\n\\hat \\beta_0 = \\bar Y - \\hat \\beta_1 \\bar X\n$$\n\n

and

\n\n\n$$\n\\hat \\beta_1 = Cor(Y, X) \\frac{Sd(Y)}{Sd(X)}.\n$$\n\n

\nReview

Let’s review some important components of statistical inference.\nConsider statistics like the following:

\n\n\n$$\n\\frac{\\hat \\theta - \\theta}{\\hat \\sigma_{\\hat \\theta}}\n$$\n\n

where \$\\hat \\theta \$ is an estimate of interest, \$\\theta\$\nis the estimand of interest and \$ \\hat \\sigma_{\\hat \\theta} \$ is\nthe standard error of \$\\hat \\theta\$. We see that in many cases\nsuch statistics often have the following properties:

\n\n

They are normally distributed and have a finite sample\n Student’s T distribution under normality assumptions.
They can be used to test \$H_0 : \\theta = \\theta_0\$ versus\n \\(H_a : \\theta >, \n
They can be used to create a confidence interval for \$\\theta\$ via\n\$\\hat \\theta \\pm Q_{1-\\alpha/2} \\hat \\sigma_{\\hat \\theta}\$\n where \$Q_{1-\\alpha/2}\$ is the relevant quantile from either a\n normal or T distribution.

\n\n

In the case of regression with iid Gaussian sampling assumptions on the errors,\nour inferences will follow very similarly to what you saw in your inference class.

\n\n

We won’t cover asymptotics for regression analysis, but suffice it to say that under assumptions\non the ways in which the \$X\$ values are collected, the iid sampling model and mean model,\nthe normal results hold to create intervals and confidence intervals

\n\n

\nResults for the regression parameters

\n\n

First, we need standard errors for our regression parameters. These are given by:

\n\n\n$$\n\\sigma_{\\hat \\beta_1}^2 = Var(\\hat \\beta_1) = \\sigma^2 / \\sum_{i=1}^n (X_i - \\bar X)^2\n$$\n\n

and

\n\n\n$$\n\\sigma_{\\hat \\beta_0}^2 = Var(\\hat \\beta_0) = \\left(\\frac{1}{n} + \\frac{\\bar X^2}{\\sum_{i=1}^n (X_i - \\bar X)^2 }\\right)\\sigma^2\n$$\n\n

In practice, \$\\sigma^2\$ is replaced by its residual variance estimate covered\nin the last chapter.

\n\n

Given how often this came up in inference,\nit’s probably not surprising that under iid Gaussian errors

\n\n\n$$\n\\frac{\\hat \\beta_j - \\beta_j}{\\hat \\sigma_{\\hat \\beta_j}}\n$$\n\n

follows a t distribution with n-2 degrees of freedom and a normal distribution for large n.\nThis can be used to create confidence intervals and perform hypothesis tests.

\n\n

\nExample diamond data set

\n\n

\n\n

Let’s go through a didactic example using our diamond pricing data.\nFirst, let’s define our outcome, predictor and estimate all of the parameters.\n(Note, again we’re hard coding these results, but lm will give it to us automatically).

\n\n

library(UsingR); data(diamond)\ny <- diamond$price; x <- diamond$carat; n <- length(y)\nbeta1 <- cor(y, x) * sd(y) / sd(x)\nbeta0 <- mean(y) - beta1 * mean(x)\ne <- y - beta0 - beta1 * x\nsigma <- sqrt(sum(e^2) / (n-2))\nssx <- sum((x - mean(x))^2)\n

\n\n\n\n

Now let’s calculate the standard errors for our regression coefficients\nand the t statistic. The natural null\nhypotheses are \$H_0: \\beta_j = 0\$. So our t statistics are\njust the estimates divided by their standard errors.

\n\n

seBeta0 <- (1 / n + mean(x) ^ 2 / ssx) ^ .5 * sigma\nseBeta1 <- sigma / sqrt(ssx)\ntBeta0 <- beta0 / seBeta0\ntBeta1 <- beta1 / seBeta1\n

\n\n\n\n

Now let’s consider getting P-values. Recall that P-values are the probability\nof getting a statistic as or larger than was actually obtained, where the\nprobability is calculated under the null hypothesis. Below I also do some\nformatting to get it to look like the output from lm.

\n\n

> pBeta0 <- 2 * pt(abs(tBeta0), df = n - 2, lower.tail = FALSE)\n> pBeta1 <- 2 * pt(abs(tBeta1), df = n - 2, lower.tail = FALSE)\n> coefTable <- rbind(c(beta0, seBeta0, tBeta0, pBeta0), c(beta1, seBeta1, tBeta1, pB\\\neta1))\n> colnames(coefTable) <- c(\"Estimate\", \"Std. Error\", \"t value\", \"P(>|t|)\")\n> rownames(coefTable) <- c(\"(Intercept)\", \"x\")\n> coefTable\n            Estimate Std. Error t value   P(>|t|)\n(Intercept)   -259.6      17.32  -14.99 2.523e-19\nx             3721.0      81.79   45.50 6.751e-40\n

\n\n\n\n

The first column is the actual estimates. The second is the standard\nerrors, the third is the t value (the first divided by the second) and the\nfinal is the t probability of getting an unsigned statistic that large under\nthe null hypothesis (the P-value for the two sided test).\nOf course, we don’t actually go through this exercise every time; lm does\nit for us.

\n\n

> fit <- lm(y ~ x);\n> summary(fit)$coefficients\n            Estimate Std. Error t value  Pr(>|t|)\n(Intercept)   -259.6      17.32  -14.99 2.523e-19\nx             3721.0      81.79   45.50 6.751e-40\n

\n\n\n\n

Remember, we reject if our P-value is less than our desired type I error\nrate. In both cases the test is for whether or not the parameter is zero.\nThis is almost always of interest for the slope, but frequently a zero\nintercept isn’t of interest so that P-value is often disregarded.

\n\n

For the slope, a value of zero represents no linear relationship between\nthe predictor and response. So, the P-value is for performing a test of whether\nany (linear) relationship exist or not.

\n\n

\nGetting a confidence interval

\n\n

Recall from your inference class, a fair number of confidence intervals\ntake the form of an estimate plus or minus a t quantile times a standard\nerror. Let’s use that formula to create confidence intervals for our\nregression parameters. Let’s first do the intercept.

\n\n

> sumCoef <- summary(fit)$coefficients\n> sumCoef[1,1] + c(-1, 1) * qt(.975, df = fit$df) * sumCoef[1, 2]\n[1] -294.5 -224.8\n

\n\n\n\n

Now let’s do the slope:

\n\n

> (sumCoef[2,1] + c(-1, 1) * qt(.975, df = fit$df) * sumCoef[2, 2]) / 10\n[1] 355.6 388.6\n

\n\n\n\n

So, we would interpret this as: “with 95% confidence, we estimate that a 0.1 carat increase in\ndiamond size results in a 355.6 to 388.6 increase in price in (Singapore) dollars”.

\n\n

\nPrediction of outcomes

\n\n

Watch this before beginning.

\n\n

Finally, let’s consider prediction again.\nConsider the problem of predicting Y at a value of X. In our example, this is\npredicting the price of a diamond given the carat.

\n\n

We’ve already covered that the estimate for prediction at point \$x_0\$ is:

\n\n\n$$\n\\hat \\beta_0 + \\hat \\beta_1 x_0\n$$\n\n

A standard error is needed to create a prediction interval. This is important,\nsince predictions by themselves don’t convey anything about how accurate we\nwould expect the prediction to be. Take our diamond example. Because the model\nfits so well, we would be surprised if we tried to sell a diamond and the offers\nwere well off our model prediction (since it seems to fit quite well).

\n\n

There’s a subtle, but important, distinction between intervals for the regression\nline at point \$x_0\$ and the prediction of what a \$y\$ would be\nat point \$x_0\$. What differs is the standard error:

\n\n

For the line at \$x_0\$ the standard error is,

\n\n\n$$\n\\hat \\sigma\\sqrt{\\frac{1}{n} + \\frac{(x_0 - \\bar X)^2}{\\sum_{i=1}^n (X_i - \\bar X)^2}}\n$$\n\n

For the prediction interval at \$x_0\$ the standard error is

\n\n\n$$\n\\hat \\sigma\\sqrt{1 + \\frac{1}{n} + \\frac{(x_0 - \\bar X)^2}{\\sum_{i=1}^n (X_i - \\bar X)^2}}\n$$\n\n

Notice that the prediction interval standard error is a little larger than the error for a line.\nThink of it this way. If we want to predict a Y value at a particular X value,\nand we knew the actual true slope and intercept, there would still be error. However,\nif we only wanted to predict the value at the line at that X value, there would be\nno variance, since we already know the line.

\n\n

Thus, the variation for the line only considers how hard it is to estimate the\nregression line at that X value. The prediction interval includes that variation,\nas well as the extra variation unexplained by the relationship between Y and X.\nSo, it has to be a little wider.

\n\n

For the diamond example, here’s both the mean value and prediction interval.\n(code and plot). Notice that to get the various intervals, one has to\nuse one of the options interval=\"confidence\" or interval=\"prediction\" in the\nprediction function.

\n\n

library(ggplot2)\nnewx = data.frame(x = seq(min(x), max(x), length = 100))\np1 = data.frame(predict(fit, newdata= newx,interval = (\"confidence\")))\np2 = data.frame(predict(fit, newdata = newx,interval = (\"prediction\")))\np1$interval = \"confidence\"\np2$interval = \"prediction\"\np1$x = newx$x\np2$x = newx$x\ndat = rbind(p1, p2)\nnames(dat)[1] = \"y\"\ng = ggplot(dat, aes(x = x, y = y))\ng = g + geom_ribbon(aes(ymin = lwr, ymax = upr, fill = interval), alpha = 0.2)\ng = g + geom_line()\ng = g + geom_point(data = data.frame(x = x, y=y), aes(x = x, y = y), size = 4)\ng\n

\n\n\n\n\n

\n\n\n

\nSummary notes

Both intervals have varying widths.\n
- Least width at the mean of the Xs.
\n
We are quite confident in the regression line, so that\ninterval is very narrow.\n
- If we knew \$\\beta_0\$ and \$\\beta_1\$ this interval would have zero width.
\n
The prediction interval must incorporate the variability\nin the data around the line.\n
- Even if we knew \$\\beta_0\$ and \$\\beta_1\$ this interval would still have width.\n*
\n

\n\n

\nExercises

Test whether the slope coefficient for the father.son data is different from zero (father as predictor, son as outcome). Watch a video solution.\n
Refer to question 1. Form a confidence interval for the slope coefficient. Watch a video solution\n
Refer to question 1. Form a confidence interval for the intercept (center the fathers’ heights first to get an intercept that is easier to interpret). Watch a video solution.\n
Refer to question 1. Form a mean value interval for the expected son’s height at the average father’s height. Watch a video solution.\n
Refer to question 1. Form a prediction interval for the son’s height at the average father’s height. Watch a video solution.\n
Load the mtcars dataset. Fit a linear regression with miles per gallon as the outcome and horsepower as the predictor. Test whether or not the horsepower power coefficient is statistically different from zero. Interpret your test.
Refer to question 6. Form a confidence interval for the slope coefficient.
Refer to question 6. Form a confidence interval for the intercept (center the HP variable first).
Refer to question 6. Form a mean value interval for the expected MPG for the average HP.
Refer to question 6. Form a prediction interval for the expected MPG for the average HP.
Refer to question 6. Create a plot that has the fitted regression line plus curves at the expected value and prediction intervals.

\n\n\n

\nMultivariable regression analysis

\n\n

In this chapter we extend linear regression so that our models\ncan contain more variables. A natural first approach is to assume\nadditive effects, basically extending our linear model to a plane\nor hyperplane. This technique represents one of the most widely\nused and successful methods in statistics.

\n\n

\nMultivariable regression analyses: adjustment

If I were to present evidence of a relationship between\nbreath mint usage (mints per day, X) and pulmonary function\n(measured in FEV), you would be skeptical.\nLikely, you would say, ‘smokers tend to use more breath mints than non smokers, smoking is related to a loss in pulmonary function. That’s probably the culprit.’\nIf asked what would convince you, you would likely say, ‘If non-smoking breath mint users had lower lung function than non-smoking non-breath mint users and, similarly, if smoking breath mint users had lower lung function than smoking non-breath mint users, I’d be more inclined to believe you’.\nIn other words, to even consider my results, I would have to demonstrate that they hold while holding smoking status fixed.

\n\n

This is one of the main uses of multivariate regression, to consider a relationship\nbetween a predictor and response, while accounting for other variables.

\n\n

\nMultivariable regression analyses: prediction

An insurance company is interested in how last year’s claims can predict a person’s time in the hospital this year. They want to use an enormous amount of data contained in claims to predict a single number. Simple linear regression is not equipped to handle more than one predictor.\nHow can one generalize SLR to incorporate lots of regressors for the purpose of prediction?\nWhat are the consequences of adding lots of regressors? Surely there must be consequences to throwing variables in that aren’t related to Y? Surely there must also be consequences to omitting variables that are?

\n\n

\nThe linear model

The general linear model extends simple linear regression (SLR)\nby adding terms linearly into the model.

\n\n\n$$\nY_i = \\beta_1 X_{1i} + \\beta_2 X_{2i} + \\ldots +\n\\beta_{p} X_{pi} + \\epsilon_{i}\n= \\sum_{k=1}^p X_{ki} \\beta_k + \\epsilon_{i}\n$$\n\n

Here \$X_{1i}=1\$ typically, so that an intercept is included.\nLeast squares (and hence maximum likelihood estimates under iid Gaussianity\nof the errors) minimizes:

\n\n\n$$\n\\sum_{i=1}^n \\left(Y_i - \\sum_{k=1}^p X_{ki} \\beta_k\\right)^2\n$$\n\n

Note, the linearity referred to in these models is linearity in the coefficients.\nThus

\n\n\n$$\nY_i = \\beta_1 X_{1i}^2 + \\beta_2 X_{2i}^2 + \\ldots +\n\\beta_{p} X_{pi}^2 + \\epsilon_{i}\n$$\n\n

is still a linear model. We’ve just squared the elements of the\npredictor variables.

\n\n

\nEstimation

\n\n

Watch this video demonstration.

\n\n

Recall, the LS estimate for regression through the origin is,

\n\n

\$E[Y_i]=X_{1i}\\beta_1\$, was \$\\sum X_i Y_i / \\sum X_i^2.\$

\n\n

Let’s consider two regressors, \$E[Y_i] = X_{1i}\\beta_1 + X_{2i}\\beta_2 = \\mu_i\$.\nLeast squares tries to minimize:

\n\n\n$$\n\\sum_{i=1}^n (Y_i - X_{1i} \\beta_1 - X_{2i} \\beta_2)^2\n$$\n\n

We describe fitting with two regressors using residuals, since it will help\nus to understand how multivariable regression adjusts an effect for another\nvariable. The result is that the estimate for \$\\beta_1\$ is:

\n\n\n$$\n\\hat \\beta_1 = \\frac{\\sum_{i=1}^n e_{i, Y | X_2} e_{i, X_1 | X_2}}{\\sum_{i=1}^n e_{i, X_1 | X_2}^2},\n$$\n\n

where, \$e_{i, Y|X_2}\$ is the residual having fit \$Y\$ on \$X_2\$\nand \$e_{i, X_1|X_2}\$ is the residual having fit \$X_1\$ on \$X_2\$.\nThat is, the regression estimate for \$\\beta_1\$ is the regression\nthrough the origin estimate having regressed \$X_2\$ out of both\nthe response and the predictor. Similarly, the regression estimate for\n\$\\beta_2\$ is the regression through the origin estimate having regressed\n\$X_1\$ out of both the response and the predictor.

\n\n

More generally, multivariate regression estimates are exactly those having removed the linear relationship of the other variables from both the regressor and response.\nThis demonstrates the sense in which multivariate regression variables adjust\nfor the effect of the other variables.

\n\n

\nExample with two variables, simple linear regression

\n\n

We already have one of the most important examples of two variables down, linear\nregression. The linear regression model is:

\n\n\n$$\nY_{i} = \\beta_1 X_{1i} + \\beta_2 X_{2i}\n$$\n\n

where \$X_{2i} = 1\$ is an intercept term. Let’s double check our\nrule, since we already know what the least squares estimates are in this case.

\n\n

Notice the fitted coefficient of \$X_{2i}\$ on \$Y_{i}\$ is\n\$\\bar Y\$, the mean of the Ys. Then the residuals are\n\$e_{i, Y | X_2} = Y_i - \\bar Y\$.

\n\n

Similarly, the fitted coefficient of \$X_{2i}\$ on \$X_{1i}\$\nis \$\\bar X_1\$. Then, the residuals are\n\$e_{i, X_1 | X_2}= X_{1i} - \\bar X_1\$.

\n\n

Thus let’s work out the estimate for \$\\beta_1\$:

\n\n\n$$\n\\hat \\beta_1 = \\frac{\\sum_{i=1}^n e_{i, Y | X_2} e_{i, X_1 | X_2}}{\\sum_{i=1}^n e_{i, X_1 | X_2}^2} = \\frac{\\sum_{i=1}^n (X_i - \\bar X)(Y_i - \\bar Y)}{\\sum_{i=1}^n (X_i - \\bar X)^2}\n= Cor(X, Y) \\frac{Sd(Y)}{Sd(X)}\n$$\n\n

This agrees with the estimate that we came up with before.

\n\n

\nThe general case

In the general case with \$p\$ regressors,\nleast squares solutions have to minimize:

\n\n\n$$\n\\sum_{i=1}^n (Y_i - X_{1i}\\beta_1 - \\ldots - X_{pi}\\beta_p)^2\n$$\n\n

The least squares estimate for the coefficient of a multivariate regression model is exactly regression through the origin with the linear relationships with the other regressors removed from both the regressor and outcome by taking residuals.\nIn this sense, multivariate regression “adjusts” a coefficient for the linear impact of the other variables.

\n\n

Of course, we don’t fit multivariate regression models in this way\never in practice, we rely on software like lm. In fact, the\nprograms that fit multivariate regression don’t do it this way either. However,\nthinking in the terms of residuals is the most conceptually useful way to think about what multivariate regression is accomplishing.

\n\n

\nSimulation demonstrations

\n\n

\n\n

Let’s do some simulation exercises to convince ourselves how multivariable\nregression works.

\n\n

\nLinear model with three variables

\n\n

> n = 100; x = rnorm(n); x2 = rnorm(n); x3 = rnorm(n)\n## Generate the data\n> y = 1 + x + x2 + x3 + rnorm(n, sd = .1)\n## Get the residuals having removed X2 and X3 from X1 and Y\n> ey = resid(lm(y ~ x2 + x3))\n> ex = resid(lm(x ~ x2 + x3))\n## Fit regression through the origin with the residuals\n> sum(ey * ex) / sum(ex ^ 2)\n[1] 1.009\n## Double check with lm\n> coef(lm(ey ~ ex - 1))\n    ex\n1.009\n## Fit the full linear model to show that it agrees\ncoef(lm(y ~ x + x2 + x3))\n(Intercept)           x          x2          x3\n     1.0202      1.0090      0.9787      1.0064\n

\n\n\n

You can see that the estimate, 1.009, is obtained by the recipe we outlined;\nand it agrees with the function lm produces for us.

\n\n

\nInterpretation of the coefficients

\n\n

Let’s go through the interpretation of the coefficients. Consider the\npredicted mean for a given set of values of the regressors:

\n\n\n$$\nE[Y | X_1 = x_1, \\ldots, X_p = x_p] = \\sum_{k=1}^p x_{k} \\beta_k.\n$$\n\n

Now consider incrementing \$X_1\$ (and only \$X_1\$) by 1.

\n\n\n$$\nE[Y | X_1 = x_1 + 1, \\ldots, X_p = x_p] = (x_1 + 1) \\beta_1 + \\sum_{k=2}^p x_{k} \\beta_k\n$$\n\n

Now let’s subtract these two equations:

\n\n\n$$\nE[Y | X_1 = x_1 + 1, \\ldots, X_p = x_p] - E[Y | X_1 = x_1, \\ldots, X_p = x_p]\n= (x_1 + 1) \\beta_1 + \\sum_{k=2}^p x_{k} \\beta_k - \\sum_{k=1}^p x_{k} \\beta_k = \\beta_1\n$$\n\n

Thus, the interpretation of a multivariate regression coefficient is the expected change in the response per unit change in the regressor, holding all of the other regressors fixed.\nThe latter part of the phrase is important, by holding the other regressors\nconstant, we are investigating an adjusted effect, just like we described in the\nsmoking and breath mint usage example at the beginning of the chapter.

\n\n

In the next chapter, we’ll do examples and go over context-specific\ninterpretations.

\n\n

\nFitted values, residuals and residual variation

All of our simple linear regression quantities can be extended to linear models.\nHere we list them out in one place. Our statistical model is:

\n\n\n$$\nY_i = \\sum_{k=1}^p X_{ki} \\beta_{k} + \\epsilon_{i}\n$$\n\n

where \$\\epsilon_i \\sim N(0, \\sigma^2)\$.\nOur fitted responses are:

\n\n\n$$\n\\hat Y_i = \\sum_{k=1}^p X_{ki} \\hat \\beta_{k}.\n$$\n\n

We can define our residuals exactly as in linear regression:

\n\n\n$$\ne_i = Y_i - \\hat Y_i\n$$\n\n

Our variance estimate is.

\n\n\n$$\n\\hat \\sigma^2 = \\frac{1}{n-p} \\sum_{i=1}^n e_i ^2\n$$\n\n

To get predicted responses at new values, \$x_1, \\ldots, x_p\$,\nsimply plug them into the linear model \$\\sum_{k=1}^p x_{k} \\hat \\beta_{k}\$.

\n\n

Coefficients have standard errors, we can label them as\n\$\\hat \\sigma_{\\hat \\beta_k}\$, and

\n\n\n$$\n\\frac{\\hat \\beta_k - \\beta_k}{\\hat \\sigma_{\\hat \\beta_k}}\n$$\n\n

follows a t distribution with \$n-p\$ degrees of freedom.\nPredicted responses have standard errors and we can calculate predicted\nand expected response intervals.

\n\n

\nSummary notes on linear models

Linear models are the single most important applied statistical and machine\nlearning technique, by far.
Some amazing things that you can accomplish with linear models\n
- Decompose a signal into its harmonics.
- Flexibly fit complicated functions.
- Fit factor variables as predictors.
- Uncover complex multivariate relationships with the response.
- Build accurate prediction models.
\n

\n\n

\nExercises

\n\n

Load the dataset Seatbelts as part of the datasets package via data(Seatbelts). Use\nas.data.frame to convert the object to a dataframe. Fit a linear model of driver deaths\nwith kms and PetrolPrice as predictors. Interpret your results. Watch a video Solution\n
Predict the number of driver deaths at the average kms and petrol levels. Watch a video solution.\n
Take the residual for DriversKilled having regressed out kms and an intercept Watch a video solution.\nand the residual for PetrolPrice having regressed out kms and an intercept. Fit a regression\nthrough the origin of the two residuals and show that it is the same as your coefficient\nobtained in question 1.
Take the residual for DriversKilled having regressed out PetrolPrice and an intercept.\nTake the residual for kms having regressed out PetrolPrice and an intercept. Fit a regression\nthrough the origin of the two residuals and show that it is the same as your coefficient\nobtained in question 1.

\n\n\n

\nMultivariable examples and tricks

\n\n

Watch this before beginning.

\n\n

In this chapter we cover a few examples of multivariate regression in order\nto get a hands on sense of the basics.

\n\n

\nData set for discussion

\n\n

We’ll start with the Swiss dataset that is part of the\ndatasets package. This can be loaded in R with:

\n\n

> require(datasets); data(swiss); ?swiss\n

\n\n\n\n

Standardized fertility measure and socio-economic indicators for each of 47\nFrench-speaking provinces of Switzerland at about 1888.

\n\n

A data frame with 47 observations on 6 variables, each of which is in percent,\ni.e., in [0, 100].

\n\n

[,1]   Fertility          a common standardized fertility measure\n[,2]   Agriculture        percent of males involved in agriculture as occupation\n[,3]\t Examination        percent draftees receiving highest on army examination\n[,4]\t Education          percent education beyond primary school for draftees\n[,5]\t Catholic           percent catholic (as opposed to protestant)\n[,6]\t Infant.Mortality   live births who live less than 1 year\n

\n\n\n\n

All variables but Fertility give percentages of the population.

\n\n\n

Let’s see the result of calling lm on this data set.

\n\n

> summary(lm(Fertility ~ . , data = swiss))\n                 Estimate Std. Error t value  Pr(>|t|)\n(Intercept)       66.9152   10.70604   6.250 1.906e-07\nAgriculture       -0.1721    0.07030  -2.448 1.873e-02\nExamination       -0.2580    0.25388  -1.016 3.155e-01\nEducation         -0.8709    0.18303  -4.758 2.431e-05\nCatholic           0.1041    0.03526   2.953 5.190e-03\nInfant.Mortality   1.0770    0.38172   2.822 7.336e-03\n

\n\n\n\n

Agriculture is expressed in percentages (0 - 100), representing the percentage\nof the male population involved in agriculture.
\nThe regression slope estimate for this variable is -0.1721. We interpret\nthis coefficient as follows:

\n\n

Our model estimates an expected 0.17 decrease in standardized\n fertility for every 1% increase in percentage of males involved in\n agriculture, holding the remaining variables constant.

\n\n

Note that the the t-test for \$H_0: \\beta_{Agri} = 0\$ versus\n\$H_a: \\beta_{Agri} \\neq 0\$ is significant since 0.0187 is less\nthat typical benchmarks (0.05, for example). Note that, by default, R is\nreporting the P-value for the two sided test. If you want the one sided test,\ncalculate it directly using the T-statistic and the degrees of freedom.\n(You can figure it out from the two sided P-value, but it’s easy to get tripped\nup with signs.)

\n\n

Interestingly, the unadjusted estimate is

\n\n

summary(lm(Fertility ~ Agriculture, data = swiss))$coefficients\n            Estimate Std. Error t value  Pr(>|t|)\n(Intercept)  60.3044    4.25126  14.185 3.216e-18\nAgriculture   0.1942    0.07671   2.532 1.492e-02\n

\n\n\n\n

Notice that the sign of the slope estimate reversed! This is an\nexample of so-called “Simpson’s Paradox”. This purported paradox (which is actually not\na paradox at all) simply points out that unadjusted and adjusted effects can\nbe the reverse of each other. Or in other words, the apparent relationship between X and Y may\nchange if we account for Z. Let’s explore multivariate adjustment and\nsign reversals with simulation.

\n\n

\nSimulation study

\n\n

Below we simulate 100 random variables with a linear relationship between X1, X2 and Y.\nNotably, we generate X1 as a linear function of X2. In this simulation, X1 has a negative\nadjusted effect on Y while X2 has a positive adjusted effect (adjusted referring to the effect\nincluding both variables). However, X1 is related to X2. Notice our unadjusted effect of\nX1 is of the opposite sign (and way off), while the adjusted one is about right. What’s happening?\nOur unadjusted model is picking up the effect X2 as it’s represented in X1. Play around\nwith the generating coefficients to see how you can make the estimated relationships very\ndifferent than the generating ones. More than anything, this illustrates that multivariate\nmodeling is hard stuff.

\n\n

> n = 100; x2 <- 1 : n; x1 = .01 * x2 + runif(n, -.1, .1); y = -x1 + x2 + rnorm(n, s\\\nd = .01)\n> summary(lm(y ~ x1))$coef\n\n            Estimate Std. Error t value  Pr(>|t|)\n(Intercept)    1.454      1.079   1.348 1.807e-01\nx1            96.793      1.862  51.985 3.707e-73\n>summary(lm(y ~ x1 + x2))$coef\n             Estimate Std. Error  t value   Pr(>|t|)\n(Intercept)  0.001933  0.0017709    1.092  2.777e-01\nx1          -1.020506  0.0163560  -62.393  4.211e-80\nx2           1.000133  0.0001643 6085.554 1.544e-272\n

\n\n\n\n

To confirm what’s going on, let’s look at some plots. In the left\npanel, we plot Y versus X1. Notice the positive relationship. However,\nif we look at X2 (the color) notice that it clearly goes up with Y\nas well. If we adjust both the X1 and Y variable by taking the\nresidual after having regressed X2, we get the actual correct\nrelationship between X1 and Y.

\n\n\n

\nBack to this data set

In our data set, the sign reverses itself with the inclusion of Examination and Education. However, the percent of males in the province working in agriculture is negatively related to educational attainment (correlation of -0.6395) and Education and Examination (correlation of 0.6984) are obviously measuring similar things.

\n\n

So, given now that we know including correlated variables with our variable\nof interest into our regression relationship can drastically change things\nwe have to ask: “Is the positive marginal an artifact for not having accounted for, say, Education level? (Education does have a stronger effect, by the way.)” At the minimum, anyone claiming that provinces that are more agricultural have higher fertility rates would immediately be open to criticism that the real effect is due to Education levels.

\n\n

You might think then, why don’t I just always include all variables\nthat I have into my regression model to avoid incorrectly adjusted effects?\nOf course, it’s not this easy and there’s negative consequences to\nincluding variables unnecessarily. We’ll talk more about model building\nand the process of choosing which variables to include or not in the\nnext chapter.

\n\n

\nWhat if we include a completely unnecessary variable?

Next chapter we’ll discuss working with a collection of\ncorrelated predictors. But you might wonder, what happens\nif you include a predictor that’s completely unnecessary.\nLet’s try some computer experiments with our fertility data.\nIn the code below,\nz adds no new linear information, since it’s a linear\ncombination of variables already included. R just drops\nterms that are linear combination of other terms.

\n\n

> z <- swiss$Agriculture + swiss$Education\nlm(Fertility ~ . + z, data = swiss)\\\n> Call:\nlm(formula = Fertility ~ . + z, data = swiss)\n\nCoefficients:\n     (Intercept)       Agriculture       Examination         Education          Cath\\\nolic  \n          66.915            -0.172            -0.258            -0.871             0\\\n.104  \nInfant.Mortality                 z  \n           1.077                NA  \n

\n\n\n\n

This is a fundamental point of multivariate regression: regression\nmodels fit the linear space of the regressors. Therefore, any linear\nreorganization of the regressors will result in an equivalent fit,\nwith different covariates of course. However, the percentage of\nthe variance explained in the response will be constant. It is only\nthrough adding variables that are not perfectly explained by the\nexisting ones that one can explain more variation in the response.\nSo, for example, models with covariates i) X and Z, ii) X+Z and X-Z and\niii) 2X and 4Z will\nall explain the same amount of variation in Y. A third variable, W say,\nwill only explain more variation in Y if it’s not perfectly explained\nby X and Z. R lets you know when you’ve done this by putting redundant\nvariables as having NA coefficients.

\n\n

\nDummy variables are smart

\n\n

\n\n

It is interesting to note that models with\nfactor variables as predictors are simply special cases of regression\nmodels. As an example, consider the linear model:

\n\n\n$$\nY_i = \\beta_0 + X_{i1} \\beta_1 + \\epsilon_{i}\n$$\n\n

where each \$X_{i1}\$ is binary so that it is a 1 if measurement\n\$i\$ is in a group and 0 otherwise. As an example, consider\na variable as treated versus not in a clinical trial. Or, in a more\ndata science context, consider an A/B test comparing two ad campaigns\nwhere Y is the click through rate.

\n\n

Refer back to our model. For people in the group, {$$}E[Y_i] = \\beta_0\n+ \\beta_1{/$$}, and for people not in the group, {$$}E[Y_i] =\n\\beta_0{/$$}. The least squares fits work out to be {$$}\\hat \\beta_0\n+ \\hat \\beta_1{/$$} as the mean for those in the group and {$$}\\hat\n\\beta_0{/$$} as the mean for those not in the group. The variable\n\$\\beta_1\$ is interpreted as the increase or decrease in the\nmean comparing those in the group to those not. The T-test for that\ncoefficient is exactly the two group T test with a common variance.

\n\n

Finally, note including a binary variable that is 1 for those not in the group would be redundant,\nit would create three parameters to describe two means. Moreover, we know from the last section that\nincluding redundant variables will result in R just setting one of them to NA. We know that the intercept\ncolumn is a column of ones, the group variable is one for those in the group while a variable for\nthose not in the group would just be the subtraction of these two. Thus, it’s linearly redundant and unnecessary.

\n\n

\nMore than two levels

Consider a multilevel factor level. For didactic reasons, let’s say a three level factor. As an example\nconsider a variable for US political party affiliation: Republican, Democrat, Independent/other. Let’s use the model:

\n\n\n$$\nY_i = \\beta_0 + X_{i1} \\beta_1 + X_{i2} \\beta_2 + \\epsilon_i.\n$$\n\n

Here the variable \$X_{i1}\$ is 1 for Republicans and 0 otherwise,\nthe variable \$X_{i2}\$ is 1 for Democrats and 0 otherwise. As before, we don’t\nneed an \$X_{i3}\$ for Independent/Other, since it would be redundant.

\n\n

So now consider the implications of more model. If person\n\$i\$ is Republican then \$E[Y_i] = \\beta_0 +\\beta_1\$. On the other hand,\nIf person \$i\$ is Democrat then \$E[Y_i] = \\beta_0 + \\beta_2\$.\nFinally, if \$i\$ is Independent/Other \$E[Y_i] = \\beta_0\$.

\n\n

So, we can interpret our coefficients as follows.\n\$\\beta_1\$ compares the mean for Republicans to that of Independents/Others.\n\$\\beta_2\$ compares the mean for Democrats to that of Independents/Others.\n\$\\beta_1 - \\beta_2\$ compares the mean for Republicans to that of Democrats.\nNotice the coefficients are all comparisons to the category that we left out, Independents/Others.\nIf one category is an obvious reference category, chose that one to leave out.\nIn R, if our variable is a factor variable, it will create the dummy variables for us and pick one of the\nlevels to be the reference level. Let’s go through an example to see.

\n\n

\nInsect Sprays

\n\n

Let’s consider a model with factors. Consider the InsectSprays dataset in R. The data\nmodels the number of dead insects from different pesticides. Since it’s not clear from the documentation,\nlet’s assume (probably accurately)\nthat these were annoying bad insects, like fleas, mosquitoes or cockroaches, and not good ones like butterflies\nor ladybugs. After getting over that mental hurdle, let’s plot the data.

\n\n

require(datasets);data(InsectSprays); require(stats); require(ggplot2)\ng = ggplot(data = InsectSprays, aes(y = count, x = spray, fill  = spray))\ng = g + geom_violin(colour = \"black\", size = 2)\ng = g + xlab(\"Type of spray\") + ylab(\"Insect count\")\ng\n

\n\n\n\n

Here’s the plot. There are probably better ways to model this data, but let’s use a\nlinear model just to illustrate factor variables.

\n\n\n

\n\n\n

First, let’s set Spray A as the reference (the default, since it has the lowest alphanumeric factor level).

\n\n

> summary(lm(count ~ spray, data = InsectSprays))$coef\n            Estimate Std. Error t value  Pr(>|t|)\n(Intercept)  14.5000      1.132 12.8074 1.471e-19\nsprayB        0.8333      1.601  0.5205 6.045e-01\nsprayC      -12.4167      1.601 -7.7550 7.267e-11\nsprayD       -9.5833      1.601 -5.9854 9.817e-08\nsprayE      -11.0000      1.601 -6.8702 2.754e-09\nsprayF        2.1667      1.601  1.3532 1.806e-01\n

\n\n\n\n

Therefore, 0.8333 is the estimated mean comparing Spray B to Spray A (as B - A),\n-12.4167 compares Spray C to Spray A (as C - A) and so on. The inferential statistics: standard\nerrors, t value and P-value all correspond to those comparisons. The intercept, 14.5, is the\nmean for Spray A. So, its inferential statistics are testing whether or not the mean for\nSpray A is zero. As is often the case, this test isn’t terribly informative and often yields\nextremely small statistics (since we know the spray kills some bugs). The estimated mean for Spray B\nis its effect plus the intercept (14.5 + 0.8333); the estimated mean for Spray C is 14.5 - 12.4167 (its\n effect plus the intercept) and so on for the rest of the factor levels.

\n\n

Let’s hard code the factor levels so we can directly see what’s going on. Remember, we simply\nleave out the dummy variable for the reference level.

\n\n

> summary(lm(count ~\n             I(1 * (spray == 'B')) + I(1 * (spray == 'C')) +\n             I(1 * (spray == 'D')) + I(1 * (spray == 'E')) +\n             I(1 * (spray == 'F'))\n           , data = InsectSprays))$coef\n                      Estimate Std. Error t value  Pr(>|t|)\n(Intercept)            14.5000      1.132 12.8074 1.471e-19\nI(1 * (spray == \"B\"))   0.8333      1.601  0.5205 6.045e-01\nI(1 * (spray == \"C\")) -12.4167      1.601 -7.7550 7.267e-11\nI(1 * (spray == \"D\"))  -9.5833      1.601 -5.9854 9.817e-08\nI(1 * (spray == \"E\")) -11.0000      1.601 -6.8702 2.754e-09\nI(1 * (spray == \"F\"))   2.1667      1.601  1.3532 1.806e-01\n

\n\n\n\n

Of course, it’s identical. You might further ask yourself, what would\nhappen if I included a dummy variable for Spray A? Would the world implode?\nNo, it just realizes that one of the dummy variables is redundant and drops\nit.

\n\n

> summary(lm(count ~\n   I(1 * (spray == 'B')) + I(1 * (spray == 'C')) +  \n   I(1 * (spray == 'D')) + I(1 * (spray == 'E')) +\n   I(1 * (spray == 'F')) + I(1 * (spray == 'A')), data = InsectSprays))$coef\n\n                      Estimate Std. Error t value  Pr(>|t|)\n(Intercept)            14.5000      1.132 12.8074 1.471e-19\nI(1 * (spray == \"B\"))   0.8333      1.601  0.5205 6.045e-01\nI(1 * (spray == \"C\")) -12.4167      1.601 -7.7550 7.267e-11\nI(1 * (spray == \"D\"))  -9.5833      1.601 -5.9854 9.817e-08\nI(1 * (spray == \"E\")) -11.0000      1.601 -6.8702 2.754e-09\nI(1 * (spray == \"F\"))   2.1667      1.601  1.3532 1.806e-01\n

\n\n\n\n

However, if we drop the intercept, then the Spray A term is\nno longer redundant. Then each coefficient is the mean for\nthat Spray.

\n\n

> summary(lm(count ~ spray - 1, data = InsectSprays))$coef\n       Estimate Std. Error t value  Pr(>|t|)\n\nsprayA   14.500      1.132  12.807 1.471e-19\nsprayB   15.333      1.132  13.543 1.002e-20\nsprayC    2.083      1.132   1.840 7.024e-02\nsprayD    4.917      1.132   4.343 4.953e-05\nsprayE    3.500      1.132   3.091 2.917e-03\nsprayF   16.667      1.132  14.721 1.573e-22\n

\n\n\n\n

So, for example, 14.5 is the mean for Spray A (as we already knew),\n15.33 is the mean for Spray B (14.5 + 0.8333 from our previous model formulation),\n2.083 is the mean for Spray C (14.5 - 12.4167 from our previous model formulation) and so on.\nThis is a nice trick if you want your model formulated in the terms of the group means,\nrather than the group comparisons relative to the reference group.

\n\n

Also, if there are no other covariates, the estimated coefficients for this model are\nexactly the empirical means of the groups. We can use dplyr to check this really\neasily and grab the mean for each group.

\n\n

> library(dplyr)\n> summarise(group_by(InsectSprays, spray), mn = mean(count))\nSource: local data frame [6 x 2]\n\n  spray     mn\n1     A 14.500\n2     B 15.333\n3     C  2.083\n4     D  4.917\n5     E  3.500\n6     F 16.667\n

\n\n\n\n

Often your lowest alphanumeric level isn’t the level that you’re most\ninterested in as a reference group. There’s an easy fix for that with\nfactor variables; use the relevel function. Here we give a simple\nexample. We created a variable spray2 that has Spray C as the reference\nlevel.

\n\n

> spray2 <- relevel(InsectSprays$spray, \"C\")\n> summary(lm(count ~ spray2, data = InsectSprays))$coef\n            Estimate Std. Error t value  Pr(>|t|)\n(Intercept)    2.083      1.132  1.8401 7.024e-02\nspray2A       12.417      1.601  7.7550 7.267e-11\nspray2B       13.250      1.601  8.2755 8.510e-12\nspray2D        2.833      1.601  1.7696 8.141e-02\nspray2E        1.417      1.601  0.8848 3.795e-01\nspray2F       14.583      1.601  9.1083 2.794e-13\n

\n\n\n\n

Now the intercept is the mean for Spray C and all of the coefficients are interpreted\nwith respect to Spray C. So, 12.417 is the comparison between Spray A and Spray C (as A - C)\nand so on.

\n\n

\nSummary of dummy variables

If you haven’t seen this before, it might seem rather strange. However, it’s essential to understand how\ndummy variables are treated, as otherwise huge interpretation errors can be made. Here we give a brief bullet summary\nof dummy variables to help solidify this information.

\n\n

If we treat a variable as a factor, R includes an intercept and omits the alphabetically first level of the factor.\n
- The intercept is the estimated mean for the reference level.
- The intercept t-test tests for whether or not the mean for the reference level is 0.
- All other t-tests are for comparisons of the other levels versus the reference level.
- Other group means are obtained the intercept plus their coefficient.
\n
If we omit an intercept, then it includes terms for all levels of the factor.\n
- Group means are now the coefficients.
- Tests are tests of whether the groups are different than zero.
\n
If we want comparisons between two levels, neither of which is the reference level, we could refit the model with one of them as the reference level.

\n\n

\nOther thoughts on this data

We don’t suggest that this is in anyway a thorough analysis of this data. For example,\nthe data are counts which are bounded from below by 0. This clearly violates the assumption of normality of the errors.\nAlso there are counts near zero, so both the actual assumption and the intent of this assumption are violated.\nFurthermore, the variance does not appear to be constant (look back at the violin plots).\nPerhaps taking logs of the counts would help. But, since there are 0 counts, maybe log(Count + 1).\nAlso, we’ll cover Poisson GLMs for fitting count data.

\n\n

\nFurther analysis of the `swiss` dataset

\n\n

\n\n

Then watch this video.

\n\n

Let’s create some dummy variables in the swiss dataset to illustrate them in a more\nmultivariable context. Just to remind ourselves of the dataset, here’s the first few\nrows.

\n\n

> library(datasets); data(swiss)\n> head(swiss)\n\n          Fertility Agriculture Examination Education Catholic Infant.Mortality\nCourtelary     80.2        17.0          15        12     9.96             22.2\nDelemont       83.1        45.1           6         9    84.84             22.2\nFranches-Mnt   92.5        39.7           5         5    93.40             20.2\nMoutier        85.8        36.5          12         7    33.77             20.3\nNeuveville     76.9        43.5          17        15     5.16             20.6\nPorrentruy     76.1        35.3           9         7    90.57             26.6\n

\n\n\n\n

Let’s create a binary variable out of the variable Catholic to illustrate dummy\nvariables in multivariable models. However, it should be noted that this isn’t\npatently absurd, since the variable is highly bimodal anyway. Let’s just split\nat majority Catholic or not:

\n\n

> library(dplyr)\n> swiss = mutate(swiss, CatholicBin = 1 * (Catholic > 50))\n

\n\n\n\n

Since we’re interested in Agriculture as a variable and Fertility as an outcome, let’s plot\nthose two color coded by the binary Catholic variable:

\n\n

g = ggplot(swiss, aes(x = Agriculture, y = Fertility, colour = factor(CatholicBin)))\ng = g + geom_point(size = 6, colour = \"black\") + geom_point(size = 4)\ng = g + xlab(\"% in Agriculture\") + ylab(\"Fertility\")\ng\n

\n\n\n\n\n

\n\n\n

Our model is:

\n\n\n$$\nY_i = \\beta_0 + X_{i1} \\beta_1 + X_{i2} \\beta_2 + \\epsilon_i\n$$\n\n

where \$Y_i\$ is Fertility, \$X_{i1}\$ is ‘Agriculture and\n\$X_{i2}\$ is CatholicBin. Let’s first fit the model with \$X_{i2}\$\nremoved.

\n\n

> summary(lm(Fertility ~ Agriculture, data = swiss))$coef\n\n            Estimate Std. Error t value  Pr(>|t|)\n(Intercept)  60.3044    4.25126  14.185 3.216e-18\nAgriculture   0.1942    0.07671   2.532 1.492e-02\n

\n\n\n\n

This model just assumes one line through the data (linear regression).\nNow let’s add our second variable. Notice that the model is

\n\n\n$$\nY_i = \\beta_0 + X_{i1} \\beta_1 + \\epsilon_{i}\n$$\n\n

when \$X_{i2} = 0\$ and

\n\n\n$$\nY_i = (\\beta_0 +\\beta_2) + X_{i1} \\beta_1 + \\epsilon_{i}\n$$\n\n

when \$X_{i2} = 1\$. Thus, the coefficient in front of the binary\nvariable is the change in the intercept between non-Catholic and Catholic\nmajority provinces. In other words, this model fits parallel lines\nfor the two levels of the factor variable. If the factor variable had\n4 levels, it would fit 4 parallel lines, where the coefficients for\nthe factors are the change in the intercepts to the reference level.

\n\n

## Parallel lines\nsummary(lm(Fertility ~ Agriculture + factor(CatholicBin), data = swiss))$coef\n\n                     Estimate Std. Error t value  Pr(>|t|)\n(Intercept)           60.8322     4.1059  14.816 1.032e-18\nAgriculture            0.1242     0.0811   1.531 1.329e-01\nfactor(CatholicBin)1   7.8843     3.7484   2.103 4.118e-02\n

\n\n\n\n

Thus, 7.8843 is the estimated\nchange in the intercept in the expected relationship between Agriculture and\nFertility going from a non-Catholic majority province to a Catholic majority.

\n\n

Often, however, we want both a different intercept and slope. This is easily\nobtained with an interaction term

\n\n\n$$\nY_i = \\beta_0 + X_{i1} \\beta_1 + X_{i2} \\beta_2 + X_{i1}X_{i2} \\beta_3 + \\epsilon_i .\n$$\n\n

Now consider with \$X_{i2} = 0\$, the model reduces to:

\n\n\n$$\nY_i = \\beta_0 + X_{i1} \\beta_1 + \\epsilon_{i}.\n$$\n\n

When \$X_{i2} = 1\$ the model is

\n\n\n$$\nY_i = (\\beta_0 + \\beta_2) + X_{i1} (\\beta_1 + \\beta_3) + \\epsilon_{i}.\n$$\n\n

Thus, the coefficient in front of the main effect \$X_{i2}\$, labeled \$\\beta_2\$ in our model, is the\nchange in the intercept, while the coefficient in front of\ninteraction term \$X_{i2}X_{i1}\$, labeled \$\\beta_3\$ in our model, is the\nchange in the slope. Let’s try it:

\n\n

> summary(lm(Fertility ~ Agriculture * factor(CatholicBin), data = swiss))$coef\n\n                                 Estimate Std. Error t value  Pr(>|t|)\n(Intercept)                      62.04993    4.78916 12.9563 1.919e-16\nAgriculture                       0.09612    0.09881  0.9727 3.361e-01\nfactor(CatholicBin)1              2.85770   10.62644  0.2689 7.893e-01\nAgriculture:factor(CatholicBin)1  0.08914    0.17611  0.5061 6.153e-01\n

\n\n\n\n

Thus, 2.8577 is the estimated change in the intercept of the linear relationship between Agriculture and\nFertility going from non-Catholic majority to\nCatholic majority to Catholic majority provinces. The interaction term, 0.08914, is the estimate change in\nthe slope. The estimated intercept in non-Catholic provinces is 62.04993 while the estimated intercept in\nCatholic provinces is 62.04993 + 2.85770. The estimated slope in non-Catholic majority provinces is 0.09612\nwhile it is 0.09612 + 0.08914 for Catholic majority provinces. If the factor has more than two levels, all\nof the main effects are change in the intercepts from the reference level while all of the interaction terms\nare changes in slope (again compared to the reference level).

\n\n

Homework exercise, plot both lines on the data to see the fit!

\n\n

\nExercises

\n\n

Do exercise 1 of the previous chapter if you have not already.\nLoad the dataset Seatbelts as part of the datasets package via data(Seatbelts). Use\nas.data.frame to convert the object to a dataframe. Fit a linear model of driver deaths\nwith kms and PetrolPrice as predictors. Interpret your results.
Repeat question 1 for the outcome being the log of the count of driver deaths. Interpret your coefficients. Watch a video solution.\n
Refer to question 1. Add the dummy variable law and interpret the results. Repeat this\nquestion with a factor variable that you create called lawFactor that takes the levels No and Yes.\nChange the reference level from No to Yes. Watch a video solution.\n
Discretize the PetrolPrice variable into four factor levels. Fit the linear model with this factor\nto see how R treats multiple level factor variables. Watch a video solution.\n
Perform the plot requested at the end of the last chapter.

\n\n\n

B0;256;0c# Adjustment

\n\n

\n\n

Adjustment, is the idea of putting regressors into a linear model to\ninvestigate the role of a third variable on the relationship between\nanother two. Since it is often the case that a third variable can\ndistort, or confound, the relationship between two others.

\n\n

As an example, consider looking at lung cancer rates and breath mint\nusage. For the sake of completeness, imagine if you were looking\nat forced expiratory volume (a measure of lung function) and breath\nmint usage. If you found a statistically significant regression relationship, it\nwouldn’t be wise to rush off to the newspapers with the headline\n“Breath mint usage causes shortness of breath!”, for a variety of reasons.\nFirst off, even if the association is sound, you don’t know that it’s\ncausal. But, more importantly in this case, the likely culprit\nis smoking habits. Smoking rates are likely related to both breath mint\nusage rates and lung function. How would you defend your finding\nagainst the accusation that it’s just variability in smoking habits?

\n\n

If your finding held up among non-smokers and smokers analyzed\nseparately, then you might have something. In other words, people\nwouldn’t even begin to believe this finding unless it held up\nwhile holding smoking status constant. That is the idea of\nadding a regression variable into a model as adjustment. The\ncoefficient of interest is interpreted as the effect of the\npredictor on the response, holding the adjustment variable\nconstant.

\n\n

In this chapter, we’ll use simulation to\ninvestigate how adding a regressor into a model addresses\nthe idea of adjustment.

\n\n

\nExperiment 1

\n\n

Let’s first generate some data. Consider the model

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X + \\tau T + \\epsilon_i\n$$\n\n

We’re interested in the relationship between our binary treatment,\n\$T\$, and\n\$Y\$. However, we’re concerned that the relationship\nmay depend on the continuous variable, \$X\$.

\n\n

Let’s simulate some data.

\n\n

n <- 100; t <- rep(c(0, 1), c(n/2, n/2)); x <- c(runif(n/2), runif(n/2));\nbeta0 <- 0; beta1 <- 2; tau <- 1; sigma <- .2\ny <- beta0 + x * beta1 + t * tau + rnorm(n, sd = sigma)\n

\n\n\n\n

Let’s plot the data. Below I give the code for the first plot; the rest\nof the code for plots throughout this chapter is omitted. (However,\nyou can see the course git repository for the rest of the code.)

\n\n

Simulation 1

plot(x, y, type = \"n\", frame = FALSE)\nabline(lm(y ~ x), lwd = 2)\nabline(h = mean(y[1 : (n/2)]), lwd = 3)\nabline(h = mean(y[(n/2 + 1) : n]), lwd = 3)\nfit <- lm(y ~ x + t)\nabline(coef(fit)[1], coef(fit)[2], lwd = 3)\nabline(coef(fit)[1] + coef(fit)[3], coef(fit)[2], lwd = 3)\npoints(x[1 : (n/2)], y[1 : (n/2)], pch = 21, col = \"black\", bg = \"lightblue\", cex = \\\n2)\npoints(x[(n/2 + 1) : n], y[(n/2 + 1) : n], pch = 21, col = \"black\", bg = \"salmon\", c\\\nex = 2)\n

\n\n\n\n\n

\n\n\n

Looking at this plot, notice that the X variable is\nunrelated to treatment/group status (color). In addition, the X variable\nis clearly linearly related to Y, but the intercept\nof this relationship depends on group status. The treatment variable is also\nrelated to Y; especially look at the horizontal lines\nwhich connect the group means to the Y axis. The third line is the\nwhat you would get if you just fit X and ignored group.\nFurthermore, notice that the relationship between group status and Y is constant depending on X.\nIn other words, both the apparent relationship and our estimated model have parallel lines. (Remember, our\n model, by not including an interaction term, did not allow for estimated non-parallel lines.)

\n\n

Finally, notice that the estimated relationship between the group variable and the outcome doesn’t\nchange much, regardless of whether X is accounted for or not. You can see this by comparing the\ndistance between the horizontal lines and the distance between the intercepts of the fitted lines. The horizontal lines are the group averages (disregarding X).\nThat the relationship doesn’t change much is ultimately a statement about balance. The nuisance variable (X) is well balanced between\nlevels of the group variable. So, whether you account for X or not, you get about the same answer. Moreover, we have lots of data at every level of\nX to make a direct comparison of the group on Y.\nOne way to try to achieve such balance with high probability is to randomize the group variable. This is especially useful,\nof course, when one doesn’t get to observe the nuisance covariate. Though be careful that as the number\nof unobserved covariates

\n\n

Now let’s consider less ideal settings.

\n\n

\nExperiment 2

\n\n\n

In this experiment, the X variable is highly related to group status.\nThat is, if you know the X variable, you could very easily predict\nwhich group they belonged to. If we disregard X, there’s an apparent\nstrong relationship between the group variable and Y. However, if\nwe account for X, there’s basically none. In this case, the apparent\neffect of group on Y is entirely explained by X. Our regression\nmodel would likely have a strong significant effect if group was\nincluded by itself and this effect would vanish if X was included.

\n\n

Further notice, there are no data to directly compare the groups\nat any particular value of X. (There’s no vertical overlap\nbetween the blue and red points.) Thus the adjusted effect\nis entirely based on the model, specifically the assumption\nof linearity. Try to draw curves on this plot assuming\nnon-linear relationships outside of their cloud of points\nfor the blue and red groups. You quickly will conclude\nthat many relationships are possible that would differ\nfrom this model’s conclusions. Worse still, you have no data to\ncheck the assumptions. Of course, R will churn forward without\nany complaints fitting this model and reporting no significant\ndifference between the groups.

\n\n

It’s worth noting at this point, that our experiments\njust show how the data can arrive at different effects\nwhen X is included or not. In a real application,\nit may be the case that X should be\nincluded and maybe that it shouldn’t be.

\n\n

For example,\nconsider an example that I was working on a few years ago. Imagine\nif group was whether or not the subject was\ntaking blood pressure medication and X was systolic blood pressure\n(ostensibly, the two variable giving the same information). It may\nnot make sense to adjust for blood pressure when looking at blood\npressure medication on the outcome.

\n\n

On the other hand consider another setting I ran into. A colleague was\nstudying chemical brain measurements of patients with a severe mental\ndisorder versus controls post mortem. However, the time since death\nwas highly related to the time the brain was stored since death,\nperhaps due to the differential patient sources of the two groups. The\ntime since death was strongly related to the outcome we were studying.\nIn this case, it is very hard to study the groups as they were so\ncontaminated by this nuisance covariate.

\n\n

Thus we arrive at the conclusion that whether or not to include\na covariate is a complex process relying on both the statistics and\na careful investigation into the subject matter being studied.

\n\n

\nExperiment 3

\n\n\n

In this experiment, we simulated data where the marginal (ignoring X)\nand conditional (using X) associations differ. First note that\nif X is ignored, one would estimate a higher marginal mean for Y\nof the red group over the blue group. However, if we look at the\nintercept in the fitted model, the blue group has a higher\nintercept. In other words, if you were to fit this linear\nmodel as lm(Y ~ Group) you would get one answer and\nlm(Y ~ Group + X) would give you the exact opposite answer,\nand in both cases the group effect would be highly statistically\nsignificant!

\n\n

Also in this settings, there isn’t a lot of overlap\nbetween the groups for any given X. That means there\nisn’t a lot of direct evidence to compare the groups\nwithout relying heavily on the model. In other words,\ngroup status is related to X quite strongly (though not\n as strongly as in the previous example). The adjusted\nrelationship suggest that the blue group is larger\nthan the red group. However, the reversal\nof the effect comes as\nbigger X means more likely red and bigger X means\nhigher Y.

\n\n

Let’s concoct an example around a way this data could\nhave occurred. Suppose\nthat you’re comparing two ad campaigns (labeled\nblue and red). Y is the sales from the ad (suppose you can measure this)\nand X is time of day that the ad is shown. Ads\nshown later on in the day do better than ads\nshown earlier. However, the blue ad campaign\ntended to get run in the morning while the\nred one tended to get run in the evening. So,\nignoring time of day leads to the erroneous\nconclusion that the the red ad did better. Again\nrandomization of the ads to time slots would\nlikely have eliminated this problem.

\n\n

\nExperiment 4

\n\n\n

Now that you’ve gotten the hang of it. You can\nsee how marginal and conditional associations\ncan differ. Experiment 4 is a case where the marginal association\nis minimal yet the conditional association is large.\nIn this case, by adding X to the model, the group\neffect became more statistically significant.

\n\n

\nExperiment 5

\n\n\n

\n\n\n

Let’s look at a weird one. In this case,\nthe best fitting model has both a group main\neffect and interaction with X. The main point here\nis that there is no meaningful group effect, the\neffect of group depends on what level of X you’re\nat. At a small value of X, the red group is here\nand at a large value of X, the blue group is higher;\nat intermediate values, they’re the same. Thus, it\nmakes no sense to talk about a group effect in this\nexample; group and X are intrinsically linked in their\nimpact on Y.

\n\n

As an example, imagine if Y is health outcome, X is time\nand group is two medications. One makes you much better\nright away then much worse as time goes on and the other\ndoesn’t do much at the start but steadily improves symptoms\nover time. Of course, most examples seen in practice aren’t that\nextreme. Still even with a slight departure in constant slopes,\nthe meaning of a main group effect goes away.

\n\n

\nSome final thoughts

\n\n

Nothing we’ve discussed is intrinsic to having a discrete group and\ncontinuous X. One, the other, both or neither could be discrete. What\nthis reinforces is that modeling multivariable relationships is hard.\nYou should continue to play around with simulations to see how the\ninclusion or exclusion of another variable can change apparent\nrelationships.

\n\n

We should also caution that our discussion only dealt with\nassociations. Establishing causal or truly mechanistic\nrelationships requires quite a bit more thinking.

\n\n

\nExercises

Load the dataset Seatbelts as part of the datasets package via data(Seatbelts). Use\nas.data.frame to convert the object to a dataframe. Fit a linear model of driver deaths\nwith kms and PetrolPrice as predictors. Interpret your results.
Compare the kms coefficient with and without the inclusion of the PetrolPrice variable in the model. Watch a video solution.\n
Compare the PetrolPrice coefficient with and without the inclusion of the kms variable in the model.

\n\n\n

\nResiduals, variation, diagnostics

\n\n

\n\n

\nResiduals

\n\n

Recall from Chapter 6 that the vertical distances between the observed data points and the fitted\nregression line are called residuals. We can generalize this idea to the vertical distances between\nthe observed data and the fitted surface in multivariable settings.

\n\n

To be specific, recall our linear model, which was\nspecified as \$Y_i = \\sum_{k=1}^p X_{ki} \\beta_j + \\epsilon_{i}\$. Throughout this lecture,\nwe’ll also assume that \$\\epsilon_i \\stackrel{iid}{\\sim} N(0, \\sigma^2)\$, even though this assumption isn’t\nnecessary for the definition of the residuals.

\n\n

We define the residuals as:

\n\n\n$$\ne_i = Y_i - \\hat Y_i = Y_i - \\sum_{k=1}^p X_{ki} \\hat \\beta_j.\n$$\n\n

This definition is identical (\$Y_i - \\hat Y_i\$) to our definition in the linear regression case.\nThe residuals are the vertical distances between the observed data points and the fitted regression surface.\nJust like in linear regression, our estimate of residual variation is basically the average of the squared\nresiduals. Specifically, \$\\hat \\sigma^2 = \\frac{\\sum_{i=1}^n e_i^2}{n-p}\$. Just like the before, we divide\nby \$n-p\$ rather than \$n\$ so that the estimate is unbiased, \$E[\\hat \\sigma^2] = \\sigma^2\$.

\n\n

Obtaining and plotting residuals in R is particularly easy. The function resid will return the residuals of a model\nfit with lm. Some useful plots, including a residual plot, can be performed with the plot function on the output\nof a lm fit. Consider the swiss dataset from previous chapters.

\n\n

data(swiss); par(mfrow = c(2, 2))\nfit <- lm(Fertility ~ . , data = swiss); plot(fit)\n

\n\n\n\n\n

\"The — The result of the method `plot` on the `swiss` dataset.

\n\n\n

Consider the upper left hand plot of the residuals (\$e_i\$) versus\nthe fitted values (\$\\hat Y_i\$). Often, a horizontal reference line at 0 is drawn\nsince (whenever an intercept is included) the residuals must sum to 0 and\nso will lie above and below the zero. Just like in our previous residual\nplots, one should look for any systematic patterns or large outlying observations.

\n\n

Note that this is one of many residual plots that one may be interested in\nperforming. For example, one might want to look at plots of residuals by\nindividual predictors or, as is done by plot, versus leverage (defined\n later in this chapter).

\n\n

\nInfluential, high leverage and outlying points

\n\n

As previously mentioned, it is a good idea to check our data for outliers.\nWe may want to refer back to these points to see if we can ascertain how\nthey became outliers, such as a misrecording. In addition, we may want\nto fit the models with and without those points included in order to\nascertain their influence on the model fit and inferential goals.

\n\n

Outliers can results for a variety of reasons. They can be real, but\ninconvenient, data. They could arise from spurious processes like\nprocessing or recording errors. They can have varying degrees of\ninfluence on our model. Thus, calling a point an outlier is vague and\nwe need a more precise language to discuss points that fall outside of\nour cloud of data. The plot below is useful for understanding different\nsorts of outliers.

\n\n\n

The lower left hand point is not an outlier having neither leverage\nnor influence on our fitted model. The upper left hand point is an\noutlier in the Y direction, but not in the X. It will have little\nimpact on our fitted model, since there’s lots of X points nearby to\ncounteract its effect. This point is said to have low leverage and\ninfluence. The upper right hand point is outside of the range of X\nvalues and Y values, but conforms nicely to the regression\nrelationship. This point has little effect on the fitted model. It has\nhigh leverage, but chooses not to exert it, and thus has low\ninfluence. The lower right hand point is outside of the range of X\nvalues, but not the Y values. However, it does not conform to the\nrelationship of the remainder of the points at all. This outlier has\nhigh leverage and influence.

\n\n

From this discussion you can maybe guess at the formal definition of\ntwo important terms: leverage and influence. Leverage discusses how\noutside of the norm a point’s X values are from the cloud of other X\nvalues. A point with high leverage has the opportunity to dramatically\nimpact the regression model. Whether or not it does so depends on how\nclosely it conforms to the fit.

\n\n

The other concept, influence, is a measure of how much impact a point has\non the regression fit. The most direct way to measure influence is\nfit the model with the point included and excluded.

\n\n

\nResiduals, Leverage and Influence measures

\n\n

\n\n

Now that we understand the three concepts of residuals, leverage and influence, we present a laundry list\nof probes. Do ?influence.measures to see the full suite of influence measures in stats.

\n\n

First consider residuals. We already know if fit is the output of\nlm (as in fit = lm(y ~ x1 + x2)), then resid(fit) returns the\nordinary residuals. A problem, though, is that these are defined as\n\$Y_i - \\hat Y_i \$ and thus have the units of the outcome. This\nisn’t great for comparing residual values across different analyses\nwith different experiments. So, some efforts to standardize residuals\nhave been made. Two of the most common are:

\n\n

\nrstandard - residuals divided by their standard deviations
\nrstudent - residuals divided by their standard deviations, where the ith data point was deleted in the calculation of the standard deviation for the residual to follow a t distribution

\n\n

Both of these endeavor to create T-like (as in Student’s T distribution)\nstatistics so that one can threshold residuals using T cutoffs. This is\nwhy these sorts of residuals are called studentized. The rstudent\nresiduals are exactly T distributed while the rstandard is not. The\nrstandard residuals are sometimes called internally standardized while the\nrstudent are called externally. The distinction between the residuals\nis mostly for establishing probability based cutoffs. Instead, we recommend\nlooking at the residuals as a collective and using the cutoffs loosely.\nUnder this way of thinking,\nthe distinctions over which of these two kinds of standardization are used\nis more academic than practical.

\n\n

A common use for residuals is to diagnose normality of the errors. This\nis often done by plotting the residual quantiles versus normal quantiles.\nThis is called a residual QQ plot. Your residuals should fall roughly\non a line if plotted in a normal QQ plot. There is of course noise and\na perfect fit would not be expected even if the model held.

\n\n

Leverage is largely measured by one quantity, so called hat diagonals, which can be obtained in R by the function hatvalues. The\nhat values are necessarily between 0 and 1 with larger values indicating\ngreater (potential for) leverage.

\n\n

After leverage, there are quite a few ways to probe for influence. These\nare:

\n\n

\ndffits - change in the predicted response when the \$i^{th}\$ point is deleted in fitting the model.
\ndfbetas - change in individual coefficients when the \$i^{th}\$ point is deleted in fitting the model.
\ncooks.distance - overall change in the coefficients when the \$i^{th}\$ point is deleted.

\n\n

In other words, the dffits check for influence in the fitted values,\ndfbetas check for influence in the coefficients individually and cooks.distance checks for influence in the coefficients as a collective.

\n\n

Finally, there’s a residual measure that’s also an influence measure.\nParticularly, consider resid(fit) / (1 - hatvalues(fit)) where fit is the linear model fit. This is the so-called PRESS residuals. These\nare the residual error from leave one out cross validation. That is, the difference in the response and the predicted response at data point \$i\$, where it was not included in the model fitting.

\n\n

\nHow do I use all of these things?

\n\n

First of all, be wary of simplistic rules for diagnostic plots and measures. The use of these tools is context specific. It’s better to understand what they are trying to accomplish and use them judiciously. Not all diagnostics measures have meaningful absolute scales. You can look at them relative to the values across the data. Even for the ones with known\nexact distributions to establish cutoffs,\nthose distributions (like the externally studentized residual) have degrees of freedom that depend on the sample size, so a single threshold can’t\nbe used across all settings.

\n\n

A better way to think about these tool is as diagnostics, like a\nphysician diagnosing a health issue. These tools probe your data in different ways to diagnose different problems. Some examples include:

\n\n

Patterns in your residual plots generally indicate some poor aspect of model fit.
Heteroskedasticity (non constant variance).
Missing model terms.
Temporal patterns (plot residuals versus collection order).
Residual QQ plots investigate normality of the errors.
Leverage measures (hat values) can be useful for diagnosing data entry errors and points that have a high potential for influence.
Influence measures get to the bottom line, ‘how does deleting or including this point impact a particular aspect of the model’.

\n\n

Let’s do some experiments to see how these measure hold up.

\n\n\n

\nSimulation examples

\n\n

\nCase 1

\n\n

In what follows, we’re going to try several simulation settings and see the values of some\non the residuals, influence measures and leverage. In our first\ncase, we simulate 100 points. The 101st point, c(10, 10),\nhas created a strong regression relationship where there shouldn’t be one.\nNote we prepend this point at the beginning of the Y and X vectors.

\n\n

n <- 100; x <- c(10, rnorm(n)); y <- c(10, c(rnorm(n)))\nplot(x, y, frame = FALSE, cex = 2, pch = 21, bg = \"lightblue\", col = \"black\")\nabline(lm(y ~ x))\n

\n\n\n\n

\n\n

First, let’s look at the dfbetas. Note the dfbetas are 101 by 2 dimensional, since there’s a dfbeta\nfor both the intercept and the slope. Let’s just output the first 10 for the slope.

\n\n

> fit <- lm(y ~ x)\n> round(dfbetas(fit)[1 : 10, 2], 3)\n     1      2      3      4      5      6      7      8      9     10\n 6.007 -0.019 -0.007  0.014 -0.002 -0.083 -0.034 -0.045 -0.112 -0.008\n

\n\n\n\n

Clearly the first point has a much, much larger dfbeta for the slope than\nthe other points. That is, the slope changes dramatically from when we include\nthis point to not including it. Try it out with Cook’s distance and the\ndffits. Let’s look at leverage.

\n\n

round(hatvalues(fit)[1 : 10], 3)\n    1     2     3     4     5     6     7     8     9    10\n0.445 0.010 0.011 0.011 0.030 0.017 0.012 0.033 0.021 0.010\n

\n\n\n\n

Again, this point has a much higher leverage value than that of\nthe other points. By having a large leverage value and\ndfbeta, we’re seeing that this point has a high potential for influence, and\ndecided to exert it.

\n\n

\nCase 2

\n\n

Consider a second case where the point lies on a natural\nline defined by the data, but is well outside of the cloud\nof X values. Since the code is so similar, I don’t show it.\nBut, as always, it can be found in the github repo for the\ncourses.

\n\n\n

Now let’s consider the dfbetas and the leverage for the\nfirst 10 observations.

\n\n

> round(dfbetas(fit2)[1 : 10, 2], 3)\n     1      2      3      4      5      6      7      8      9     10\n-0.072 -0.041 -0.007  0.012  0.008 -0.187  0.017  0.100 -0.059  0.035\n> round(hatvalues(fit2)[1 : 10], 3)\n    1     2     3     4     5     6     7     8     9    10\n0.164 0.011 0.014 0.012 0.010 0.030 0.017 0.017 0.013 0.021\n

\n\n\n\n

As we would expect, the dfbeta value for the first point is well with\nthe range of the other points. The leverage is much larger than the others.\nIn this case, the point has high leverage, but chooses not to exert it as influence.

\n\n

Play around with more simulation examples to get a feeling for what these\nmeasures do. This will help more than anything in understanding their\nvalue.

\n\n

\nExample described by Stefanski

\n\n

\n\n

We end with a really fun example from Stefanski in TAS 2007 volume 61.\nThis paper illustrates how a residual plot can unveil hidden treasures that would be nearly\nimpossible to detect with other kinds of plots. He has several examples on his website\nand we go through one. First, let’s read in the data and do a scatterplot matrix.

\n\n

dat <- read.table('http://www4.stat.ncsu.edu/~stefanski/NSF_Supported/Hidden_Images/\\\norly_owl_files/orly_owl_Lin_4p_5_flat.txt', header = FALSE)\npairs(dat)\n

\n\n\n\n\n

\n\n\n

It looks like a big mess of nothing. We can fit a model and get\nthat all of the variables are highly significant

\n\n

> summary(lm(V1 ~ . -1, data = dat))$coef\n   Estimate Std. Error t value  Pr(>|t|)\nV2   0.9856    0.12798   7.701 1.989e-14\nV3   0.9715    0.12664   7.671 2.500e-14\nV4   0.8606    0.11958   7.197 8.301e-13\nV5   0.9267    0.08328  11.127 4.778e-28\n

\n\n\n\n

Can we call it a day? Let’s check a residual plot.

\n\n

fit <- lm(V1 ~ . - 1, data = dat); plot(predict(fit), resid(fit), pch = '.')\n

\n\n\n\n\n

\n\n\n

There appears to be a pattern. The moral of the story here is that residual plots\ncan really hone in on systematic patterns in the data that are completely non-apparent\nfrom other plots.

\n\n

\nBack to the Swiss data

\n\n\n

\nExercises

\n\n

Load the dataset Seatbelts as part of the datasets package via data(Seatbelts). Use\nas.data.frame to convert the object to a dataframe. Fit a linear model of driver deaths\nwith kms, PetrolPrice and law as predictors.
Refer to question 1. Directly estimate the residual variation via the function resid. Compare with R’s residual variance estimate. Watch a video solution.\n
Refer to question 1. Perform an analysis of diagnostic measures including, dffits, dfbetas, influence and hat diagonals. Watch a video solution.\n

\n\n

LocalWords: shouldn prepend lang rnorm cex pch bg Stefanski TAS scatterplot\n LocalWords: dat lm V1 coef V2 989e V3 500e V4 301e V5 778e resid swiss

\n\n\n

\nMultiple variables and model selection

\n\n

This chapter represents a challenging question: “How do we choose what\nto variables to include in a regression model?”. Sadly, no single easy\nanswer exists and the most reasonable answer would be “It depends.”.\nThese concepts bleed into ideas of machine learning, which is largely\nfocused on high dimensional variable selection and weighting. In this\nchapter we cover some of the basics and, most importantly, the\nconsequences of over- and under-fitting a model.

\n\n

\nMultivariable regression

\n\n

In our Coursera Data Science Specialization, we have an entire class\non prediction and machine learning. So, in this class, our focus will\nbe on modeling. That is, our primary concern is winding up with an\ninterpretable model, with interpretable coefficients. This is a very\ndifferent process than if we only care about prediction or machine\nlearning. Prediction has a different set of criteria, needs for\ninterpretability and standards for generalizability. In modeling, our\ninterest lies in parsimonious, interpretable representations of the\ndata that enhance our understanding of the phenomena under study.

\n\n

Like nearly all aspects of statistics, good modeling decisions are\ncontext dependent. Consider a good model for prediction, versus one\nfor studying mechanisms, versus one for trying to establish causal\neffects. There are, however, some principles to help you guide your\nway.

\n\n

Parsimony is a core concept in model selection. The idea of\nparsimony is to keep your models as simple as possible (but no\nsimpler). This builds on the idea of Occam’s\nrazor, in that all else\nbeing equal simpler explanations are better than complex ones. Simpler\nmodels are easier to interpret and are less finicky. Complex models\noften have issues with fitting and, especially, overfitting. (To see\na counterargument, consider Andrew Gelman’s\nblog..)

\n\n

Another principle that I find useful for looking at statistical models\nis to consider them as lenses through which to look at your data. (I\nattribute this quote to the great statistician Scott Zeger.) Under\nthis philosophy, what’s the right model - whatever one connects the\ndata to a true, parsimonious statement about what you’re\nstudying. Unwin and authors have formalized these ideas more into\nsomething they call exploratory model\nanalysis\nI like this, as it turns our focus away from trying to get a single,\nbest, true model and instead focuses on utilizing models as ways to\nprobe data. This is useful, since all models are wrong in some fashion.
\nKeep this in mind as we focus on variable inclusion and exclusion in this chapter.

\n\n

\nThe Rumsfeldian triplet

\n\n

Before we begin, I’d like to give a quote from Donald Rumsfeld, the\ncontroversial Secretary of Defense of the US during the start of the\nAfghanistan and second Iraq wars. He gave this quote regarding weapons\nof mass destruction (read more about it\nhere):

\n\n

“There are known knowns. These are things we know that we know. There\nare known unknowns. That is to say, there are things that we know we\ndon’t know. But there are also unknown unknowns. There are things we\ndon’t know we don’t know.” - Donald Rumsfeld

\n\n

This quote, widely derided for its intended purpose, is quite\ninsightful in the unintended context of regression model\nselection. Specifically, in our context “Known Knowns” are regressors\nthat we know we should check for inclusion in the model. The\n“Known Unknowns” are regressors that we would like to include in the\nmodel, but don’t have. The “Unknown Unknowns” are regressors that we\ndon’t even know about that we should have included in the model.

\n\n

In this chapter, we’ll talk about Known Knowns; variables that are\npotentially of interest in our model that we have. Known Unknowns and\nUnknown Unknowns (especially) are more challenging to deal with. A\ncentral method for dealing with Unknown Unknowns is randomization. If\nyou’d like to compare a treatment to a control, or perform an A/B test\nof two advertising strategies, randomization will help insure that\nyour treatment is balanced across levels of the Unknown Unknowns with\nhigh probability. (Of course, being unobserved, you can never know\nwhether or not the randomization was effective.)

\n\n

For Known Unknowns, those variables we wish we had collected but did\nnot, there are several strategies. For example, a proxy variable\nmight be of use. As an example, we had some brain volumetric\nmeasurements via MRIs and really wished we had done the processing to\nget intra-cranial volume (head size). The need for this variable was\nbecause we didn’t want to compare brain volumetric measurements and\nconclude that bigger people with bigger heads have more brain\nmass. This would be a useless conclusion. For example, killer whales\nhave bigger brains than dolphins, but that doesn’t tell you much about\nkiller whales or dolphins. More interesting would be if killer whales\nwho were exposed to toxic chemicals had lower brain volume relative to\ntheir intra-cranial volume than whales who weren’t. In our case, (we\nwere studying humans), we used height, gender and other\nanthropomorphic measurements to get a good guess of intra-cranial\nvolume.

\n\n

For the rest of the lecture, let’s discuss the known knowns and what\ntheir unnecessary inclusion and exclusion implies in our analysis.

\n\n

\nGeneral rules

Here we state a couple of general rules regarding model selection for our known knowns.

\n\n

Omitting variables results in bias in the coefficients of interest - unless the regressors are uncorrelated with the omitted ones.

\n\n

I want to reiterate this point: if the omitted variable is uncorrelated with the included variables, its omission has no impact\non estimation. It might explain some residual variation, thus it could have an impact on inference.\nAs previously mentioned, this lack of impact of uncorrelated variables is why we randomize treatments;\nrandomization attempts to disassociate our treatment indicator with variables that we don’t have to put in the model. Formal\ntheories of inference can be designed around the use of randomization. However,\nin a practical sense, if there’s too many unobserved confounding variables, even randomization won’t help you, since\nwith high probability one will stay correlated with the treatment.

\n\n

In most cases we won’t have randomization. So, to avoid bias, why don’t we throw everything into the regression model?\nThe following rule prevents us from doing that:

\n\n

Including variables that we shouldn’t have increases standard errors of the regression variables.

\n\n

Actually, including any new variables increases the actual (not estimated) standard errors of other regressors.\nSo we don’t want to idly throw variables into the model. In addition the model must tend toward perfect\nfit as the number of non-redundant regressors approaches the sample size. Our \$R^2\$\nincreases monotonically as more regressors are included, even unrelated white noise.

\n\n

\nR squared goes up as you put regressors in the model

\n\n

Let’s try a simulation. In this simulation, no regression relationship exists. We simulate data and \$p\$ regressors as random normals.\nThe plot is of the \$R^2\$.

\n\n

n <- 100\nplot(c(1, n), 0 : 1, type = \"n\", frame = FALSE, xlab = \"p\", ylab = \"R^2\")\ny <- rnorm(n); x <- NULL; r <- NULL\nfor (i in 1 : n){\n   x <- cbind(x, rnorm(n))\n   r <- c(r, summary(lm(y ~ x))$r.squared)\n}\nlines(1 : n, r, lwd = 3)\nabline(h = 1)\n

\n\n\n\n\n

\n\n\n

Notice that the \$R^2\$ goes up, monotonically, as the number of regressors\nis increased. This reminds us of a couple of things. First, irrelevant variables\nexplain residual variation by chance. And, when evaluating fit, we have to take\ninto account the number of regressors included. The adjusted \$R^2\$ is better\nfor these purposes than \$R^2\$ since it accounts for the number of variables\nincluded in the model. In R, you can get the adjusted \$R^2\$ very easily\nby grabbing summary(fitted_model)$adj.r.squared instead of summary(fitted_model)$r.squared.

\n\n

\nSimulation demonstrating variance inflation

\n\n

Now let’s use simulation to demonstrate variation inflation. In this case,\nwe’re going to simulate three regressors, x1, x2 and x3. We then repeatedly\ngenerate data from a model, where y only depends on x1. We fit three models,\ny ~ x1, y ~ x1 + x2, and y ~ x1 + x2 + x3. We do this over and over\nagain and look at the standard deviation of the x1 coefficient.

\n\n

> n <- 100; nosim <- 1000\n> x1 <- rnorm(n); x2 <- rnorm(n); x3 <- rnorm(n);\n> betas <- sapply(1 : nosim, function(i){\n  y <- x1 + rnorm(n, sd = .3)\n  c(coef(lm(y ~ x1))[2],\n    coef(lm(y ~ x1 + x2))[2],\n    coef(lm(y ~ x1 + x2 + x3))[2])\n  })\n> round(apply(betas, 1, sd), 5)\n\n     x1      x1      x1\n0.02839 0.02872 0.02884\n

\n\n\n\n

Notice that the standard error for the x1 coefficient goes up\nas more regressors are included (left to right in our vector output). It’s important to note that\nthese are the actual standard errors (obtained by repeatedly simulating\nthe data). These aren’t obtainable in a single dataset since we only get one realization.\nThe estimated standard errors, the ones we have access to in a data analysis,\nmay not go up as you include more regressors.

\n\n

Now let’s see if we can make the variance inflation worse.\nIn this case, I’ve made x2 and x3 correlated with x1.

\n\n

> n <- 100; nosim <- 1000\n> x1 <- rnorm(n); x2 <- x1/sqrt(2) + rnorm(n) /sqrt(2)\n> x3 <- x1 * 0.95 + rnorm(n) * sqrt(1 - 0.95^2);\n> betas <- sapply(1 : nosim, function(i){\n  y <- x1 + rnorm(n, sd = .3)\n  c(coef(lm(y ~ x1))[2],\n    coef(lm(y ~ x1 + x2))[2],\n    coef(lm(y ~ x1 + x2 + x3))[2])\n})\n> round(apply(betas, 1, sd), 5)\n     x1      x1      x1\n0.03131 0.04270 0.09653\n

\n\n\n\n

Notice that the variance inflation goes up quite a bit more.\nThis is an issue with including variables that are highly correlated\nwith the ones that we are interested in. In the first simulation,\nthe regressors were simulated independently, and the variance inflation\nwasn’t bad. In the second, they were correlated and it was much worse.

\n\n

\nSummary of variance inflation

Notice variance inflation was much worse when we included a variable that\nwas highly related to x1.
We don’t know \$\\sigma\$, the residual variance, so we can’t know the actual variance inflation amount.
However, \$\\sigma\$ drops out of the ratio of the standard errors.\nThus, if one sequentially adds variables, one can check the variance (or sd) inflation for including each one.
When the other regressors are actually orthogonal (correlation 0) to the regressor of interest, then there is no variance inflation.
The variance inflation factor (VIF) is the increase in the variance for the ith regressor compared to the ideal setting where it is orthogonal to the other regressors.\n
- The square root of the VIF is the increase in the sd instead of variance.
\n
Remember, variance inflation is only part of the picture. We want to include certain variables, even if they dramatically inflate our variance.

\n\n

Let’s revisit our previous simulation to show how one can estimate the relative\nincrease in variance. Let’s simulate a single dataset, and I’ll show how to get\nthe relative increase in variance for including x2 and x3. All you need to do is\ntake the ratio of the variances for that coefficient. If you don’t exactly\nunderstand the code, don’t worry. The idea is that we can obtain these from an observed\ndata set.

\n\n

> y <- x1 + rnorm(n, sd = .3)\n> a <- summary(lm(y ~ x1))$cov.unscaled[2,2]\n> c(summary(lm(y ~ x1 + x2))$cov.unscaled[2,2],\n  summary(lm(y~ x1 + x2 + x3))$cov.unscaled[2,2]) / a\n[1] 1.895 9.948\n

\n\n\n\n

Now let’s check it by referring to our previous simulation and see what\nthe relative variance for x1 is when including the x2 and x2 plus x3 models.

\n\n

> temp <- apply(betas, 1, var); temp[2 : 3] / temp[1]\n   x1    x1\n1.860 9.506\n

\n\n\n\n

Notice that it’s the same (about). In other words, from a single observed dataset\nwe can perfectly estimate the relative variance inflation caused by adding a regressor.

\n\n

\nSwiss data revisited

\n\n

> data(swiss);\n> fit1 <- lm(Fertility ~ Agriculture, data = swiss)\n> a <- summary(fit1)$cov.unscaled[2,2]\n>fit2 <- update(fit1, Fertility ~ Agriculture + Examination)\n> fit3 <- update(fit1, Fertility ~ Agriculture + Examination + Education)\n> c(summary(fit2)$cov.unscaled[2,2],\n    summary(fit3)$cov.unscaled[2,2]) / a\n[1] 1.892 2.089\n

\n\n\n\n

Thus inclusion of Examination increases the variance of the Agriculture effect by 89.2% while\nfurther adding Examination and Education causes a 108.9% increase. Again,\nthe observed standard errors won’t follow these percentages. These are the\nincreases if we actually knew \$\\sigma^2\$.

\n\n

Let’s look at the variance inflation factors. These measure how much\nvariance inflation the variable causes relative to the setting where\nit was orthogonal to the other regressors. This is nice because it\nhas a well contained interpretation within a single model fit. Also,\none doesn’t have to do all of the model refitting we did above to\nexplore variance inflation. So, in general, the VIFs are the most\nconvenient entity to work with.

\n\n

> library(car)\n> fit <- lm(Fertility ~ . , data = swiss)\n> vif(fit)\n     Agriculture      Examination        Education         Catholic Infant.Mortality\n           2.284            3.675            2.775            1.937            1.108\n> sqrt(vif(fit)) #If you prefer sd inflation\n     Agriculture      Examination        Education         Catholic Infant.Mortality\n           1.511            1.917            1.666            1.392            1.052\n

\n\n\n\n

\nImpact of over- and under-fitting on residual variance estimation

\n\n

\n\n

Assuming that the model is linear with additive iid errors,\nwe can mathematically describe the impact of omitting necessary variables or including unnecessary ones. These\ntwo rules follow:

\n\n

If we underfit the model, that is omit necessary covariates, the variance estimate is biased.
If we correctly, or overfit, the model, including all necessary covariates and possibly some unnecessary ones,\nthe variance estimate is unbiased. However, the variance of the variance is larger if we include unnecessary variables.

\n\n

These make sense. If we’ve omitted important variables, we’re attributing residual variation that is really systematic\nvariation explainable by those omitted covariates. Therefore, we would expect a variance estimate that is systematically\noff (biased). We would also expect absence of bias when we throw the kitchen sink at the model and include everything (necessary and unnecessary).\nHowever, then our variance estimate is unstable (the variance of the variance estimate is larger).

\n\n

\nCovariate model selection

Ideally, you include only the necessary variables in a regression model.\nHowever, it’s impossible to know in practice which ones are necessary\nand which ones are not. Thus we have to discuss variable selection a little bit.\nAutomated covariate selection is a difficult topic. It depends heavily on how rich of a covariate space one wants to explore.\nThe space of models explodes quickly as you add interactions and polynomial terms.

\n\n

In the Data Science Specialization prediction class, we’ll cover many modern methods for traversing large model spaces for the purposes of\nprediction. In addition, principal components or factor analytic models on covariates are often useful for reducing complex covariate spaces.

\n\n

It should also be noted that careful design can often eliminate the\nneed for complex model searches at the analyses stage. For example,\nrandomized designs, randomized block designs, crossover designs, clinical\ntrials, A/B testing are all examples of designs where randomization,\nbalance and stratification are used to create data sets that have more\ndirect analyses. However, control over the design is often limited in\ndata science.

\n\n

I’ll give my favorite approach for model selection when I’m trying to get a parsimonious explanatory model. (I would use a different strategy\n for prediction.) Given a coefficient that I’m interested in, I like to use covariate adjustment and multiple models to probe that effect to evaluate it for robustness and to see what other covariates knock it out or amplify it. In other words, if I have an effect, or absence of an effect, that I’d like to\n report, I try to first come up with criticisms of that effect and then use models to try to answer those criticisms.

\n\n

As an example, if I had a significant effect of lead exposure on brain size I would think about the following criticism. Were the high exposure\npeople smaller than the low exposure people. To address this, I would consider adding head size (intra-cranial volume). If the lead exposed were\nmore obese than the non-exposed, I would put a model with body mass index (BMI) included. This isn’t a terribly systematic approach,\nbut it tends to teach you a lot about the the data as you get your hands dirty. Most importantly, it makes you think hard about the questions\nyour asking and what are the potential criticisms to your results. Heading those criticisms off at the pass early on is a good idea.

\n\n

\nHow to do nested model testing in R

\n\n

One particular model selection technique is so useful I’ll cover it since it likely wouldn’t be covered in a machine learning or prediction class.\nIf the models of interest are nested and without lots of parameters differentiating them,\nit’s fairly uncontroversial to use nested likelihood ratio tests for model selection. Consider the following example:

\n\n

> fit1 <- lm(Fertility ~ Agriculture, data = swiss)\n> fit3 <- update(fit1, Fertility ~ Agriculture + Examination + Education)\n> fit5 <- update(fit1, Fertility ~ Agriculture + Examination + Education + Catholic \\\n+ Infant.Mortality)\n> anova(fit1, fit3, fit5)\n\nAnalysis of Variance Table\n\nModel 1: Fertility ~ Agriculture\nModel 2: Fertility ~ Agriculture + Examination + Education\nModel 3: Fertility ~ Agriculture + Examination + Education + Catholic +\n    Infant.Mortality\n  Res.Df  RSS Df Sum of Sq    F  Pr(>F)\n1     45 6283\n2     43 3181  2      3102 30.2 8.6e-09 ***\n3     41 2105  2      1076 10.5 0.00021 ***\n---\nSignif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n

\n\n\n\n

Notice how the three models I’m interested in are nested. That is, Model 3 contains all of the Model 2 variables which contains all of the Model 1\nvariables. The P-values are for a test of whether all of the new variables are all zero or not (i.e. whether or not they’re necessary). So this model would conclude that all of the added Model 3 terms are necessary over Model 2 and all of the Model 2 terms are necessary over Model 1. So, unless there were some other\ncompelling reasons, we’d pick Model 3.\nAgain, you don’t want to blindly follow a model selection procedure, but when the models are naturally nested, this is a reasonable approach.

\n\n

\nExercises

Load the dataset Seatbelts as part of the datasets package via data(Seatbelts). Use\nas.data.frame to convert the object to a dataframe. Fit a linear model of driver deaths\nwith kms, PetrolPrice and law as predictors.
Perform a model selection exercise to arrive at a final model. Watch a video solution.\n

\n\n

LocalWords: knowns regressors volumetric MRIs intra

\n\n\n

\nGeneralized Linear Models

\n\n

\n\n

Generalized linear models (GLMs) were a great advance in statistical\nmodeling. The original manuscript with the GLM framework was from\nNelder and Wedderburn in 1972.\nin the Journal of the Royal Statistical Society. The McCullagh and\nNelder book¹ is the famous standard treatise on the subject.

\n\n

Recall linear models. Linear models are the most useful applied\nstatistical technique. However, they are not without their\nlimitations. Additive response models don’t make much sense if the\nresponse is discrete, or strictly positive. Additive error models\noften don’t make sense, for example, if the outcome has to be\npositive. Transformations, such as taking a cube root of a count\noutcome, are often hard to interpret.

\n\n

In addition, there’s value in modeling the data on the scale that it\nwas collected. Particularly interpretable transformations, natural\nlogarithms in specific, aren’t applicable for negative or zero values.

\n\n

The generalized linear model is a family of models that includes linear\nmodels. By extending the family, it handles many of the issues with\nlinear models, but at the expense of some complexity and loss of some\nof the mathematical tidiness. A GLM involves three components:

\n\n

An exponential family model for the response.
A systematic component via a linear predictor.
A link function that connects the means of the response to the linear predictor.

\n\n

The three most famous cases of GLMs are: linear models, binomial and\nbinary regression and Poisson regression. We’ll go through the GLM\nmodel specification and likelihood for all three. For linear models,\nwe’ve developed them throughout the book. The next two chapters will\nbe devoted to binomial and Poisson regression. We’ll only focus on the\nmost popular and useful link functions.

\n\n

\nExample, linear models

Let’s go through an example. Assume that our response is {$$}Y_i \\sim\nN(\\mu_i, \\sigma^2){/$$}. The Gaussian distribution is an exponential\nfamily distribution. Define the linear predictor to be

\n\n

\$\\eta_i = \\sum_{k=1}^p X_{ki} \\beta_k\$.

\n\n

The link function as \$g\$ so that \$g(\\mu) = \\eta\$.\nFor linear models \$g(\\mu) = \\mu\$ so that \$\\mu_i = \\eta_i\$.\nThis yields the same likelihood model as our additive error Gaussian linear model

\n\n\n$$\nY_i = \\sum_{k=1}^p X_{ki} \\beta_k + \\epsilon_{i}\n$$\n\n

where \$\\epsilon_i \\stackrel{iid}{\\sim} N(0, \\sigma^2)\$. So, we’ve specified our\nmodel as a GLM above and with a more traditional linear model specification below.\nLet’s try an example where the GLM is more necessary.

\n\n

\nExample, logistic regression

\n\n

Assume that our outcome is a 0, 1 variable. Let’s model \$Y_i \\sim Bernoulli(\\mu_i)\$\nso that \$E[Y_i] = \\mu_i\$ where \$0\\leq \\mu_i \\leq 1\$.

\n\n

Linear predictor: \$\\eta_i = \\sum_{k=1}^p X_{ki} \\beta_k\$
Link function \$g(\\mu) = \\eta = \\log\\left( \\frac{\\mu}{1 - \\mu}\\right)\$

\n\n

In this case, \$g\$ is the (natural) log odds, referred to as the logit.\nNote then we can invert the logit function as:

\n\n\n$$\n\\mu_i = \\frac{\\exp(\\eta_i)}{1 + \\exp(\\eta_i)} ~~~\\mbox{and}~~~\n1 - \\mu_i = \\frac{1}{1 + \\exp(\\eta_i)}\n$$\n\n

Some people like to call this the expit function. The logit is useful as it\nconverts probabilities which lie in [0,1] into the whole real line, a more\nnatural space for the linear part of the model to live. Notice further,\nwe’re not transforming the outcome (Y). Instead, we’re modeling our Y\nas if it were a collection of coin flips and applying the transformation\nto the probability of a head.

\n\n

To get the estimates we maximize the likelihood. We can write out the likelihood as:

\n\n\n$$\n\\prod_{i=1}^n \\mu_i^{y_i} (1 - \\mu_i)^{1-y_i}\n= \\exp\\left(\\sum_{i=1}^n y_i \\eta_i \\right)\n\\prod_{i=1}^n (1 + \\eta_i)^{-1}\n$$\n\n

\nExample, Poisson regression

Let’s consider a problem with count data.\nAssume that :

\n\n

\$Y_i \\sim Poisson(\\mu_i)\$ so that \$E[Y_i] = \\mu_i\$ where \$0\\leq \\mu_i\$.
Linear predictor \$\\eta_i = \\sum_{k=1}^p X_{ki} \\beta_k\$.
Link function \$g(\\mu) = \\eta = \\log(\\mu)\$

\n\n

Recall that \$e^x\$ is the inverse of \$\\log(x)\$ so that we have:

\n\n\n$$\n\\mu_i = e^{\\eta_i}\n$$\n\n

Thus, the likelihood is:

\n\n\n$$\n\\prod_{i=1}^n (y_i !)^{-1} \\mu_i^{y_i}e^{-\\mu_i}\n\\propto \\exp\\left(\\sum_{i=1}^n y_i \\eta_i - \\sum_{i=1}^n \\mu_i\\right)\n$$\n\n

\nHow estimates are obtained

\n\n

For GLMs, estimates have to be obtained numerically through an iterative algorithm. The\nalgorithms are very well behaved, so convergence is usually not a problem unless you have\na lot of data on a boundary, such as a lot of 0 counts in binomial or Poisson data.\nThe standard errors are obtained also numerically, and are usually based on large sample\ntheory. The exact equation that gets solved is the so-called normal equations

\n\n\n$$\n0=\\sum_{i=1}^n \\frac{(Y_i - \\mu_i)}{Var(Y_i)}W_i\n$$\n\n

The variance differs by the model. The \$W_i\$ are derivative terms that we won’t deal with.

\n\n

For the linear model \$Var(Y_i) = \\sigma^2\$ (is constant).
For Bernoulli case \$Var(Y_i) = \\mu_i (1 - \\mu_i)\$
For the Poisson case \$Var(Y_i) = \\mu_i\$.

\n\n

In the latter two cases, it is often relevant to have a more flexible variance model, even if it doesn’t correspond to an actual likelihood.\nWe might make the following changes:

\n\n\n$$\n0=\\sum_{i=1}^n \\frac{(Y_i - \\mu_i)}{\\phi \\mu_i (1 - \\mu_i ) } W_i ~~~\\mbox{and}~~~\n0=\\sum_{i=1}^n \\frac{(Y_i - \\mu_i)}{\\phi \\mu_i} W_i\n$$\n\n

These are called ‘quasi-likelihood’ normal equations. R offers these as an option in the glm function\nas the quasipoisson and quasibinomial options. These offer more flexible variance options than straight\nPoisson and binomial models.

\n\n

\nOdds and ends

At this point, let’s do some bookkeeping before we work through examples.

\n\n

The normal equations have to be solved iteratively, resulting in\n\$\\hat \\beta_k\$ and, if included, \$\\hat \\phi\$.
Predicted linear predictor responses can be obtained as \$\\hat \\eta = \\sum_{k=1}^p X_k \\hat \\beta_k\$
Predicted mean responses as \$\\hat \\mu = g^{-1}(\\hat \\eta)\$
Coefficients are interpreted as

\n\n\n$$\ng(E[Y | X_k = x_k + 1, X_{\\sim k} = x_{\\sim k}]) - g(E[Y | X_k = x_k, X_{\\sim k}=x_{\\sim k}]) = \\beta_k\n$$\n\n

or the change in the link function of the expected response per unit change in \$X_k\$ holding other regressors constant.

\n\n

Variations on Newon/Raphson’s algorithm are used to do it.
Asymptotics are used for inference usually (but not always).
Many of the ideas from linear models can be brought over to GLMs.

\n\n

\nExercises

True or false, generalized linear models transform the observed outcome. (Discuss.) Watch a video solution.\n
True or false, the interpretation of the coefficients in a GLM are on the scale of the link function. (Discuss.) Watch a video solution.\n
True or false, the generalized linear model assumes an exponential family for the outcome. (Discuss.) Watch a video solution.\n
True or false, GLM estimates are obtained by maximizing the likelihood. (Discuss.) Watch a video solution.\n
True or false, some GLM distributions impose restrictions on the relationship between the mean and the variance. (Discuss.) Watch a video solution.\n

\n\n\n

\nBinary GLMs

\n\n

\n\n

Binary GLMs come from trying to model outcomes that can take only two\nvalues. Some examples include: survival or not at the end of a study,\nwinning versus losing of a team and success versus failure of a treatment or\nproduct. Often these outcomes are called Bernoulli outcomes, from the\nBernoulli distribution named after the famous probabilist and mathematician.

\n\n

If we happen to have several exchangeable binary outcomes for the same level\nof covariate values, then that is binomial data and we can aggregate the 0’s and\n1’s into the count of 1’s. As an example, imagine if we sprayed insect pests\nwith 4 different pesticides and counted whether they died or not. Then for\neach spray, we could summarize the data with the count of dead and total number\nthat were sprayed and treat the data as binomial rather than Bernoulli.

\n\n

\nExample Baltimore Ravens win/loss

\n\n

The Baltimore Ravens are an American Football team in\nthe US’s National Football League.² The data contains the wins\nand losses of the Ravens by the number of points that they scored. (In American football, the\ngoal is to score more points than your opponent.) It should be clear\nthat there would be a positive relationship between the number of\npoints scored and the probability of winning that particular game.

\n\n

Let’s load the data and use head to look at the first few rows.

\n\n

> download.file(\"https://dl.dropboxusercontent.com/u/7710864/data/ravensData.rda\"\n              , destfile=\"./data/ravensData.rda\",method=\"curl\")\n> load(\"./data/ravensData.rda\")\n> head(ravensData)\n  ravenWinNum ravenWin ravenScore opponentScore\n1           1        W         24             9\n2           1        W         38            35\n3           1        W         28            13\n4           1        W         34            31\n5           1        W         44            13\n6           0        L         23            24\n

\n\n\n\n

A linear regression model would look something like this:

\n\n\n$$\nY_i = \\beta_0 + \\beta_1 X_i + e_i\n$$\n\n

Where \$Y_i\$ is a binary indicator of whether or\nnot the Ravens won game \$i\$(1 for a win, 0 for a loss).\n\$X_i\$ is the number of points that they scored for that game.\nand \$\\epsilon_i\$ is the residual error term.

\n\n

Under this model then \$\\beta_0\$ is the expected value of\n\$Y_i\$ given a 0 point game. For a 0/1 variable, the expected\nvalue is the probability, so the intercept is the probability that\nthe Ravens win with 0 points scored. Then \$\\beta_1\$ is the\nincrease in probability of a win for each additional point.

\n\n

At this point in the book, I hope that fitting and interpreting\nthis model would be second nature.

\n\n

> lmRavens <- lm(ravensData$ravenWinNum ~ ravensData$ravenScore)\n> summary(lmRavens)$coef\n                      Estimate Std. Error t value Pr(>|t|)\n(Intercept)             0.2850   0.256643   1.111  0.28135\nravensData$ravenScore   0.0159   0.009059   1.755  0.09625\n

\n\n\n\n

There are numerous problems with this model. First, if the Ravens\nscore more than 63 points in a game, we estimate a 0.0159 *\n63, which is greater than 1, increase in the probability of them winning. This is\nan impossibility, since a probability can’t be greater than 1.\nSixty three is an unusual, but not impossible, score in American football, but\nthe principle applies broadly: modeling binary data with linear models\nresults in models that fail the basic assumption of the data.

\n\n

Perhaps less galling, but still an annoying aspect of the model, is that\nif the error is assumed to be Gaussian, then our model allows \$Y_i\$\nto be anything from minus infinity to plus infinity, when we know our\ndata can be only be 0 or 1. If we assume that our errors are discrete\nto force this, we assume a very strange distribution on the errors.

\n\n

There also aren’t any transformations to make things better. Any\none to one transformation of our outcome\nis still just going to have two values, thus the same set of problems.

\n\n

The key insight was to transform the probability of a 1 (win in our\n example) rather than the data itself. Which transformation is most\nuseful? It turns out that it involves the log of the odds, called\nthe logit.

\n\n

\nOdds

You’ve heard of odds before, most likely from discussions of gambling.\nFirst note, odds are a fraction greater than 0 but unbounded.\nThe odds are not a percentage or proportion. So, when someone says\n“The odds are fifty percent”, they are mistaking probability and odds.\nThey likely mean “The probability is fifty percent.”, or equivalently\n“The odds are one.”, or “The odds are one to one”, or “The odds are\nfifty [divided by] fifty.” The latter three odds statements\nare all the same since: 1, 1 / 1 and 50 / 50 are all the same number.

\n\n

If \$p\$ is a probability, the odds are defined as\n\$o = p/(1-p)\$. Note that we can go backwards as\n\$p = o / (1 + o)\$. Thus, if someone says the odds are 1 to 1,\nthey’re saying that the odds are 1 and thus \$p = 1 / (1 + 1) = 0.5\$.\nConversely, if someone says that the probability that something occurs\nis 50%, then they mean that \$p=0.5\$ so that the odds are\n\$o = p / (1 - p) = 0.5 / (1 - 0.5) = 1\$.

\n\n

The odds are famously derived using a standard fair game setting.\nImagine that you are playing a game where you flip a coin with success probability\n\$p\$ with the following rules:

\n\n

If it comes up heads, you win \$X\$ dollars.
If it comes up tails, you lose \$Y\$.

\n\n

What should we set \$X\$ and \$Y\$ for the game to be fair?\nFair means the expected earnings for either player is 0. That is:

\n\n\n$$\nE[earnings]= X p - Y (1 - p) = 0\n$$\n\n

This implies \$\\frac{Y}{X} = \\frac{p}{1 - p} = o\$.\nConsider setting \$X=1\$, then \$ Y = o\$.\nThus, the odds can be interpreted as\n“How much should you be willing to pay for a \$p\$ probability\nof winning a dollar?”\nIf \$p > 0.5\$ you have to pay more if you lose than you get if you win.\nIf \$p < 0.5\$ you have to pay less if you lose than you get if you win.

\n\n

So, imagine that I go to a horse race and the odds that a horse loses are\n50 to 1. They usually specify in terms of losing at horse tracks, so this\nwould be said to be 50 to 1 “against” where the against is implied and\nnot stated on the boards. The odds of the horse winning are\nthen 1/50. Thus, for a fair bet if I were to bet on the horse winning,\nthey should pay me 50 dollars if he wins and should pay 1 dollar if he\nloses. (Or any agreed upon multiple, such as\n100 dollars if he wins and 2 dollars if he loses.) The implied probability\nthat the horse loses is \$50 / (1 + 50)\$.

\n\n

It’s an interesting side note that the house sets the odds (hence the implied probability)\nonly by the bets coming in. They take a small fee for every bet win or lose\n(the rake). So, by setting the odds dynamically as the bets roll in,\nthey can guarantee that they make money (risk free) via the rake. Thus\nthe phrase “the house always wins” applies literally. Even more interesting\nis that by the wisdom of the crowd, the final probabilities tend to match the percentage\nof times that event happens. That is, house declared 50 to 1 horses\ntend to win about 1 out of 51 times and lose 50 out of 51 times,\neven though that 50 to 1 was set by a random collection\nof mostly uninformed bettors. This is why even your sports\njunkie friend with seemingly endless up to date sports knowledge\ncan’t make a killing betting on sports; the crowd is just too smart as a group\neven though most of the individuals know much less.

\n\n

Finally, and then I’ll get back on track, many of the machine learning\nalgorithms work on this principle of the wisdom of crowds:\nmany small dumb models can make really\ngood models. This is often called ensemble learning, where a lot of independent\nweak classifiers are combined to make a strong one. Random forests and boosting\ncome to mind as examples. In the Coursera Practical Machine Learning class,\nwe cover ensemble learning algorithms.

\n\n

Getting back to the issues at hand, recall that\nprobabilities are between 0 and 1. Odds are between 0 and infinity. So, it\nmakes sense to model the log of the odds (the logit), since it goes from\nminus infinity to plus infinity. The log of the odds is called the\nlogit:

\n\n\n$$\ng = \\mathrm{logit}(p) = \\log(p / (1 - p))\n$$\n\n

We can go backwards from the logit to the probability with the so-called\nexpit (inverse logit):

\n\n\n$$\n\\mathrm{expit}(g) = e^g / (1 + e^g) = 1 / (1 + e^{-g}) = p.\n$$\n\n

\nModeling the odds

\n\n

Let’s come up with notation for modeling the odds. Recall that\n\$Y_i\$ was our outcome, a 0 or 1 indicator of whether or not\nthe Raven’s won game \$i\$.

\n\n

Let

\n\n\n$$\np_i = \\mathrm{P}(Y_i = 1 ~|~ X_i = x_i, \\beta_0, \\beta_1)\n$$\n\n

be the probability of a win for number of points \$x_i\$.\nThen the odds that they win is \$p_i / (1 - p_i)\$ and the\nlog odds is \$\\log{p_i / (1 - p_i)} = \\mathrm{logit}(p_i)\$.

\n\n

Logistic regression puts the model on the log of the odds (logit)\nscale. In other words,

\n\n\n$$\n\\mathrm{logit}(p_i) = \\beta_0 + \\beta_1 x_i\n$$\n\n

Or, equivalently, we could just say

\n\n\n$$\nP(Y_i = 1 ~|~ X_i = x_i, \\beta_0, \\beta_1)\n= p_i = \\frac{ \\exp(\\beta_0 + \\beta_1 x_i)}{1 + \\exp(\\beta_0 + \\beta_1 x_i)}\n$$\n\n

\nInterpreting Logistic Regression

\n\n

Recall our model is:

\n\n\n$$\n\\mathrm{logit}\\{P(Y_i = 1 ~|~ X_i = x_i)\\} = \\beta_0 + \\beta_1 x_i\n$$\n\n

Let’s write this as:

\n\n\n$$\n\\mathrm{log}\\{O(Y_i = 1 ~|~ X_i = x_i)\\} = \\beta_0 + \\beta_1 x_i\n$$\n\n

where \$O(Y_i = 1 ~|~ X_i = x_i)\$ refers to the odds.\nInterpreting \$\\beta_0\$ is straightforward, it’s the\nlog of the odds of the Ravens winning for a 0 point game. Just\nlike in regular regression, for this to have meaning, a 0\nX value has to have meaning. In this case, there’s a structural\nconsideration that’s being ignored in that the Ravens can’t\nwin if they score 0 points (they can only tie or lose). \nThis is an unfortunate assumption of our model.

\n\n

For interpreting the \$\\beta_1\$ coefficient, consider\nthe following:

\n\n\n$$\n\\mathrm{log}\\{O(Y_i = 1 ~|~ X_i = x_i + 1)\\} -\n\\mathrm{log}\\{O(Y_i = 1 ~|~ X_i = x_i)\\}\n=\n\\mathrm{log}\\left\\{\n \\frac{O(Y_i = 1 ~|~ X_i = x_i + 1)}\n{O(Y_i = 1 ~|~ X_i = x_i)} \\right\\}\n= \\beta_1\n$$\n\n

So that \$\\beta_1\$ is the log of the relative increase\nin the odds of the Ravens winning\nfor a one point increase in score. The ratio of two odds is called,\nnot surprisingly, the odds ratio. So \$\\beta_1\$ is the log\nodds ratio of the Ravens winning associated with a one point increase\nin score.

\n\n

We can get rid of the log by exponentiating and then\nget that \$\\exp(\\beta_1)\$ is the odds ratio associated with\na one point increase in score. It’s a nifty fact that you\ncan often perform this exponentiation in your head, since for\nnumbers close to zero, exponentiation is about 1 + that number.\nSo, if you have a logistic regression slope\ncoefficient of 0.01, you know that e to that\ncoefficient is about 1.01. So you know that the coefficient estimates\na 1% increase in the odds of a success for every 1 unit increase\nin the regressor.

\n\n

\nVisualizing fitting logistic regression curves

\n\n

\n\n

Let’s visualize what the logistic regression model is
\nfitting. Consider setting \$\\beta_0\$ to 0 and\nvarying \$\\beta_1\$ for X being a regressor\nequally spaced between -10 and 10. Notice that\nthe logistic curves vary in their curvature.

\n\n

x = seq(-10, 10, length = 1000)\nbeta0 = 0; beta1s = seq(.25, 1.5, by = .1)\nplot(c(-10, 10), c(0, 1), type = \"n\", xlab = \"X\", ylab = \"Probability\", frame = FALS\\\nE)\nsapply(beta1s, function(beta1) {\n    y = 1 / (1 + exp( -1 * ( beta0 + beta1 * x ) ))\n    lines(x, y, type = \"l\", lwd = 3)\n}\n)\n

\n\n\n\n\n

\n\n\n

Try making the slope negative and see what happens. (It flips the curve\nfrom increasing to decreasing.)\nNow let’s hold \$\\beta_1\$ fixed and vary \$\\beta_0\$.

\n\n

x = seq(-10, 10, length = 1000)\nbeta0s = seq(-2, 2, by = .5); beta1 = 1\nplot(c(-10, 10), c(0, 1), type = \"n\", xlab = \"X\", ylab = \"Probability\", frame = FALS\\\nE)\nsapply(beta0s, function(beta0) {\n    y = 1 / (1 + exp( -1 * ( beta0 + beta1 * x ) ))\n    lines(x, y, type = \"l\", lwd = 3)\n}\n)\n

\n\n\n\n\n

\n\n\n

Notice that varying the intercept shifts the curve back and forth.\nLet’s superimpose some data with the fitted curve.

\n\n

x = seq(-10, 10, length = 1000)\nbeta0 = 0; beta1 = 1\np = 1 / (1 + exp(-1 * (beta0 + beta1 * x)))\ny = rbinom(prob = p, size = 1, n = length(p))\n\nplot(x, y, frame = FALSE, xlab = \"x\", ylab = \"y\")\nlines(lowess(x, y), type = \"l\", col = \"blue\", lwd = 3)\nfit = glm(y ~ x, family = binomial)\nlines(x, predict(fit, type = \"response\"), lwd = 3, col = \"red\")\n

\n\n\n\n\n

\n\n\n

The plot above shows the simulated binary data (black points), the\nfitted logistic curve (red) and a lowess smoother through the data (blue).\nThe lowess smoother shows a non-parametric estimate of the probability of a success\nat each x value. Think of it as a moving proportion. Logistic regression\ngets to move around the intercept and slope of the logistic curve to fit\nthe data well. Here the fit says that the probability of a 1 for low values of\nx is very small, the probability of a 1 for high values of x is high and\nit is intermediate at the points in the middle.

\n\n

\nRavens logistic regression

\n\n

\n\n

Now let’s run our binary regression model on the Ravens data.

\n\n

> logRegRavens = glm(ravensData$ravenWinNum ~ ravensData$ravenScore,family=\"binomial\\\n\")\n> summary(logRegRavens)\nCall:\nglm(formula = ravensData$ravenWinNum ~ ravensData$ravenScore,\n    family = \"binomial\")\n\nDeviance Residuals:\n   Min      1Q  Median      3Q     Max  \n-1.758  -1.100   0.530   0.806   1.495  \n\nCoefficients:\n                      Estimate Std. Error z value Pr(>|z|)\n(Intercept)            -1.6800     1.5541   -1.08     0.28\nravensData$ravenScore   0.1066     0.0667    1.60     0.11\n\n(Dispersion parameter for binomial family taken to be 1)\n\n    Null deviance: 24.435  on 19  degrees of freedom\nResidual deviance: 20.895  on 18  degrees of freedom\nAIC: 24.89\n\nNumber of Fisher Scoring iterations: 5\n\n## plotting the fit\n> plot(ravensData$ravenScore,logRegRavens$fitted,pch=19,col=\"blue\",xlab=\"Score\",ylab\\\n=\"Prob Ravens Win\")\n

\n\n\n\n\n

\"Fitted — Fitted model for the Ravens data.

\n\n\n

In this case, the data only covers some of the logistic curve, so that the\nfull “S” of the curve isn’t visible. To interpret our coefficients,\nlet’s exponentiate them.

\n\n

> exp(logRegRavens$coeff)\n          (Intercept) ravensData$ravenScore\n               0.1864                1.1125\n> exp(confint(logRegRavens))\n                         2.5 % 97.5 %\n(Intercept)           0.005675  3.106\nravensData$ravenScore 0.996230  1.303\n

\n\n\n\n

The first line of code shows that the exponentiated slope\ncoefficient is 1.11. Thus, we estimate a 11% increase\nin the odds of winning per 1 point increase in score. However,\nthe data are variable and the confident interval goes\nfrom 0.99 to 1.303. Since this interval contains 1 (or\ncontains 0 on the log scale), it’s not statistically significant.\n(It’s pretty close, though.)

\n\n

If we had included another variable in our model, say\nhome versus away game indicator, then our slope is\ninterpreted holding the value of the covariate held\nfixed. Just like in multivariable regression.

\n\n

We can also compare nested models using ANOVA and,\nby and large, our general model discussion carries\nover to this setting as well.

\n\n

\nSome summarizing comments

\n\n

Odds aren’t probabilities. In binary GLMs, we model the\nlog of the odds (logit) and our slope parameters are\ninterpreted as log odds ratios. Odds ratios of 1 or\nlog odds ratios of 0 are interpreted as no effect\nof the regressor on the outcome.

\n\n

\nExercises

\n\n

Load the dataset Seatbelts as part of the datasets package via data(Seatbelts). Use\nas.data.frame to convert the object to a dataframe. Create a new outcome variable\nfor whether or not greater than 119 drivers were killed that month. Fit a logistic regression GLM with this variable\nas the outcome and kms, PetrolPrice and law as predictors. Interpret your parameters. Watch a video solution.\n
Fit a binomial model with DriversKilled as the outcome and drivers as the total count\nwith kms , PetrolPrice and law as predictors, interpret your results. Watch a video solution.\n
Refer to Question 1. Use the anova function to compare models with just law, law and PetrolPrice and\nall three predictors. Watch a video solution.\n

\n\n\n

\nCount data

\n\n

Acknowledgment to Jeff Leek for much of the code and organization of this chapter.

\n\n

Many data take the form of unbounded count\ndata. For example, consider the number of calls\nto a call center or the number of flu cases\nin an area or the number of hits to a web site.

\n\n

In some of these cases the counts are clearly\nbounded. However, modeling the counts as unbounded\nis often done when the upper limit is not known\nor is very large relative to the number of events.

\n\n

If the upper bound is known, the techniques we’re\ndiscussing can be used to model the proportion or\nrate. The starting point for most count analysis is the\nthe Poisson distribution.

\n\n

\nPoisson distribution

\n\n

The Poisson distribution is the goto distribution for modeling\ncounts and rates. We’ll define a rate as a count per unit of time.\nFor example your heart rate is often expressed in beats per minute.\nSo, we might look at web hits per day, or disease cases per\nexposure time (incidence rates). Also, though not exactly a rate,\nwe can treat proportions as if rates when \$n\$ is large\nand the success probability is small.

\n\n

We would write that a random variable is Poisson,\n\$X \\sim \\mathrm{Poisson}(t\\lambda)\$, if its density function is:

\n\n\n$$\nP(X = x) = \\frac{(t\\lambda)^x e^{-t\\lambda}}{x!}\n$$\n\n

where \$x = 0, 1, \\ldots\$. The\nmean of the Poisson is \$E[X] = t\\lambda\$, thus \$E[X / t] = \\lambda\$.\nThe variance of the Poisson is \$Var(X) = t\\lambda\$.\nThe Poisson tends to a normal as \$t\\lambda\$ gets large and\napproximates a binomial with large \$n\$ and small \$p\$\nwhere we would think of \$t\\lambda\$ as \$n p \$.

\n\n

Here are some plots of the Poisson density to illustrate\nhow it closely approximates a normal.

\n\n

par(mfrow = c(1, 3))\nplot(0 : 10, dpois(0 : 10, lambda = 2), type = \"h\", frame = FALSE)\nplot(0 : 20, dpois(0 : 20, lambda = 10), type = \"h\", frame = FALSE)\nplot(0 : 200, dpois(0 : 200, lambda = 100), type = \"h\", frame = FALSE)\n

\n\n\n\n\n

\n\n\n

\nPoisson distribution

\n\n

Let’s analyze some data using the Poisson distribution.\nConsider the daily counts to Jeff Leek’s web site: http://biostat.jhsph.edu/~jleek/

\n\n

We’re interested in the number of hits per day.\nSince the unit of time is always one day, set \$t = 1\$ and then\nthe Poisson mean is interpreted as web hits per day.\nIf we set \$t = 24\$ then our Poisson rate\nwould be interpreted as web hits per hour.

\n\n

Let’s load the data:

\n\n

> download.file(\"https://dl.dropboxusercontent.com/u/7710864/data/gaData.rda\",destfi\\\nle=\"./data/gaData.rda\",method=\"curl\")\n> load(\"./data/gaData.rda\")\n> gaData$julian = julian(gaData$date)\n> head(gaData)\n        date visits simplystats julian\n1 2011-01-01      0           0  14975\n2 2011-01-02      0           0  14976\n3 2011-01-03      0           0  14977\n4 2011-01-04      0           0  14978\n5 2011-01-05      0           0  14979\n6 2011-01-06      0           0  14980\n\n> plot(gaData$julian,gaData$visits,pch=19,col=\"darkgrey\",xlab=\"Julian\",ylab=\"Visits\")\n

\n\n\n\n\n

\n\n\n

\nLinear regression

\n\n

We could try to fit the data with linear regression.\nThis is an often reasonable thing to do. Specifically,\nwe would start with the model

\n\n\n$$\nY_i = \\beta_0 + beta_1 x_i + \\epsilon_i\n$$\n\n

where \$Y_i\$ is the number of web hits on day i\nand \$x_i\$ is the day (expressed as a Julian\ndate, the number of days since January 1st, 1970).

\n\n

This model isn’t anywhere near as objectionable as when applied in\nthe binary case. Two sticking points remain. First, the\nresponse is a count and thus is non-negative, while our model\ndoesn’t acknowledge that. Secondly, the errors are typically\nassumed Gaussian, which is not an accurate approximation for\nsmall counts. As the counts get larger, straight application of\nlinear or multivariable regression models becomes more compelling.

\n\n

1 > plot(gaData$julian,gaData$visits,pch=19,col=\"darkgrey\",xlab=\"Julian\",ylab=\"Visits\")\n2 > lm1 = lm(gaData$visits ~ gaData$julian)\n3 > abline(lm1,col=\"red\",lwd=3)\n

\n\n\n\n\n

\n\n\n

The non-negativity could be handled by a (natural)\nlog transformation of\nthe outcome. The log has a special interpretation with respect\nto regression, so we cover it here. First, there is the issue of\nzero counts (which can’t be logged). Often a +1 is added to\nall data before taking the log, a reasonable solution but\none that harms the nice interpretation properties of the log.\nSecondly, a square root or cube root\ntransformation is often applied (which works just fine on zeros).\nWhile correcting nicely for skewness, this approach creates the issue\nof losing the nice interpretation of logs. Thus for the time being,\nlet’s assume no zero counts in the discussion.

\n\n

Consider now the model:

\n\n\n$$\n\\log(Y_i) = \\beta_0 + \\beta_1 x_i + \\epsilon_i\n$$\n\n

The quantity \$e^{E[\\log(Y)]}\$ estimates the geometric mean of \$Y\$.\nWhen you take the natural log of outcomes and fit a regression model, your exponentiated coefficients estimate things about geometric means.\nThus \$e^{\\beta_0}\$ is the geometric mean of hits on\nday 0, while \$e^{\\beta_1}\$ is the relative increase or\ndecrease in hits going from one day to the next.

\n\n

Let’s try this briefly with Jeff’s data.

\n\n

> round(exp(coef(lm(I(log(gaData$visits + 1)) ~ gaData$julian))), 5)\n  (Intercept) gaData$julian\n        0.000         1.002\n

\n\n\n\n

Thus our model estimates a 0.2% increase in geometric mean\ndaily web hits each day. What’s nice about the geometric mean is it’s\na multiplicative quantity. In this case it make sense to think\nmultiplicatively, as we would very naturally think in the terms\nof percent increases or decreases in the daily rate of web traffic.

\n\n

\nPoisson regression

\n\n

Notes