Week 5. Testing differences between groups

Written by Padraic Monaghan

5.1 Overview

This week, there are two mini lectures, and the practical workbook working with R-studio.

Before the practical on Tuesday, please try to work through the practical workbook in your group.

Bring your questions (and/or answers) to the practical.

5.2 Learning Goals

• Understand when t-tests are appropriately applied to data

• Understand the distinction between paired and independent t-tests

• Interpret p-values from t-tests

• Determine how a paired t-test is calculated

• Determine how an independent t-test is calculated

• Understand effect sizes for t-tests

• Be able to effectively interpret t-test results

• Be able to accurately present t-test results in research reports

5.3 Lectures and slides

5.3.1 Lectures

Watch Lecture week 5 part 1:

Watch Lecture week 5 part 2:

WarningCorrection: effect-size conventions (part 2, slide 30, 22 minutes into video)

Slide 30 gives 0.1, 0.3 and 0.5 as small, medium and large effects. These are Cohen’s conventions for correlations (r). For Cohen’s d, the conventions are 0.2 (small), 0.5 (medium) and 0.8 (large) (Cohen, 1992). Please choose the appropriate effect size for your test (d for t tests) and interpret it in the context of the measure; the labels are only a rough guide.

Take the quiz (not assessed) on the lecture materials.

5.3.2 Slides

Download the lecture slides for:

• part 1

• part 2

5.4 Practical Materials

5.4.1 Workbook

• Part 1 is some revision from the last 4 weeks.
• Part 2 covers running an independent t-test.
• Part 3 covers running a paired t-test.
• Part 4 provides more practice running paired t-tests.
• Part 5 presents some extras - exploring different datasets and running t-tests on those data.

5.4.1.1 Part 1: Revision

Task 1: Checklist: What I can now do

You should be able to answer yes to all the following. If you can’t yet, go back to the previous workbooks and repeat your working until you can answer yes, being able to type in and run the commands without referring to your notes.

1. I can open R-studio

2. I can open new libraries using library()

3. I can make an R script file

4. I can input a file into an object in R-studio using read_csv()

5. I can join two files together using inner_join()

6. I can select certain variables from an object using select()

7. I can select subsets of data using filter() (e.g., I can select participants in two conditions from a data set containing participants in four conditions)

8. I can make new variables using mutate()

9. I can arrange data according to subsets using group_by()

10. I can change format of data from wide to long format using pivot_longer

11. I can change format of data from long to wide format using pivot_wider

12. I can produce summaries of means and standard deviations for subsets of data after applying group_by() using summarise()

13. I can draw histograms of single variables, point plots of two ratio/interval/ordinal variables, bar plots of counts, and box plots of one categorical and one ratio/interval/ordinal variable using ggplot()

14. I can run a Chi-squared test and Cramer’s V test using chisq.test() and cramersV()

15. I can interpret the results of a Chi-squared test and Cramer’s V test and write up a simple report of the results.

16. I can save an R script file.

NoteANSWER

Here are some examples of the commands/functions in use:

rm(list=ls())
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# from week4 practical
dat <- read_csv("data/wk4/PSYC7142-shipley-scores-anonymous-17_26_year_language.csv") 

Rows: 325 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): english_status
dbl (6): subject_ID, Age, Shipley_Voc_Score, Gent_1_score, Gent_2_score, aca...

ℹ 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.

dat$academic_year <- as.factor(dat$academic_year)
#View(dat)

# from week2 practical
dat2 <- read_csv("data/wk2/ahicesd.csv")

Rows: 992 Columns: 50
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (50): id, occasion, elapsed.days, intervention, ahi01, ahi02, ahi03, ahi...

ℹ 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.

pinfo <- read_csv("data/wk2/participantinfo.csv")

Rows: 295 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (6): id, intervention, sex, age, educ, income

ℹ 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.

