Statistical tests vs. linear regression

One-sample t-test vs. linear regression

Since statistical tests family, from a simple t-test to the MANCOVA are all wrapped up around linear model concept in R, linear regression are able to reproduce the results of even the simplest statistical tests.

Run the first two lines of code below and see the output of both for yourself, if you wish. The rest of the code simply packs the results in a nice looking table: (1) tidy(object, conf.int = T) from broom package was designed to produce nice model outputs, (2) kable from knitr package makes the table look publication like, (3) %>% and select are data-manipulation functions from tidyverse package-family and (4) rbind simply puts rows beneath each other, while cbind puts the columns near each other. Since select is used by several packages, we need to tell R to use select from dplyr:: package.

# statistical test (st)
st = t.test(mtcars$mpg)

# equivalent intercept-only linear regression (lr) 
lr = lm(mtcars$mpg ~ 1)

# nice looking table
kable(
  rbind(
    tidy(st)[c(1:3, 5:7)], 
    cbind(
      tidy(lr, conf.int = T)[c(2,4:7)], method = "Simple Linear Model") ) %>% 
  dplyr::select(method, everything()) )

You see: the results are identical!

Paired sample t-test vs. linear regression

To make the comparison less abstract, let’s make up (simulate) and visualise two samples.

# stabilize randomness for reproducibility
set.seed(1)

# tibble is simply a new table
data <- tibble(
  a = rnorm(n = 10, mean = 0,   sd = 9),
  b = rnorm(n = 10, mean = 5,   sd = 7) )

# ggplot is also included into tidyverse package-family
ggplot(data %>% gather(), aes(x = factor(key), y = value))+
  geom_boxplot()+
  geom_dotplot(binaxis ='y', stackdir = 'center', binwidth = 1, alpha = .3)+
  xlab("")+
  theme_minimal()

Since, paired t-test calculates the difference between samples and then compares this difference to zero (making it a one-sample t-test, where one-sample is our difference), in the regression model we simply subtract one sample from the other and perform intercept only model.

st = t.test(data$a, data$b, paired = TRUE)
lr = lm(data$a - data$b ~ 1)

kable(
  rbind(
    tidy(st)[c(1:3, 5:7)], 
    cbind(
      tidy(lr, conf.int = T)[c(2,4:7)], method = "Simple Linear Model")) %>% 
  dplyr::select(method, everything()) )  

Also here are all results identical. Thus, mean of the differences in paired t.test() is indeed the same as the (Intercept) in the linear regression.

Unpaired sample t-test vs. linear regression

Since statistical tests compare means of groups (samples), let’s reshape our data so that groups are in one column and values of our samples are in the second column. gather() function from tidyr package makes this easy to do. Reshaping the data will help us to understand tilde operator ~, which is often used in tests like ANOVA and always used in the linear models to define a model formula. The possibility of using ~ within statistical tests also indicates the affinity of statistical tests with models. This affinity is the reason you often see “ANOVA model” in the literature instead of “ANOVA test”.

paired = FALSE in the code below is actually unnecessary, since t-test uses it per default, but it underlines our intention to use unpaired t-test here.

data2 <- data %>% gather() %>% rename(group = key, value = value)

st = t.test(value ~ group, data = data2, paired = FALSE, var.equal = TRUE)
lr = lm(value ~ group, data = data2)

kable( tidy(st)[c(1:4, 6:8)] )

kable( tidy(lr, conf.int = T) )

Here we can see that the mean of the first sample “a” mean(data$a) = 1.189825, which is the estimate1 in the t-test is identical to the (Intercept) of the linear model. The second estimate of the t-test (6.741915) is identical to the sum of estimates in the linear model (1.189825+5.552090 = 6.741915). This makes sense, because the (Intercept) is always the mean of the first sample, while the estimate for other groups are always a transition from the first (reference) sample to the particular one (e.g. groupb). See, same thing, displayed a bit differently. But why? Well, while t-test shows you two means, the transition from one mean to the other in linear model shows you the slope - rate of change of the second group compared to the first, which will be very useful when we start adding more variables to our models.

