Lab 5

Try it yourself…

1. Simulate how power changes across sample sizes for a two-sample t-test.
2. Plot your results (sample size on the x axis, power on the y axis).

3. Increase or decrease your effect size and re-run your simulation/plot. How do sample size and effect size affect power?

   experiment <- function(group_size) {
     # sample data of group size from each of two groups with different means
     # return p-value of t-test comparing them
   }
   
   power <- function(group_size) {
     # run experiment 100 times
     # return the percentage of experiments that have significant p-values
   }
   
   power_levels <- data_frame(group_size = 2:200) %>%
     mutate(sample_size = group_size * 2,
            power = map_dbl(group_size, power))
   
   ggplot(power_levels) # add plot structure

   experiment <- function(group_size) {
     sample1 <- rnorm(group_size, 20, 10)
     sample2 <- rnorm(group_size, 25, 10)
     p_value <- t.test(sample1, sample2)$p.value
     return(p_value)
   }
   
   power <- function(group_size) {
     p_values <- replicate(100, experiment(group_size))
     num_sig <- sum(p_values < 0.05) / 100
     return(num_sig)
   }
   
   power_levels <- data_frame(group_size = 2:200) %>%
     mutate(sample_size = group_size * 2,
            power = map_dbl(group_size, power))
   
   ggplot(power_levels, aes(x = sample_size, y = power)) +
     geom_hline(yintercept = 0.8, colour = "grey", linetype = "dotted") +
     geom_point() +
     geom_smooth()

   

Try it yourself…

1. Model the relationship between RTlexdec and WrittenSpokenFrequencyRatio, plot it, and summarize in words what the result means.

2. Now look at the relationship between RTnaming and WrittenSpokenFrequencyRatio. Do you notice anything different?

   ratio_lm <- lm(RTlexdec ~ WrittenSpokenFrequencyRatio, data = rts)
   ratio_tidy <- tidy(ratio_lm)
   ratio_aug <- augment(ratio_lm)
   summary(ratio_lm)

   ## 
   ## Call:
   ## lm(formula = RTlexdec ~ WrittenSpokenFrequencyRatio, data = rts)
   ## 
   ## Residuals:
   ##      Min       1Q   Median       3Q      Max 
   ## -0.34617 -0.12579  0.00009  0.10287  0.63242 
   ## 
   ## Coefficients:
   ##                             Estimate Std. Error  t value Pr(>|t|)    
   ## (Intercept)                 6.546463   0.002684 2438.935  < 2e-16 ***
   ## WrittenSpokenFrequencyRatio 0.005364   0.001992    2.693  0.00711 ** 
   ## ---
   ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   ## 
   ## Residual standard error: 0.1568 on 4566 degrees of freedom
   ## Multiple R-squared:  0.001586,   Adjusted R-squared:  0.001367 
   ## F-statistic: 7.251 on 1 and 4566 DF,  p-value: 0.00711

   ggplot(ratio_aug, aes(sample = .resid)) + stat_qq()

   

   ggplot(ratio_aug, aes(x = WrittenSpokenFrequencyRatio, y = RTlexdec)) +
     geom_point(color = "grey") +
     geom_abline(intercept = filter(ratio_tidy, term == "(Intercept)")$estimate,
                 slope = filter(ratio_tidy, term == "WrittenSpokenFrequencyRatio")$estimate,
                 colour = "blue") +
     geom_smooth(method = lm, colour = "red")

   

   naming_lm <- lm(RTnaming ~ WrittenSpokenFrequencyRatio, data = rts)
   naming_aug <- augment(naming_lm)
   summary(naming_lm)

   ## 
   ## Call:
   ## lm(formula = RTnaming ~ WrittenSpokenFrequencyRatio, data = rts)
   ## 
   ## Residuals:
   ##      Min       1Q   Median       3Q      Max 
   ## -0.29190 -0.17437  0.01484  0.16785  0.37799 
   ## 
   ## Coefficients:
   ##                             Estimate Std. Error  t value Pr(>|t|)    
   ## (Intercept)                 6.318707   0.003053 2069.426   <2e-16 ***
   ## WrittenSpokenFrequencyRatio 0.005606   0.002266    2.474   0.0134 *  
   ## ---
   ## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   ## 
   ## Residual standard error: 0.1784 on 4566 degrees of freedom
   ## Multiple R-squared:  0.001339,   Adjusted R-squared:  0.00112 
   ## F-statistic: 6.122 on 1 and 4566 DF,  p-value: 0.01339

   ggplot(naming_aug, aes(sample = .resid)) + stat_qq()

   

   ggplot(naming_aug, aes(x = WrittenSpokenFrequencyRatio, y = RTnaming)) +
     geom_point(color = "grey") +
     geom_smooth(method = lm, colour = "red")

   

