Categorial Predictors

Getting ready:

library(tidyverse)
library(broom)

senses <- read_csv("winter_2016_senses_valence.csv")

Explore and describe the data

First some histograms:

senses %>% ggplot(aes(x = Val))+
  geom_histogram(binwidth = .25)+
  theme_bw()+
  facet_wrap(~Modality, nrow = 1)

And some dotplots:

senses %>% ggplot(aes(x = Modality, y = Val, shape = Modality, color = Modality)) +
  geom_point(alpha = .3)+
  stat_summary(fun.y = "mean", fun.ymin = "mean", fun.ymax= "mean", size= 0.3, geom = "crossbar") +
  theme_bw()

Finally, some summary stats by modality:

senses %>% group_by(Modality) %>% summarise(mean = mean(Val),
                                            sd = sd(Val),
                                            min = min(Val),
                                            max = max(Val))

Working with a 2-level categorical variable

Let’s create a smaller dataset with only Taste and Smell modalities:

chem <- filter(senses, Modality %in% c('Taste', 'Smell'))

table(chem$Modality)

Smell Taste 
   25    47 

Now for a regression model:

chem_mdl <- lm(Val ~ Modality, data = chem)

summary(chem_mdl)

Call:
lm(formula = Val ~ Modality, data = chem)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.99315 -0.20870  0.04343  0.19115  0.62788 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.47101    0.06297  86.889  < 2e-16 ***
ModalityTaste  0.33711    0.07793   4.326 4.95e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3148 on 70 degrees of freedom
Multiple R-squared:  0.2109,    Adjusted R-squared:  0.1997 
F-statistic: 18.71 on 1 and 70 DF,  p-value: 4.951e-05

We can infer that Smell is the reference category. Let’s play with a couple different ways to do treatment coding:

chem <- chem %>% mutate(Mod01 = if_else(Modality == "Taste", 1, 0),
                        Mod10 = if_else(Modality == "Smell", 1, 0),
                        Modality = factor(Modality),
                        ModRe = relevel(Modality, ref = 'Taste'))

And now for some models…

chem_mdl01 <- lm(Val ~ Mod01, data = chem)
chem_mdl10 <- lm(Val ~ Mod10, data = chem)
chem_mdlFac <- lm(Val ~ Modality, data = chem)
chem_mdlRe <- lm(Val ~ ModRe, data = chem)

Take a peek at the summaries and coefficient values…

#fit summary
glance(chem_mdl01)$r.squared

[1] 0.2109257

glance(chem_mdl10)$r.squared

[1] 0.2109257

glance(chem_mdlFac)$r.squared

[1] 0.2109257

glance(chem_mdlRe)$r.squared

[1] 0.2109257

#coefficients
tidy(chem_mdl01)$estimate

[1] 5.4710116 0.3371123

tidy(chem_mdl10)$estimate

[1]  5.8081239 -0.3371123

tidy(chem_mdlFac)$estimate

[1] 5.4710116 0.3371123

tidy(chem_mdlRe)$estimate

[1]  5.8081239 -0.3371123

Now let’s try Modality with sum contrasts:

chem <- mutate(chem, ModSum = Modality)

contrasts(chem$ModSum) <- contr.sum(2)

chem_mdlSum <- lm(Val ~ ModSum, data = chem)

summary(chem_mdlSum)

Call:
lm(formula = Val ~ ModSum, data = chem)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.99315 -0.20870  0.04343  0.19115  0.62788 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.63957    0.03897 144.729  < 2e-16 ***
ModSum1     -0.16856    0.03897  -4.326 4.95e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3148 on 70 degrees of freedom
Multiple R-squared:  0.2109,    Adjusted R-squared:  0.1997 
F-statistic: 18.71 on 1 and 70 DF,  p-value: 4.951e-05

Note that R^2, F-stat, etc. are the same even though intercept and slope differ.

t-test

t.test(Val ~ Modality, data = chem, var.equal = T)

    Two Sample t-test

data:  Val by Modality
t = -4.3257, df = 70, p-value = 4.951e-05
alternative hypothesis: true difference in means between group Smell and group Taste is not equal to 0
95 percent confidence interval:
 -0.4925442 -0.1816804