all_dat <- inner_join(x = dat2, y = pinfo, by = c("id", "intervention"))
summarydata <- select(all_dat, id, ahiTotal, cesdTotal, age, occasion)
dat_17 <- filter(dat, academic_year == "201718")
all_dat2 <- group_by(all_dat, intervention, occasion)
summarise(all_dat2, mean(ahiTotal), sd(ahiTotal), n())

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by intervention and occasion.
ℹ Output is grouped by intervention.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(intervention, occasion))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 24 × 5
# Groups:   intervention [4]
   intervention occasion `mean(ahiTotal)` `sd(ahiTotal)` `n()`
          <dbl>    <dbl>            <dbl>          <dbl> <int>
 1            1        0             68.4           14.0    72
 2            1        1             69.5           13.4    30
 3            1        2             70.3           15.7    38
 4            1        3             75.0           12.8    29
 5            1        4             76.5           16.2    36
 6            1        5             75.5           14.5    27
 7            2        0             68.8           13.0    76
 8            2        1             71.6           12.5    48
 9            2        2             73.0           14.0    48
10            2        3             72.5           14.2    43
# ℹ 14 more rows

ggplot(dat, aes(x = Age)) + geom_histogram()

`stat_bin()` using `bins = 30`. Pick better value `binwidth`.

Warning: Removed 248 rows containing non-finite outside the scale range
(`stat_bin()`).

ggplot(dat, aes(x = Shipley_Voc_Score, y = Gent_1_score)) + geom_point()

Warning: Removed 11 rows containing missing values or values outside the scale range
(`geom_point()`).

ggplot(dat, aes(x = Shipley_Voc_Score, y = Gent_1_score)) + geom_point() + geom_smooth( method = lm)

`geom_smooth()` using formula = 'y ~ x'

Warning: Removed 11 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 11 rows containing missing values or values outside the scale range
(`geom_point()`).

all_dat$occasion <- as.factor(all_dat$occasion)
ggplot(all_dat, aes(x = occasion, y = ahiTotal)) + geom_boxplot()

dat_long <- pivot_longer(dat, names_to = "test", values_to = "score", cols = c("Gent_1_score", "Gent_2_score"))
# or:
dat_long <- pivot_longer(dat, names_to = "test", values_to = "score", cols = starts_with("Gent"))
dat_wide <- pivot_wider(dat_long, names_from = "test", values_from = "score")

library(lsr)
chisq.test(x = all_dat$occasion, y = all_dat$intervention)

    Pearson's Chi-squared test

data:  all_dat$occasion and all_dat$intervention
X-squared = 9.7806, df = 15, p-value = 0.8333

cramersV(x = all_dat$occasion, y = all_dat$intervention)

[1] 0.05732779

5.4.1.2 Part 2: Running an independent t-test

Task 2: Load, prepare, and explore the data

17. Clear out R using rm(list=ls())

18. Load again the data set on the Shipley and Gent vocabulary scores from week 4.

19. Set the research question: do people who self-identify as male or female have different scores on the Gent vocabulary test? The research hypothesis is: “People who identify as male or female have different vocabulary scores”. What is the null hypothesis?

NoteAnswer

There is no difference between people who self-identify as male or female on vocabulary scores.

