Worksheet 7

Published

February 22, 2026

Packages

library(tidyverse)
library(marginaleffects)
library(car)

Hair colour and pain tolerance

In lecture, we sped through a review example of one-way ANOVA, in which we investigated the effect of hair colour on pain tolerance. The data were in http://ritsokiguess.site/datafiles/hairpain.txt, with two columns hair (hair colour) and pain (pain tolerance, with a higher value indicating that the individual can withstand more pain). The hair colours were brown and blond, each subdivided into light and dark. The data values are separated by single spaces.

In this problem, we take a different approach to an analysis of the same data. In particular, we focus on three particular comparisons of pain tolerance:

  • light brown hair vs. dark brown hair
  • light blond hair vs. dark blond hair
  • the average of brown hair vs the average of blond hair.

These are comparisons we have chosen to make before looking at the data (because, we suppose, they are particular research questions of interest to us).

  1. Read in and display (some of) the data.

The last sentence of the first paragraph indicates that we do not have a .csv (for a change), but rather something that needs to be read in using read_delim:

my_url <- "http://ritsokiguess.site/datafiles/hairpain.txt"
hairpain <- read_delim(my_url, " ")
Rows: 19 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: " "
chr (1): hair
dbl (1): pain

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
hairpain
  1. For the analysis we are about to do, the categorical variable needs to be a factor. Create and save a dataframe for which this is the case.

The categorical variable is hair, so (overwriting the original dataframe, which appears to be safe):

hairpain %>% 
  mutate(hair = factor(hair)) -> hairpain
hairpain

The fctr (or fct) at the top of the new hair indicates that we have indeed made a factor column. (The contrasts stuff later needs a factor, hence making one now.)

  1. Find out how many observations you have for each hair type. (This has the side-effect of telling you which order the hair types are in, as far as R is concerned.)

This is most easily count:

hairpain %>% count(hair)

which reveals that the hair types are actually in alphabetical order. This is not exactly a surprise, but I figured it would make your life easier below to know exactly which hair type is first, second, etc. without you otherwise having to work it out.

  1. Set up contrasts to represent your comparisons of interest.

Something along these lines:

c_brown <- c(0, 1, 0, -1)
c_blond <- c(1, 0, -1, 0)
c_brown_blond <- c(-1/2, 1/2, -1/2, 1/2)

The names of the contrasts are chosen by you, but it is helpful to use names that say which things you are contrasting.

Pick out the two shades of brown (the second and fourth ones in question 3 with zeros elsewhere), the two shades of blond (the first and third ones in question 3 with zeros elsewhere), and the average of the two brown positive and of the two blond negative (for example; see below).

You have some alternatives:

  • you can switch all the positives and all the negatives in any contrast (for example, the first one could be c(0, -1, 0, 1)) because it is still “the same contrast”, that is, making the same comparison, this way).
  • you can also multiply or divide all the numbers in a contrast by the same thing (for example, the last one could be c(-1, 1, -1, 1)). This also makes “the same contrast”. My take when I am comparing averages is to divide by how many things the average is based on, in this case 2 (two types of brown hair and two types of blond hair). If, say, one of them had been two types and the other three, it would have been necessary to divide by the right thing, for example c(1/2, 1/2, -1/3, -1/3, -1/3), to make the coefficients of the contrast add up to zero (see below).

The four coefficients in each of your contrasts need to add up to zero (this is actually the definition of a contrast), which you can check as follows:

sum(c_brown)
[1] 0
sum(c_blond)
[1] 0
sum(c_brown_blond)
[1] 0

If they do, and the coefficients are of the right sign (relative to each other) and the right relative size, you are good.

  1. Verify that your contrasts are orthogonal.

Let’s recall our contrasts:

c_brown
[1]  0  1  0 -1
c_blond
[1]  1  0 -1  0
c_brown_blond
[1] -0.5  0.5 -0.5  0.5

To verify that two contrasts are orthogonal, say the second and third here, multiply the first coefficient in one by the first coefficient in the other, then repeat for the other coefficients, then add up the results, checking that your final result is zero. Here, that gives

\[ (1)(-0.5) + (0)(0.5) + (-1)(-0.5) + (0)(0.5) = (-0.5) + 0 + 0.5 + 0 = 0.\]

Check. Once you are happy that you understand the process, get R to do the heavy lifting the rest of the way:

sum(c_brown * c_blond)
[1] 0
sum(c_brown * c_brown_blond)
[1] 0
sum(c_blond * c_brown_blond)
[1] 0

All three pairs of contrasts are indeed orthogonal.