sample estimates:
mean in group Smell mean in group Taste 
           5.471012            5.808124 

Same p-value as regression model!

A predictor with more than 2 categories

Now we’ll take a look at the whole dataset, which has 5 categories for sense modalities.

sense_all <- lm(Val ~ Modality, data = senses)

summary(sense_all)

Call:
lm(formula = Val ~ Modality, data = senses)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.99315 -0.16482 -0.02158  0.15920  1.15734 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.57966    0.01889 295.308  < 2e-16 ***
ModalitySmell -0.10865    0.05643  -1.925   0.0549 .  
ModalitySound -0.17447    0.03758  -4.643 4.66e-06 ***
ModalityTaste  0.22846    0.04314   5.296 1.96e-07 ***
ModalityTouch -0.04523    0.03737  -1.210   0.2269    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2659 on 400 degrees of freedom
Multiple R-squared:  0.1455,    Adjusted R-squared:  0.137 
F-statistic: 17.03 on 4 and 400 DF,  p-value: 6.616e-13

To compare: ANOVA

Now for the ANOVA approach:

aov_senses <- aov(Val ~ Modality, data = senses)

summary(aov_senses)

             Df Sum Sq Mean Sq F value   Pr(>F)    
Modality      4  4.815  1.2036   17.03 6.62e-13 ***
Residuals   400 28.274  0.0707                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Same model fit - F-stat and p.val are identical.

getting at the contrasts…

The easystats package includes the modelbased package, which is helpful for extracting model contrasts of interest when working with categorical variables.

library(easystats)

We’ll first use the estimate_means() function to get the model-based estimates for each category. Thse are useful to report as descriptive statistics, or for plotting.

sense_all_means <- estimate_means(sense_all, by = "Modality")

sense_all_means

Estimated Marginal Means

Modality | Mean |   SE |       95% CI | t(400)
----------------------------------------------
Sight    | 5.58 | 0.02 | [5.54, 5.62] | 295.31
Smell    | 5.47 | 0.05 | [5.37, 5.58] | 102.89
Sound    | 5.41 | 0.03 | [5.34, 5.47] | 166.41
Taste    | 5.81 | 0.04 | [5.73, 5.88] | 149.77
Touch    | 5.53 | 0.03 | [5.47, 5.60] | 171.66

Variable predicted: Val
Predictors modulated: Modality

plot(sense_all_means)

We can use the estimate_contrasts() function to examine pairwise differences among categorial variables.

sense_all_contrasts <- estimate_contrasts(sense_all, contrast = "Modality",
                                          p_adjust = "tukey")

sense_all_contrasts

Marginal Contrasts Analysis

Level1 | Level2 | Difference |   SE |         95% CI | t(400) |      p
----------------------------------------------------------------------
Smell  | Sight  |      -0.11 | 0.06 | [-0.22,  0.00] |  -1.93 |  0.306
Sound  | Sight  |      -0.17 | 0.04 | [-0.25, -0.10] |  -4.64 | < .001
Taste  | Sight  |       0.23 | 0.04 | [ 0.14,  0.31] |   5.30 | < .001
Touch  | Sight  |      -0.05 | 0.04 | [-0.12,  0.03] |  -1.21 |  0.745
Sound  | Smell  |      -0.07 | 0.06 | [-0.19,  0.06] |  -1.06 |  0.829
Taste  | Smell  |       0.34 | 0.07 | [ 0.21,  0.47] |   5.12 | < .001
Touch  | Smell  |       0.06 | 0.06 | [-0.06,  0.19] |   1.02 |  0.846
Taste  | Sound  |       0.40 | 0.05 | [ 0.30,  0.50] |   7.97 | < .001
Touch  | Sound  |       0.13 | 0.05 | [ 0.04,  0.22] |   2.82 |  0.040
Touch  | Taste  |      -0.27 | 0.05 | [-0.37, -0.17] |  -5.43 | < .001

Variable predicted: Val
Predictors contrasted: Modality
p-value adjustment method: Tukey

You can do the same thing with the ANOVA model, too!