20. To test the research hypothesis, we will filter people who self-identify as male or female from the data set. To be inclusive, additional research questions would be part of your research project to analyse also people who self-identify as other gender. Run this command to extract a subset of the data (note that the | stands for “or”, and means Gender matches male or gender matches female:

dat <- read_csv("data/wk4/PSYC7142-shipley-scores-anonymous-17_26_gender.csv")

Rows: 325 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Gender
dbl (4): subject_ID, Shipley_Voc_Score, Gent_1_score, Gent_2_score

ℹ 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.

dat2 <- filter(dat, Gender == 'Male' | Gender == 'Female')

21. Draw a box plot of Gent vocabulary test 1 scores by gender. For a box plot, note that we need data in “long format”, where each observation is on one line, and we have a column that indicates which condition (in this case Gender) the participant is in. Does it look like there might be a gender effect? What is the direction of the effect?

NoteAnswer

ggplot(dat2, aes(x = Gender, y = Gent_1_score)) + geom_boxplot()

Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_boxplot()`).

The graph indicates that maybe males are slightly higher, but looks like a lot of overlap.

22. Note that unless we had filtered the data, the box plot would contain ‘NA’ as well, which stands for missing data. In a data set it’s always a good idea to call missing data ‘NA’ rather than just leaving them blank because this could be interpreted as a zero or as an error of filling in data. Missing values make things untidy, so it’s good practice to focus only on the variables we need for the t-test and remove all other missing values. Use select() to get just the Gender and Gent_1_score variables, and put this in a new object called ‘dat3’.

NoteAnswer

dat3 <- select(dat2, Gender, Gent_1_score)

23. Next, in order to run a t-test we have to remove any rows of data which contain a ‘NA’ - either in the Gender or the Gent_1_score variables. We do this using drop_na(dat3), put the result in a new object called ‘dat4’. Run this command:

dat4 <- drop_na(dat3)

24. Now, redraw the box plot from Step 21. Check there are just two groups still.

NoteAnswer

ggplot(dat2, aes(x = Gender, y = Gent_1_score)) + geom_boxplot()

Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Yes, two groups.

25. Compute mean and SDs for people who self-identify as male or female on Gent vocabulary test 1 scores.

WarningHint

Use group_by() and summarise().

NoteAnswer

dat5 <- group_by(dat4, Gender)
summarise(dat5, mean(Gent_1_score), sd(Gent_1_score))

# A tibble: 2 × 3
  Gender `mean(Gent_1_score)` `sd(Gent_1_score)`
  <chr>                 <dbl>              <dbl>
1 Female                 59.4               13.3
2 Male                   62.4               13.5

# Or, if you know how to use the %>% (pipe) ...:
dat4 %>% group_by(Gender) %>% summarise(mean(Gent_1_score), sd(Gent_1_score))

# A tibble: 2 × 3
  Gender `mean(Gent_1_score)` `sd(Gent_1_score)`
  <chr>                 <dbl>              <dbl>
1 Female                 59.4               13.3
2 Male                   62.4               13.5

Task 3: Run the independent t-test and measure effect size

26. Conduct an independent t-test using this command:

t.test(Gent_1_score ~ Gender, data = dat4 ) 

• ‘Gent_1_score ~ Gender’ : the ~ can be interpreted as ‘by’, i.e., compute Gent_1_score by Gender

27. The results should look like this, do yours?

    Welch Two Sample t-test

data:  Gent_1_score by Gender
t = -1.6353, df = 106.13, p-value = 0.1049
alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
95 percent confidence interval:
 -6.673859  0.640606
sample estimates:
mean in group Female   mean in group Male 
            59.35102             62.36765 

28. The key part of the results to look at is the one that has t = -1.6353, df = 106.13, p-value = 0.1049. This is the result that you report: t(106.13) = -1.64, p = .105.

The value is negative because the function includes Female before Male - and Female score is lower than Male score. What matters is how far away from zero the t-test is (either positively or negatively). The df value is slightly odd because the t.test() function figures out degrees of freedom in a technical way which takes into account differences in variance in the data between the two groups. We can just use the value that the t.test() function gives us.

29. Is this a significant difference?

NoteAnswer

No, it isn’t. The p-value is greater than 0.05.

30. Now we need to compute the effect size, using Cohen’s d. You need to load the library lsr then use this command:

cohensD(Gent_1_score ~ Gender, method = "unequal", data = dat4)

• It’s pretty much the same as the t-test() command except that we use ‘method = ’unequal’. For a paired t-test you would use ‘method = ’paired’

31. What is the effect size? Make a brief report of the results - reporting means and SDs, the t-test, p-value, and Cohen’s d. Discuss your brief report in your group.

NoteAnswer

d = 0.22

An ideal brief report of the results would state what the research question/hypothesis is, describe the means and SDs of the groups being compared, and then report the t-test statistic, with p-value and Cohen’s d effect size, then provide a brief interpretation of what the results mean.

Based on Hyde and Linn (1988), we hypothesised that people who self-identify as female may score slightly higher than males in terms of vocabulary scores. Males (mean = 62.37, SD = 13.50) scored higher than females (mean = 59.35, SD = 13.30) on the first time participants attempted the Gent vocabulary test, however this difference was not significant, t(106.13) = -1.64, p = .105, Cohen’s d = 0.22. We did not find evidence to support the hypothesis.

Reference

Hyde, J. S., & Linn, M. C. (1988). Gender differences in verbal ability: A meta-analysis. Psychological Bulletin, 104(1), 53–69. https://doi.org/10.1037/0033-2909.104.1.53

32. Make sure all commands are in the source window, save them as a new R script file.

Task 4: Practise running another independent t-test

33. Next research question: do people who are native English speakers have different vocabulary scores than those who learned English as a second language? What is the research hypothesis and the null hypothesis?

NoteAnswer

Research Hypothesis: People who speak English as a native language have higher vocabulary scores than those with English as a second language.

Null Hypothesis: There is no difference in vocabulary scores between native English and second language English speakers.

34. Repeat the Steps 22-30 in Tasks 2 and 3 except using english_status in place of Gender throughout.

NoteAnswer

dat <- read_csv("data/wk4/PSYC7142-shipley-scores-anonymous-17_26_year_language.csv")

Rows: 325 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): english_status
dbl (6): subject_ID, Age, Shipley_Voc_Score, Gent_1_score, Gent_2_score, aca...

ℹ 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.

dat2 <- select(dat, english_status, Gent_1_score)
dat3 <- drop_na(dat2)
ggplot(dat3, aes(x = english_status, y = Gent_1_score)) + geom_boxplot()

dat4 <- group_by(dat3, english_status)
summarise(dat4, mean(Gent_1_score), sd(Gent_1_score), n())

# A tibble: 2 × 4
  english_status `mean(Gent_1_score)` `sd(Gent_1_score)` `n()`
  <chr>                         <dbl>              <dbl> <int>
1 ESL                            50.8              15.2    112
2 native                         65.2               8.87   201

t.test(Gent_1_score ~ english_status, data = dat3 )

    Welch Two Sample t-test

data:  Gent_1_score by english_status
t = -9.1635, df = 154.21, p-value = 2.963e-16
alternative hypothesis: true difference in means between group ESL and group native is not equal to 0
95 percent confidence interval:
 -17.42638 -11.24535
sample estimates:
   mean in group ESL mean in group native 
            50.84821             65.18408 

cohensD(Gent_1_score ~ english_status, method = "unequal", data = dat3 )

[1] 1.153389

35. Write a brief report of the results, including means and SDs for native speakers and ESL speakers, t-test, p-value, and Cohen’s d. Discuss your report in your group.

NoteAnswer

We hypothesised that English as an additional language would result in lower overall vocabulary scores than English as a first language, because of the relation between length of time learning a language and language skills (Davies et al., 2017). Native English speakers (mean = 65.18, SD = 8.87) scored significantly higher on the Gent vocabulary test than speakers of English as a second language (mean = 50.85, SD = 15.17), t(154.21) = -9.16, p < .001, Cohen’s d = 1.15, a large effect size. Thus, the hypothesis was supported, with native English speakers scoring significantly higher on the Gent vocabulary test than those with English as an additional language.

Reference

Davies, R. A. I., Arnell, R., Birchenough, J., Grimmond, D., & Houlson, S. (2017). Reading Through the Life Span: Individual Differences in Psycholinguistic Effects. Journal of Experimental Psychology: Learning, Memory, and Cognition, 43(8), 1298-1338. https://doi.org/10.1037/xlm0000366

36. Save your R script file.

5.4.1.3 Part 3: Conducting a paired t-test

Task 5: Conducting a paired t-test

37. Clear out R-studio before we get started again using rm(list=ls())

38. We are going to investigate again the data from this paper: Woodworth, R.J., O’Brien-Malone, A., Diamond, M.R. and Schuez, B., 2017. Data from, “Web-based Positive Psychology Interventions: A Reexamination of Effectiveness”. Journal of Open Psychology Data, 6(1).

Our research question is whether happiness scores are affected by the interventions. We will look at the pre-test (occasion 0) and the first test after the intervention (occasion 1).

39. What is the research hypothesis and what is the null hypothesis?

NoteAnswer

RH: happiness scores change (increase) from the first to the second occasion of testing.

NH: happiness scores do not change across occasions.

40. For a paired t-test we can only include data from people who have produced scores at both occasions of testing. So, we need a slightly different version of the data, which you can download here for the ahicesd.csv file and here for the participantinfo2.csv file.

Remind yourself what these data mean.

41. Once again, join the ahicesd.csv and participantinfo2.csv data in R-studio by aligning the names for the participant numbers in these two data sets (see week 2 workbook for reminders about this).

NoteAnswer

library(tidyverse)
dat <- read_csv("data/wk5/ahicesd.csv")

Rows: 992 Columns: 50
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (50): id, occasion, elapsed.days, intervention, ahi01, ahi02, ahi03, ahi...

ℹ 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.

pinfo <- read_csv("data/wk5/participantinfo2.csv")

New names:
Rows: 147 Columns: 7
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(7): ...1, id, intervention, sex, age, educ, income
ℹ 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.
• `` -> `...1`

all_dat <- inner_join(x = dat, y = pinfo, by = c("id", "intervention"))

42. Let’s select only the relevant variables. Use select() to select only id, ahiTotal, and occasion variables, and save this as a new object called ‘summarydata’

NoteAnswer

summarydata <- select(all_dat, id, ahiTotal, occasion)

43. Use filter to pull out only occasion == 0 or occasion == 1 scores

WarningHint

use occasion == 0 | occasion == 1, save this as a new object called summarydata2

NoteAnswer

summarydata2 <- filter(summarydata, occasion == 0 | occasion == 1)

44. Here is where we would usually remove all the NA values, but there aren’t any in this file (so we don’t need drop_na()).

45. Now, we need to make sure occasion is treated as a categorical variable, rather than a continuous variable, so we need to convert it to a factor:

summarydata2$occasion <-as.factor(summarydata2$occasion)

46. Now, draw a box plot of ahiTotal scores by occasion (why do we use a box plot?)

NoteAnswer

ggplot(summarydata2, aes(x = occasion, y = ahiTotal)) + geom_boxplot()

We use a boxplot because we have one nominal/categorical variable and one interval/ratio/ordinal measure.

47. Compute mean and SD for each occasion

NoteAnswer

summarydata3 <- group_by(summarydata2, occasion)
summarise(summarydata3, mean(ahiTotal), sd(ahiTotal), n())

# A tibble: 2 × 4
  occasion `mean(ahiTotal)` `sd(ahiTotal)` `n()`
  <fct>               <dbl>          <dbl> <int>
1 0                    69.3           12.3   147
2 1                    72.4           12.6   147

# or:
group_by(summarydata2, occasion) %>% summarise(mean(ahiTotal), sd(ahiTotal), n())

# A tibble: 2 × 4
  occasion `mean(ahiTotal)` `sd(ahiTotal)` `n()`
  <fct>               <dbl>          <dbl> <int>
1 0                    69.3           12.3   147
2 1                    72.4           12.6   147

48. In order to run the paired t-test, we first need to make sure that the paired values (in this case the measures from the same person) are on the same row of the data. So, let’s use pivot_wider to put the ahiTotal scores for occasion 0 in one column, and the ahiTotal scores for occasion 1 in another column:

summarydata2_wide <- pivot_wider(summarydata2, names_from = occasion, values_from = ahiTotal, names_prefix = "occasion_")

NoteNote

Note that we use the names_prefix = “occasion_” as an extra argument to make the column names “occasion_0” and “occasion_1”.

What would happen if we didn’t use names_prefix? It would still work, but it’s a bit more awkward to work with the output.

49. Then, we can run the paired t-test in the following way:

t.test(summarydata2_wide$occasion_0, summarydata2_wide$occasion_1, paired = TRUE)

Is the result significant?

NoteAnswer

yes, because p < .05.

50. Now run Cohen’s d: it’s similar to the way we do the paired t-test:

cohensD( summarydata2_wide$occasion_0, summarydata2_wide$occasion_1, method =  "paired")

What is the value for Cohen’s d?

NoteAnswer

d = 0.4059904

51. Write up a brief report of the result and discuss in your group.

NoteAnswer

We tested whether participants’ happiness scores at first testing after the interventions were different than their scores prior to the interventions. We found that, prior to the intervention, scores were significantly lower (M = 69.3, SD = 12.3) than they were immediately after the interventions (M = 72.4, SD = 12.6), t(146) = -4.92, p < .001, Cohen’s d = 0.41, a small effect. The results indicate that the intervention had a positive effect on happiness scores.

52. Save your R script file.

5.4.1.4 Part 4: More practice running a paired t-test

We are going to figure out whether people have different scores the first and second time they take the Gent vocabulary test.

53. Go back to the vocabulary scores data. Load the data into dat, and make another object dat2 that contains only the subject_ID, Gent_1_score and Gent_2_score.

54. Some people did not do all the tests - look at participant 46 for instance. To do a t-test we need data where the person does both tests. We can filter out the scores where there are no NAs by repeating the drop_na we did at step 23, above. Call the new data object dat3.

NoteAnswer

dat <- read_csv("data/wk4/PSYC7142-shipley-scores-anonymous-17_26_year_language.csv")

Rows: 325 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): english_status
dbl (6): subject_ID, Age, Shipley_Voc_Score, Gent_1_score, Gent_2_score, aca...

ℹ 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.

dat$academic_year <- as.factor(dat$academic_year)
dat2 <- select(dat, subject_ID, Gent_1_score, Gent_2_score)
dat3 <- drop_na(dat2)

55. Run the paired t-test.

Each row contains a score for each person for each Gent test, and so we are ready to run the paired t-test:

NoteAnswer

t.test(dat3$Gent_1_score, dat3$Gent_2_score, paired = TRUE)

    Paired t-test

data:  dat3$Gent_1_score and dat3$Gent_2_score
t = -2.3561, df = 301, p-value = 0.01911
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -2.1451382 -0.1926101
sample estimates:
mean difference 
      -1.168874 

cohensD(dat3$Gent_1_score, dat3$Gent_2_score, method =  "paired")

[1] 0.1355798

56. Get mean and SD of the scores for the tests.

NoteAnswer

summarise(dat3, mean(Gent_1_score), sd(Gent_1_score), mean(Gent_2_score), sd(Gent_2_score))

# A tibble: 1 × 4
  `mean(Gent_1_score)` `sd(Gent_1_score)` `mean(Gent_2_score)`
                 <dbl>              <dbl>                <dbl>
1                 60.0               13.2                 61.1
# ℹ 1 more variable: `sd(Gent_2_score)` <dbl>

57. In order to draw a box plot of the Gent vocabulary scores taken at the first and second occasion, we will need the data in “long format”, so that there is a column saying which test it is and a column reporting the scores on the Gent tests. Apply the pivot_longer command and then draw the box plot.

NoteAnswer

dat4 <- pivot_longer(dat3, names_to = "test", values_to = "score", cols = c("Gent_1_score", "Gent_2_score"))
ggplot(dat4, aes(x = test, y = score)) + geom_boxplot()

5.4.1.5 Part 5: Extras

58. In the Gent vocabulary scores data, is there a significant difference between native and ESL speakers for your academic year group?

NoteAnswer

dat <- read_csv("data/wk4/PSYC7142-shipley-scores-anonymous-17_26_year_language.csv")

Rows: 325 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): english_status
dbl (6): subject_ID, Age, Shipley_Voc_Score, Gent_1_score, Gent_2_score, aca...

ℹ 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.

dat2 <- filter(dat, academic_year == "202627")
ggplot(dat2, aes(x = english_status, y = Gent_1_score)) + geom_boxplot()

dat3 <- group_by(dat2, english_status)
summarise(dat3, mean(Gent_1_score), sd(Gent_1_score), mean(Gent_2_score), sd(Gent_2_score), n())

# A tibble: 2 × 6
  english_status `mean(Gent_1_score)` `sd(Gent_1_score)` `mean(Gent_2_score)`
  <chr>                         <dbl>              <dbl>                <dbl>
1 ESL                            59                11.3                  61.2
2 native                         65.5               6.48                 65.8
# ℹ 2 more variables: `sd(Gent_2_score)` <dbl>, `n()` <int>

# or:
dat2 %>% group_by(english_status) %>% summarise(mean(Gent_1_score), sd(Gent_1_score), mean(Gent_2_score), sd(Gent_2_score), n())

# A tibble: 2 × 6
  english_status `mean(Gent_1_score)` `sd(Gent_1_score)` `mean(Gent_2_score)`
  <chr>                         <dbl>              <dbl>                <dbl>
1 ESL                            59                11.3                  61.2
2 native                         65.5               6.48                 65.8
# ℹ 2 more variables: `sd(Gent_2_score)` <dbl>, `n()` <int>

t.test(Gent_1_score ~ english_status, data = dat2 )

    Welch Two Sample t-test

data:  Gent_1_score by english_status
t = -1.708, df = 11.776, p-value = 0.1138
alternative hypothesis: true difference in means between group ESL and group native is not equal to 0
95 percent confidence interval:
 -14.912747   1.821838
sample estimates:
   mean in group ESL mean in group native 
            59.00000             65.54545 

cohensD(Gent_1_score ~ english_status, method = "unequal", data = dat2 )

[1] 0.7104731

t.test(Gent_2_score ~ english_status, data = dat2 )

    Welch Two Sample t-test

data:  Gent_2_score by english_status
t = -1.5188, df = 15.201, p-value = 0.1493
alternative hypothesis: true difference in means between group ESL and group native is not equal to 0
95 percent confidence interval:
 -10.982478   1.837023
sample estimates:
   mean in group ESL mean in group native 
            61.20000             65.77273 

cohensD(Gent_2_score ~ english_status, method = "unequal", data = dat2 )

[1] 0.5963303

t.test(Gent_2_score ~ english_status, data = dat2 )

    Welch Two Sample t-test

data:  Gent_2_score by english_status
t = -1.5188, df = 15.201, p-value = 0.1493
alternative hypothesis: true difference in means between group ESL and group native is not equal to 0
95 percent confidence interval:
 -10.982478   1.837023
sample estimates:
   mean in group ESL mean in group native 
            61.20000             65.77273 

cohensD(Gent_2_score ~ english_status, method = "unequal", data = dat2 )

[1] 0.5963303

nothing significant…

59. Are there significant differences for the Shipley vocabulary test measures between males and females, or between those with English as first or second language?

NoteAnswer

dat3 <- group_by(dat2, english_status)
summarise(dat3, mean(Shipley_Voc_Score, na.rm=TRUE), sd(Shipley_Voc_Score, na.rm=TRUE), n())

# A tibble: 2 × 4
  english_status mean(Shipley_Voc_Score, na.rm = …¹ sd(Shipley_Voc_Score…² `n()`
  <chr>                                       <dbl>                  <dbl> <int>
1 ESL                                          31.2                   5.65    10
2 native                                       34.3                  13.0     22
# ℹ abbreviated names: ¹​`mean(Shipley_Voc_Score, na.rm = TRUE)`,
#   ²​`sd(Shipley_Voc_Score, na.rm = TRUE)`

# or:
dat2 %>% group_by(english_status) %>% summarise(mean(Shipley_Voc_Score, na.rm=TRUE), sd(Shipley_Voc_Score, na.rm=TRUE), n())

# A tibble: 2 × 4
  english_status mean(Shipley_Voc_Score, na.rm = …¹ sd(Shipley_Voc_Score…² `n()`
  <chr>                                       <dbl>                  <dbl> <int>
1 ESL                                          31.2                   5.65    10
2 native                                       34.3                  13.0     22
# ℹ abbreviated names: ¹​`mean(Shipley_Voc_Score, na.rm = TRUE)`,
#   ²​`sd(Shipley_Voc_Score, na.rm = TRUE)`

t.test(Shipley_Voc_Score ~ english_status, data = dat2)

    Welch Two Sample t-test

data:  Shipley_Voc_Score by english_status
t = -0.93043, df = 29.993, p-value = 0.3596
alternative hypothesis: true difference in means between group ESL and group native is not equal to 0
95 percent confidence interval:
 -9.817334  3.671880
sample estimates:
   mean in group ESL mean in group native 
            31.20000             34.27273 

cohensD(Shipley_Voc_Score ~ english_status, method = "unequal", data = dat2 )

[1] 0.3060571

For the eagle-eyed among you, you might notice that some of the scores are outside the possible range - these typos often occur in data, and should be filtered out. What is the result if you do the filtering?

NoteAnswer

Shipley vocabulary scores cannot be higher than 40, so remove any scores above 40:

dat4 <- filter(dat2, Shipley_Voc_Score <= 40)
dat5 <- group_by(dat4, english_status)
summarise(dat5, mean(Shipley_Voc_Score), sd(Shipley_Voc_Score), n())

# A tibble: 2 × 4
  english_status `mean(Shipley_Voc_Score)` `sd(Shipley_Voc_Score)` `n()`
  <chr>                              <dbl>                   <dbl> <int>
1 ESL                                 31.2                    5.65    10
2 native                              31.7                    5.19    21

t.test(Shipley_Voc_Score ~ english_status, data = dat4)

    Welch Two Sample t-test

data:  Shipley_Voc_Score by english_status
t = -0.24305, df = 16.474, p-value = 0.811
alternative hypothesis: true difference in means between group ESL and group native is not equal to 0
95 percent confidence interval:
 -4.989396  3.960825
sample estimates:
   mean in group ESL mean in group native 
            31.20000             31.71429 

cohensD(Shipley_Voc_Score ~ english_status, method = "unequal", data = dat4)

[1] 0.09479229

Removing the one impossible score (88) brings the native speakers’ SD down from 13.0 to 5.2, and the difference between the groups all but disappears.

60. The data from this paper are called ClassDraw.csv available here. The data are on osf as well, but they’re not formatted in a very accessible way, so I adjusted them and put them here.

Jalava, S. T., Wammes, J. D., & Cheng, K. (2023). Drawing your way to an A: Long-lasting improvements in classroom quiz performance following drawing. Psychonomic Bulletin & Review, 30, 1939–1945. https://doi.org/10.3758/s13423-023-02294-2

There are some useful tips in the results of this study about the benefit of doodling…

My challenge to you:

• Can you make a ggplot that looks a bit like Figure 2 from this study?

NoteAnswer

jalava <- read_csv("data/wk5/ClassDraw.csv")

Rows: 168 Columns: 9
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Gender
dbl (8): ID, Age, write_1, write_2, draw_1, draw_2, exam_draw, exam_write

ℹ 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.

jalava1 <- pivot_longer(jalava, names_to = c("test",".value"), names_sep="_" , values_to = "score", cols = c("write_1","write_2","draw_1","draw_2"))
jalava2 <- pivot_longer(jalava1, names_to = "test_time", values_to = "score", cols = c(`1`,`2`))

jalava2$test_time <- as.factor(jalava2$test_time)
ggplot(jalava2, aes(x = test_time, y = score, fill = test)) + geom_boxplot()

Warning: Removed 48 rows containing non-finite outside the scale range
(`stat_boxplot()`).

# median looks very similar (boxplots very overlapping), but means are a bit different:
jalava2 %>% group_by(test, test_time) %>% summarise(mean(score, na.rm=TRUE), sd(score, na.rm=TRUE))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by test and test_time.
ℹ Output is grouped by test.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(test, test_time))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 4 × 4
# Groups:   test [2]
  test  test_time `mean(score, na.rm = TRUE)` `sd(score, na.rm = TRUE)`
  <chr> <fct>                           <dbl>                     <dbl>
1 draw  1                               0.486                     0.215
2 draw  2                               0.322                     0.216
3 write 1                               0.467                     0.212
4 write 2                               0.3                       0.212

61. If you’ve made some progress on the de Zubicaray et al. (2024) data analysis, to the point where you have replicated (more or less) the multiple regression analysis from the first study in that paper, then here is the next step.

Note, this is a major challenge, but it will be a step into new research!

The task is to take the de Zubicaray data, and add a new variable which determines whether the word is a palindrome or not. A palindrome is a word that is spelled the same backwards as forwards (e.g., abba, rotavator).

I would like to know if being a palindrome or not affects processing of the words. My prediction is that being a palindrome helps access the word.

What you will need to do is apply a function called “palindrome” to the words in the de Zubicaray dataset, to make a new variable, called or some other name then run a multiple regression with all the other variables plus this one.

Here is a script with the function called palindrome included, with some comments to show how it is put into practice. You will need to add the function into the top of an analysis script (e.g., rmd or r script) file and make sure that is loaded along with the other commands you are using.

5.4.2 Data

Data referred to in this workbook:

• PSYC7142-shipley-scores-anonymous-17_26_gender.csv
• PSYC7142-shipley-scores-anonymous-17_26_year_language.csv
• ahicesd.csv adapted for t-test
• participantinfo2.csv adapted for t-test.
• ClassDraw.csv

5.4.3 Answers

The answers to the workbook appear below each question in the workbook, above, so you can check your work after you’ve had a go at the exercises..!.

5.5 Extras

Also optionally, read the articles on the importance of statistical understanding and insights from good data visualisation:

• How scientists can be better at statistics. Note that this article is hosted on the Spectator website, and (content warning) refers to Harold Shipman, a serial killer.

• Florence Nightingale and data visualisation. Note that this is hosted on Scientific American.