Having orthogonal contrasts means that the tests we are going to do below are independent of each other, so that the result of one does not influence any of the others.

  1. Set up your contrasts as a matrix, and set this matrix as the contrasts of your categorical variable. Then run the ANOVA as a regression.

That means this:

m <- cbind(c_brown, c_blond, c_brown_blond)
contrasts(hairpain$hair) <- m

cbind makes the matrix, and the thing inside contrasts should be your dataframe and categorical variable column separated by a $.

Then, run what would be an ANOVA as a regression using lm, and display the summary:

hairpain.1 <- lm(pain ~ hair, data = hairpain)
summary(hairpain.1)

Call:
lm(formula = pain ~ hair, data = hairpain)

Residuals:
   Min     1Q Median     3Q    Max 
-11.20  -5.45  -0.50   4.30  13.60 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         47.575      1.884  25.257 1.05e-13 ***
hairc_brown         -2.550      2.741  -0.930  0.36695    
hairc_blond         -4.000      2.584  -1.548  0.14251    
hairc_brown_blond  -15.250      3.767  -4.048  0.00105 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.172 on 15 degrees of freedom
Multiple R-squared:  0.576, Adjusted R-squared:  0.4912 
F-statistic: 6.791 on 3 and 15 DF,  p-value: 0.004114

Interpretation is coming shortly, but you might note that the three rows are no longer the intercept and slopes you would have expected, but are the three contrasts you set up.

Extra: There were four groups (hair types), so in a regular ANOVA there would be three df for hair types, and your \(F\)-test would have 3 and 15 df.1 When you analyze using contrasts, each contrast has 1 df and hence can be tested with a \(t\)-test, as here.2 As a result, you can have as many contrasts as you have (numerator) df, here 3 because there are 4 groups. You have some freedom to define the contrasts according to the things you want to test, but you can only have a fixed number of them, and they have to be orthogonal.3

  1. Interpret your results.

The three tests in the summary table are labelled with hair followed by the names you gave the contrasts. Ignore the intercept:

  • there is no difference in pain tolerance between the people with light brown and dark brown hair (P-value 0.37)
  • there is no difference in pain tolerance between the people with light blond and dark blond hair (P-value 0.14)
  • there is a difference in pain tolerance between brown-haired and blond-haired people (P-value 0.0011).

Note how this is more focused than the regular ANOVA, that looked at all the differences between groups, including some that we didn’t really care about.

There are only two hair colours, brown and blond, so we don’t need to do Tukey, but we can figure out which hair colour is associated with higher pain tolerance either by looking at group means:

hairpain %>% 
  group_by(hair) %>% 
  summarize(mean_pain = mean(pain))

or by reasoning it out from the contrast. In my contrast:

c_brown_blond
[1] -0.5  0.5 -0.5  0.5

the positive coefficients went with brown hair. The Estimate in the summary was negative, which means that the pain tolerance was higher in the people with the negative coefficients, namely the blond-haired people. This is consistent with my table of means just above.

Cats

144 cats, 47 female and 97 male, had been used in an experiment with the muscle-relaxing drug digitalis. For each cat, its Sex was recorded, along with the cat’s body weight (Bwt) in kilograms and the weight of its heart (Hwt) in grams. We are interested in the relationship between the weight of a cat’s heart (response) and its body weight (explanatory), and whether that relationship is different for male and female cats. The data are in the file http://ritsokiguess.site/datafiles/cats.csv.

  1. Read in and display (some of) the data.

As ever:

my_url <- "http://ritsokiguess.site/datafiles/cats.csv"
cats <- read_csv(my_url)
Rows: 144 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Sex
dbl (2): Bwt, Hwt

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
cats

144 rows, with the columns promised. The bit about digitalis is not actually at all relevant to what we are doing.

  1. Make a suitable graph of the three variables. Add appropriate regression lines to your graph.

Two quantitative and one categorical means a scatterplot of the two quantitative variables, bearing in mind which one of them is the response, with the categorical variable’s levels being indicated by colour. Here, we are predicting heart weight from body weight:

ggplot(cats, aes(x = Bwt, y = Hwt, colour = Sex)) + geom_point() + 
  geom_smooth(method = "lm")
`geom_smooth()` using formula = 'y ~ x'

Some things to note for yourself from the graph:

  • the male cats tended to be bigger than the females (the heaviest female cat had a body weight of 3 kg, but there were many male cats heavier than that)
  • cats with a bigger body weight tended to have a bigger heart, whether male or female
  • is there a difference between male and female cats? Not clear.
  • the slopes of the two regression lines look a bit different, but whether that difference is significant remains to be seen.
  1. Fit a suitable analysis of covariance model, and display its output.

Include an interaction (and check whether it is necessary, which is coming up):

