Linear regression on continuous data

There are various ways to draw \(H_0\), but let’s focus on what it means in terms of the slope of a line. If slope is a measure of effect and there is zero effect, we’d expect zero slope, as well, or \(\beta_1 = 0\):

We can use lm() to estimate the two parts of this line, the intercept \(\beta_0\), and the slope, \(\beta_0\):

m1 <- lm(Y ~ X, line_data)
coef(m1)

##  (Intercept)            X 
## 16.648587544  0.003142142

As is clear from the graph, the line crosses the Y axis at 16.6 and has a slope dang near 0. We have no evidence on which to reject the null hypothesis.

Now consider a case in which the value of Y quite clearly increases with the value of X:

m2 <- lm(Y ~ X, line_data)
coef(m2)

## (Intercept)           X 
##  -47.038254    5.089453

This line crosses the Y axis at -47 and for every 1 unit X increases, Y increases 5 units. Hence, as X goes from 15 to 20, Y increases \(5~\times~5\) units, from 30 to 55. This is likely enough evidence to reject the null hypothesis, and the \(P\) value confirms it:

summary(m2)$coefficients

Both the intercept and slope are significantly different from 0, based on the \(t\) statistic for each term, with \(P\) values derived from the \(t\) distribution based on the model’s degrees of freedom. As with even a basic \(t\)-test, the \(t\) statistic in the summary table of the linear model is calculated as the ratio of the estimated effect (the coefficient) and the error around that estimate:

# Intercept (Beta 0)
  -47.04 / 11.38 

## [1] -4.133568

# X (Beta 1)
  5.089 / 0.6301

## [1] 8.076496

Comparing two groups

It becomes easy, then, to imagine finding the group means in the categorical data (let’s start with just two groups):

We know the old saying, two points make a line, and here we have two points… so let’s connect them with a line:

From here it is easy to see how a slope can be calculated from the difference between the means.

But slope (\(\beta_1\)) is just one of the two terms required to statistically describe a line. We also need an intercept (\(\beta_0\)).

This is the genius of ANOVA – how it gets an intercept from two categorical variables. It does so by simply assigning the first group to X = 0, making its mean by default the Y intercept. The next group is tested as X = 1, so the difference in the means is by default the slope, because when \(\Delta X = 1\),

\(\frac{\Delta Y}{\Delta X} = \frac{\Delta Y}{1} = \Delta Y\):

m3 <- lm(Y ~ Group, group_data)
coef(m3)

Indeed, \(\beta_0 = 13.85\), and \(\beta_1 = -7.07\), because \(13.85 - 6.78 = 7.07\)

summary(m3)$coefficients %>%
  pander::pander() 

Again, the model tells us whether the difference between the means is significant relative to the variance by calculating a \(t\) statistic and finding its P value from the t distribution based on the degrees of freedom.

To round out the example, let’s consider a case in which there is no difference between the means:

m4 <- lm(Y ~ Group, filter(group_data, Group !="B"))
summary(m4)$coefficients %>%
pander::pander() 

The difference in the means, \(1.2\), isn’t substantially greater than the variance, \(0.7\), so the resulting \(t\) statistic, \(1.7\), is too small to return a \(P\) value below the generally-accepted \(P = 0.05\). Thus the model confirms what we can see from the graph: the slope is too shallow, meaning there’s no meaningful difference between the group means.

t-statistics vs. F statistics

So far we’ve focused on the \(t\) statistic, which is returned by calls to summary() on lm objects. \(t\) statistics are useful because the \(t\) distribution is two-sided and symetrical, like a normal distribution. Thus, \(t\) values can either be negative or positive: the sign corresponds to the direction of the effect (increasing or decreasing), and the absolute value measures the size of the effect.

An alternative to the \(t\) statistic is the \(F\) statistic, which is based on a one-sided \(F\) distribution and is returned by an \(F\) test. The results of the \(F\) test are returned at the end of the results of summary():

summary(m3)

## 
## Call:
## lm(formula = Y ~ Group, data = filter(group_data, Group != "C"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1208 -0.8761 -0.2030  0.4395  4.1186 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.8507     0.4625   29.95  < 2e-16 ***
## GroupB       -7.0703     0.6540  -10.81 2.66e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.462 on 18 degrees of freedom
## Multiple R-squared:  0.8665, Adjusted R-squared:  0.8591 
## F-statistic: 116.9 on 1 and 18 DF,  p-value: 2.658e-09

We can also perform the \(F\) test alone – and get a full ANOVA table – by calling anova() on the lm object:

anova(m3) 

Yet again, the \(F\) statistic is obtained by comparing the term we hypothesize explains variation in the response variable (Group) against the error that remains unexplained (Residuals):

\(F = \frac{Mean~Sum~Squares_{Group}}{Mean~Sum~Squares_{Residuals}}\)

\(116.9 = \frac{249.9}{2.139}\)

The \(F\) test indicates that this is a statistically-significant effect, but because the \(F\) distribution is one-sided, we do not know whether the effect was an increase or a decrease without seeing the data or knowing the group means.

Comparing three (or more) groups

The \(F\) test is unique because it conducts a “global” test of the independent cateogrical variable. Compare the results of tests of Groups B and C against A:

But unlike a basic pairwise \(t\)-test, we can also test all three groups together:

m_all <- lm(Y ~ Group, group_data)
anova(m_all)

This is confusing, because the model returns a significant P value even though we know that one combination of means are not different (A vs. C).

When there are three or more groups in an \(F\) test, we must interpret a significant result as telling us that at least one group mean is different from at least one other group mean. Every group mean could be different from every other group mean, or only a single combination is different, or any combination in between. We just don’t know from the basic output of the \(F\) test. This is a second limitation (in addition to the one-sided \(F\) statistic).

summary(m_all)

aov(m_all) %>% TukeyHSD()