T-test indicates that groups are not significantly different (p = 0.104). The switching from sample “a” to sample “b” in the linear model (see the second row in the model output) has identical p-value = 0.104. The t-statistics and confident intervals in the second row are also identical, but with the opposite sign (which I can’t explain at the moment. So, please, leave a comment if you can). But what it the first row about? Interestingly, the first row in the model output presents the results of a one-sample t-test of the first sample “a”. Run the code below and compare estimate, t-statistic, p-values and confidence intervals with the first row of the model output above. You’ll see that the results almost identical.

t.test(data$a)

## 
##  One Sample t-test
## 
## data:  data$a
## t = 0.53557, df = 9, p-value = 0.6052
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -3.835753  6.215403
## sample estimates:
## mean of x 
##  1.189825

One-way ANOVA vs. linear regression

Let’s take the Wage dataset from ISLR package to quickly get 5 group in our categorical variable we want to compare. Then we compare results of two models, ANOVA and regression.

Similarly to tidy, glance from broom package also presents model results in a tidy way.

st = aov(wage ~ education, Wage)
lr = lm(wage ~ education, Wage)

kable( rbind(glance(st), glance(lr)) )

The results are, again, identical!

Two-way ANOVA vs. linear regression

Even if we add more predictors and interactions into our model, the results statistical tests (here represented by Two-way ANOVA) and regression will match. Even the syntaxes in R is the same. Syntaxes might differ in other programming languages, but since R is the primer statistical language and we here try to learn stats, we stick with R for a while.

st = aov(wage ~ education * jobclass, Wage)
lr = lm(wage ~ education * jobclass, Wage)

kable( rbind(glance(st), glance(lr)) )

One-way ANCOVA vs. linear regression

st = aov(wage ~ education * age, Wage)
lr = lm(wage ~ education * age, Wage)

kable( rbind(glance(st), glance(lr)) )

Any-way MANOVA vs. linear regression

The multivariate analysis of variance (MANOVA) can also be conducted by lm function:

summary.aov( aov(cbind(wage, age)~education, Wage) )

##  Response wage :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## education      4 1226364  306591  229.81 < 2.2e-16 ***
## Residuals   2995 3995721    1334                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response age :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## education      4   4608 1151.90  8.7353 5.227e-07 ***
## Residuals   2995 394941  131.87                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary.aov(  lm(cbind(wage, age)~education, Wage) )

##  Response wage :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## education      4 1226364  306591  229.81 < 2.2e-16 ***
## Residuals   2995 3995721    1334                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response age :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## education      4   4608 1151.90  8.7353 5.227e-07 ***
## Residuals   2995 394941  131.87                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Please, don’t be confused by the function summary.aov(object, ...) used for linear model. It simply computes tables for one or more fitted model objects, where object contains the results returned by a model fitting function (e.g., aov, lm or glm).

Any-way MANCOVA vs. linear regression

# library(jmv)
st <- mancova(data = Wage,
              deps = vars(wage, age),
              factors = education, covs = year)
lr = lm(cbind(wage,age) ~ education+year, Wage)

st$univar                   

## 
##  Univariate Tests                                                                                 
##  ──────────────────────────────────────────────────────────────────────────────────────────────── 
##                 Dependent Variable    Sum of Squares    df      Mean Square    F         p        
##  ──────────────────────────────────────────────────────────────────────────────────────────────── 
##    education    wage                         1226364       4         306591    230.70    < .001   
##                 age                             4608       4           1152      8.74    < .001   
##    year         wage                           16807       1          16807     12.65    < .001   
##                 age                              494       1            494      3.75     0.053   
##    Residuals    wage                         3978914    2994           1329                       
##                 age                           394446    2994            132                       
##  ────────────────────────────────────────────────────────────────────────────────────────────────

summary.aov(lr)

##  Response wage :
##               Df  Sum Sq Mean Sq F value    Pr(>F)    
## education      4 1226364  306591 230.700 < 2.2e-16 ***
## year           1   16807   16807  12.647 0.0003821 ***
## Residuals   2994 3978914    1329                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response age :
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## education      4   4608 1151.90  8.7433 5.149e-07 ***
## year           1    494  494.13  3.7507   0.05288 .  
## Residuals   2994 394446  131.75                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1