cats.1 <- lm(Hwt ~ Bwt * Sex, data = cats)
summary(cats.1)

Call:
lm(formula = Hwt ~ Bwt * Sex, data = cats)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7728 -1.0118 -0.1196  0.9272  4.8646 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9813     1.8428   1.618 0.107960    
Bwt           2.6364     0.7759   3.398 0.000885 ***
SexM         -4.1654     2.0618  -2.020 0.045258 *  
Bwt:SexM      1.6763     0.8373   2.002 0.047225 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.442 on 140 degrees of freedom
Multiple R-squared:  0.6566,    Adjusted R-squared:  0.6493 
F-statistic: 89.24 on 3 and 140 DF,  p-value: < 2.2e-16

Extra:

  • The interaction is significant (just), so we need to keep it. That means that the two lines on your graph do in fact have significantly different slopes (which may or may not be what you guessed)
  • The Estimate for SexM says how the intercept for males compares to the baseline Sex, females. If you imagine a cat with body weight 0 (which is of course impossible) by extending the lines on your graph to the left, this says that the Male line will be considerably less than the Female line by the time you get back to a Bwt of zero. Unlike the example in lecture, the Male line ends up being above the Female one eventually (once Bwt gets above about 2.5), because it has a larger slope.
  • The Estimate for Bwt is an ordinary regression slope, but since the slopes of the two lines are different, we have to say that it is the slope for the baseline Sex, females. That is, a female cat that is 1 kg heavier will have a heart that is 2.64 grams heavier, on average.
  • The Intercept is likewise the intercept for females.
  • The interaction term represents the difference in slopes for male and female cats: the slope for male cats is 1.68 bigger than the slope for female cats.
  1. Is the interaction term significant? What does your answer mean in the context of the data?

The interaction term has a P-value of 0.047, which is just less than 0.05, so it is significant (at that \(\alpha\)). This means that the two lines predicting heart weight from body weight for males and females have different slopes.

This is unlike the first example we had in class, where the interaction term and hence the two slopes were not significantly different, but like the second one about ice cream. (In all honesty, the second lecture example was more interesting than this one.)

  1. For male and female cats of body weights 2.5 and 3.5 kg (all four combinations), obtain predicted heart weights.

This uses marginaleffects in the customary way. Make a new dataframe of values to predict for first, and then do your predictions, displaying only the columns of interest.

Hence, first this:

new <- datagrid(model = cats.1, Bwt = c(2.5, 3.5), Sex = c("F", "M"))
new

and then this:

cbind(predictions(cats.1, new)) %>% 
  select(estimate, Bwt, Sex)

Add things like confidence limits for the predictions if you like (though we’re not using them).

  1. (2 points) Using your predictions, verify that the slopes for males and females are different.

Take a trip back to high school and recall that slope is rise over run (or search for it and make a note of where you found out what to do). In this case, the “run” part is 1, because the two body weights are 1 kg apart, so the slope in each case is just the “rise” part:

Slope for (baseline) females is

12.21 - 9.57
[1] 2.64

and the slope for males is

13.91 - 9.60
[1] 4.31

These are indeed different, as required. (If you forgot to add an interaction term in your ANCOVA model, they will be the same.) You don’t need to be any more accurate than this.

Extra: going back to our output from earlier:

summary(cats.1)

Call:
lm(formula = Hwt ~ Bwt * Sex, data = cats)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7728 -1.0118 -0.1196  0.9272  4.8646 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   2.9813     1.8428   1.618 0.107960    
Bwt           2.6364     0.7759   3.398 0.000885 ***
SexM         -4.1654     2.0618  -2.020 0.045258 *  
Bwt:SexM      1.6763     0.8373   2.002 0.047225 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.442 on 140 degrees of freedom
Multiple R-squared:  0.6566,    Adjusted R-squared:  0.6493 
F-statistic: 89.24 on 3 and 140 DF,  p-value: < 2.2e-16

the slope for the baseline females is indeed the same as the Estimate for Bwt. The slope for males is the baseline slope, plus the change in slope between females and males, which is the interaction term, hence

2.64 + 1.68
[1] 4.32

as we found from the predictions.

Footnotes

  1. If you have done ANOVAs by hand, you would have looked up 3 and 15 df in your tables to get the P-value for your \(F\)-test, with a result like “the P-value is between 0.01 and 0.05”.↩︎

  2. An \(F\) with 1 df on the top is the square of a \(t\) with the same df on the bottom, so such tests are equivalent.↩︎

  3. Not strictly true, but it makes things a lot more complicated if you don’t have orthogonal contrasts, and that’s beyond our scope in this course. One of the PASIAS problems has some discussion on this if you want to learn more.↩︎