Regression with categorical variables

Packages for this section

library(tidyverse)
library(broom)

Jumping rats, again

my_url <- "http://datafiles.ritsokiguess.site/jumping.txt"
rats <- read_delim(my_url," ")

The data (randomly chosen rows)

rats %>% slice_sample(n = 10)

ANOVA as regression

  • With a quantitative response and a categorical explanatory variable, this is some kind of ANOVA
  • but what if we run it as a regression?
rats.1 <- lm(density ~ group, data = rats)

The output

summary(rats.1)

Call:
lm(formula = density ~ group, data = rats)

Residuals:
   Min     1Q Median     3Q    Max 
 -47.1  -13.3   -3.3   12.5   51.9 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    601.100      6.826  88.064  < 2e-16 ***
groupHighjump   37.600      9.653   3.895 0.000584 ***
groupLowjump    11.400      9.653   1.181 0.247912    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.58 on 27 degrees of freedom
Multiple R-squared:  0.3714,    Adjusted R-squared:  0.3249 
F-statistic: 7.978 on 2 and 27 DF,  p-value: 0.001895

The ANOVA \(F\)-test

anova(rats.1)
  • but not Tukey:
TukeyHSD(rats.1)
Error in `UseMethod()`:
! no applicable method for 'TukeyHSD' applied to an object of class "lm"

So what about those “slopes”?

tidy(rats.1)
  • The first category alphabetically is Control. This is the “baseline” category.
  • The intercept 601.1 is the mean bone density for the control group.
  • Highjump term: high-jumping rats have mean bone density that is 37.6 higher than the control rats (ie., 638.7)
  • Lowjump term: low-jumping rats have mean that is 11.4 higher than control, ie., 612.5.

Group means

rats %>% 
  group_by(group) %>% 
  summarize(mean_density = mean(density))

confirming what we worked out from the slopes.

Testing a categorical variable

The summary (or tidy) output only compares each group with the baseline, so does not give an overall test for any differences among groups. To do that, use drop1:

drop1(rats.1, test = "F")

This says that the jumping groups do not all have the same mean bone density (the same thing that the ANOVA says: here it is the same test).

The crickets

  • Male crickets rub their wings together to produce a chirping sound.
  • Rate of chirping, called “pulse rate”, depends on species and possibly on temperature.
  • Sample of crickets of two species’ pulse rates measured; temperature also recorded.
  • Does pulse rate differ for species, especially when temperature accounted for?

The crickets data

Read the data:

my_url <- "http://datafiles.ritsokiguess.site/crickets2.csv"
crickets <- read_csv(my_url)
crickets %>% slice_sample(n = 10) # display sample of rows

Scatterplot

ggplot(crickets, aes(x = temperature, y = pulse_rate, 
                     colour = species)) + 
  geom_point() + geom_smooth(method = "lm")

Comments

  • Relationship between temperature and pulse rate looks linear for each species, so add regression lines to the graph.
  • As temperature increases, pulse rate also increases (for both species)
  • At the same temperature, pulse rate higher for exclamationis than for niveus
  • The fit of the observations to their line looks very good.
  • The two lines are close to parallel: the rate of increase of pulse rate with temperature is almost the same for both species.

Fit model with lm

crickets.1 <- lm(pulse_rate ~ temperature + species, 
                 data = crickets)

Can I remove anything? No:

drop1(crickets.1, test = "F")

drop1 is right thing to use in a regression with categorical (explanatory) variables in it: “can I remove this categorical variable as a whole?”

The summary

summary(crickets.1)

Call:
lm(formula = pulse_rate ~ temperature + species, data = crickets)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.0128 -1.1296 -0.3912  0.9650  3.7800 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -7.21091    2.55094  -2.827  0.00858 ** 
temperature     3.60275    0.09729  37.032  < 2e-16 ***
speciesniveus -10.06529    0.73526 -13.689 6.27e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.786 on 28 degrees of freedom
Multiple R-squared:  0.9896,    Adjusted R-squared:  0.9888 
F-statistic:  1331 on 2 and 28 DF,  p-value: < 2.2e-16

Conclusions

  • Slope for temperature says that increasing temperature by 1 degree increases pulse rate by 3.6 (same for both species)
  • Species exclamationis is baseline.
  • Slope for speciesniveus says that pulse rate for niveus about 10 lower than that for exclamationis at same temperature
  • R-squared of almost 0.99 is very high, so that the prediction of pulse rate from species and temperature is very good.
  • Model assumes that the slope for each species is the same.
    • if not? See later.

Residuals

  • We can also check residuals from this kind of model:
crickets.1 %>% augment(crickets) -> crickets.1a

and then…

Residuals vs fitted

ggplot(crickets.1a, aes(x = .fitted, y = .resid)) + 
  geom_point()

Normal quantile plot of residuals

ggplot(crickets.1a, aes(sample = .resid)) + 
  stat_qq() + stat_qq_line()

Residuals vs. \(x\)’s

If you try to do them all at once, this happens:

crickets.1a %>% 
  pivot_longer(temperature:species, 
               names_to = "xname", values_to = "xval")
Error in `pivot_longer()`:
! Can't combine `temperature` <double> and `species` <character>.

so do them separately 1/2

ggplot(crickets.1a, aes(x = temperature, y = .resid)) + 
  geom_point()

so do them separately 2/2

ggplot(crickets.1a, aes(x = species, y = .resid)) + 
  geom_point()

Comments

  • Residuals vs fitted looks random
  • Normal quantile plot is reasonably close to normal
  • Residuals vs temperature looks random
  • Residuals vs species looks random
    • looking for residuals vs each species to be scattered around zero with about equal spread.

What if we think the slopes are not equal?

  • Add an interaction term species:temperature, shorthand as below:
crickets.2 <- lm(pulse_rate ~ species * temperature, 
                 data = crickets)
  • Test everything droppable for significance using drop1:
drop1(crickets.2, test = "F")
  • Interaction term not significant, so slopes are not significantly different.

Interpreting slopes in presence of interaction 1/2

tidy(crickets.2)

Interpreting slopes in presence of interaction 2/2

  • intercept is estimated pulse rate when temperature is 0 and species is baseline (exclamationis)
  • speciesniveus term is the change in intercept when species changes from baseline to niveus
  • temperature term is change in pulse rate when temperature increases by 1 but only for baseline species exclamationis
  • speciesniveus:temperature term is change in slope as temperature increases by 1 but for species niveus.

All of that as numbers

  • The line for exclamationis has intercept \(-11.0\) and slope \(3.75\)
  • The line for niveus has intercept \(-11.0 - 4.35 = -15.35\) and slope \(3.75 - 0.234 = 3.52\).
  • The two lines have almost the same slope.
  • So now we can do predictions. Let’s do temperature 24:
  • for exclamationis:
-11.0 + 3.75 * 24
[1] 79
  • and for niveus:
-15.35 + 3.52 * 24
[1] 69.13

Compare the scatterplot