Here we looked at one continuous predictor and a continuous outcome, but regression is much more general than that. So we’re going to look at a few other ways you can use it.

Categorical predictors

Let’s say we want to look at the effect of being in one of 2 (or more) groups on some response. For example, we can look at RTlexdec as a function of the age of the subjects

rts_age <- rts %>%
  group_by(AgeSubject) %>%
  summarise(mean_rt = mean(RTlexdec))

ggplot(rts, aes(x = AgeSubject, y = RTlexdec)) +
  geom_violin() +
  geom_point(aes(y = mean_rt), data = rts_age, size = 3)

age_lm <- lm(RTlexdec ~ AgeSubject, data = rts)
summary(age_lm)

## 
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = rts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25776 -0.08339 -0.01669  0.06921  0.52685 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      6.660958   0.002324 2866.44   <2e-16 ***
## AgeSubjectyoung -0.221721   0.003286  -67.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared:  0.4992, Adjusted R-squared:  0.4991 
## F-statistic:  4552 on 1 and 4566 DF,  p-value: < 2.2e-16

This looks very familiar but if we go back to our line equation, something is weird.

\[ \text{RTlexdec} = 6.66 - 0.22 \cdot \text{AgeSubject} + \mathcal{N}(0, \sigma)\]

How do you multiply AgeSubject, which is either the word “old” or the word “young”, by a number?

Turns out R automatically turns “old” and “young” into 0s and 1s.

But you can’t see this until you convert it into a factor…

rts$AgeSubject <- factor(rts$AgeSubject)
contrasts(rts$AgeSubject)

##       young
## old       0
## young     1

It is treating “old” as 0 and “young” as 1 and it just did this alphabetically by default. This is called dummy coding.

So how do we read the output of our model?

age_tidy <- tidy(age_lm)
age_tidy

## # A tibble: 2 x 5
##              term   estimate   std.error statistic p.value
##             <chr>      <dbl>       <dbl>     <dbl>   <dbl>
## 1     (Intercept)  6.6609581 0.002323772 2866.4417       0
## 2 AgeSubjectyoung -0.2217215 0.003286310  -67.4682       0

The Intercept corresponds to the y value when the dummy coded variable is 0, so the intercept is the average RTlexdec for old people.

rts_age

## # A tibble: 2 x 2
##   AgeSubject  mean_rt
##        <chr>    <dbl>
## 1        old 6.660958
## 2      young 6.439237

AgeSubjectyoung is the coefficient that you add when the dummy coded variable is 1. So the average RTlexdec for young people is -0.22 less than the intercept.

intercept <- filter(age_tidy, term == "(Intercept)")$estimate
slope <- filter(age_tidy, term == "AgeSubjectyoung")$estimate

old_code <- contrasts(rts$AgeSubject)["old", ]
young_code <- contrasts(rts$AgeSubject)["young", ]

# old RTs
intercept + slope * old_code

## [1] 6.660958

# young RTs
intercept + slope * young_code

## [1] 6.439237

Dummy coding vs. effects coding

Another option is to do what is sometimes called effects coding. Instead of the default dummy coding that R assigns, you can change the contrast codes to something that is more meaningful to you.

contrasts(rts$AgeSubject) <- c(-0.5, 0.5)
contrasts(rts$AgeSubject)

##       [,1]
## old   -0.5
## young  0.5

So now the indicator variable for AgeSubject is -0.5 for “old” and 0.5 for “young”.

How will this affect our regression?

age_lm_2 <- lm(RTlexdec ~ AgeSubject, data = rts)
age_tidy_2 <- tidy(age_lm_2)
summary(age_lm_2)

## 
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = rts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25776 -0.08339 -0.01669  0.06921  0.52685 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.550097   0.001643 3986.29   <2e-16 ***
## AgeSubject1 -0.221721   0.003286  -67.47   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared:  0.4992, Adjusted R-squared:  0.4991 
## F-statistic:  4552 on 1 and 4566 DF,  p-value: < 2.2e-16

As you can see the slope is exactly the same – it’s still telling us that going from old to young (1 unit change) means subtracting -0.22. But note that the intercept is different. Now the intercept is the grand mean (across all conditions).

What if we rewrite our formulas from earlier:

intercept <- filter(age_tidy_2, term == "(Intercept)")$estimate
slope <- filter(age_tidy_2, term == "AgeSubject1")$estimate

old_code <- contrasts(rts$AgeSubject)["old", ]
young_code <- contrasts(rts$AgeSubject)["young", ]

# old RTs
intercept + slope * old_code

##      old 
## 6.660958

# young RTs
intercept + slope * young_code

##    young 
## 6.439237

We still get the exact same thing. In this case, where there’s just one predictor, the choice of contrast coding is entirely a matter of preference. In multiple regression, where you have many predictors potentially interacting, the right contrast coding can make a huge difference for how you interpret what’s going on in your model, specifically whether you are looking at simple effects or main effects.

There are many different coding schemes for testing various hypotheses, more examples: https://stats.idre.ucla.edu/r/library/r-%20library-contrast-coding-systems-for-%20categorical-variables/

Multiple Predictors

Now remember when we looked at the histogram, we saw that the RTs aren’t quite normally distributed in that there were two peaks. And now we know that young and old subjects respond differently. Let’s incorporate that into our regression. First some plots:

ggplot(rts, aes(x = WrittenFrequency, y = RTlexdec, colour = AgeSubject)) +
  geom_point(alpha = 0.1)

We see there are two groups here. Remember we are evaluating our model by seeing how close the fitted values are to the real data points, and we can get a lot closer if we have different lines for different age.

ggplot(rts, aes(x = WrittenFrequency, y = RTlexdec, colour = AgeSubject)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "lm")

Let’s make a regression model with multiple predictors:

\[ \text{RTlexdec} = \beta_0 + \beta_1 \cdot \text{WrittenFrequency} + \beta_2 \cdot \text{AgeSubject} + \mathcal{N}(0, \sigma)\]

freq_age_lm <- lm(RTlexdec ~ WrittenFrequency + AgeSubject, data = rts)
summary(freq_age_lm)

## 
## Call:
## lm(formula = RTlexdec ~ WrittenFrequency + AgeSubject, data = rts)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34622 -0.06029 -0.00722  0.05178  0.44999 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       6.7359314  0.0037621 1790.48   <2e-16 ***
## WrittenFrequency -0.0370103  0.0007033  -52.62   <2e-16 ***
## AgeSubject1      -0.2217215  0.0025930  -85.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08763 on 4565 degrees of freedom
## Multiple R-squared:  0.6883, Adjusted R-squared:  0.6882 
## F-statistic:  5040 on 2 and 4565 DF,  p-value: < 2.2e-16

• So now we have the intercept, which is the mean of RTlexdec for WrittenFrequency = 0 and across all AgeSubject (because of the effects coding).

• The effect of WrittenFrequency (-0.037) is pretty much the same as before which makes sense because we were seeing the same effect of Written Frequency in both groups

• The effect of age is very similar in size to when it was the sole predictor in the model as well.

Note that both the effects of Age and Frequency are significant suggesting that they are both independently contributing to explaining the variance in RTlexdec.

Critically, the \(r^2\) is much larger, suggesting we’re doing a much better job of fitting the data than with the previous model. All that extra variance is now being explained by adding the Age factor.

Linear model: \[ y_i = \beta_0 + \beta_1 \cdot x_{i1} + \beta_2 \cdot x_{i2} + \ldots + \beta_k \cdot x_{ik} + \epsilon_i \]