Simulation Practice

In this brief tutorial, we’ll review some basic R functions we’ve learned and learn how to do some simple simulation of data.

library(tidyverse)

Let’s imagine a pretest-posttest design with an experimental and control group. We’ll say that the treatment involves learning vocabulary while reading with interactive glosses. The control group gets no glosses.

Since we’re simulating, we can have lots of participants for free (yay!). Let’s make a list of 100 participants.

ID <- 1:100

Now let’s put them into groups. First we’ll make the labels using the rep() function and join them with c():

group <- c(rep("treatment", 50), rep("control", 50))

Now we’ll make a data frame with both variables together (a tibble, actually):

d <- tibble(ID, group)

Let’s imagine that at the pretest, there’s really no systematic difference between the two groups in their vocabulary knowledge. And for simplicity’s sake, we’ll imagine that our vocabulary measure is a test with 100 possible points. These are learners, so they shouldn’t have extremely high scores to start out.

set.seed(2021) #use any number - this is just for reproducibility

d$pretest <- rnorm(100, m = 30, sd = 5)

Let’s take a look at the score distribution:

d %>% ggplot(aes(x=pretest))+
  geom_histogram(binwidth = 5)+
  scale_x_continuous(limits = c(0, 100))+
  theme_bw()

Let’s see if the groups look any different:

d %>% ggplot(aes(x=pretest))+
  geom_histogram(binwidth = 5)+
  scale_x_continuous(limits = c(0, 100))+
  theme_bw()+
  facet_wrap(~group, nrow = 1)

And a numeric summary….

d %>% group_by(group) %>% summarise(mean = mean(pretest),
                                    sd = sd(pretest),
                                    range = range(pretest))

`summarise()` has grouped output by 'group'. You can override using the `.groups` argument.

Now we’ll generate a posttest score, but let’s imagine that the treatment is quite effective - students reading with interactive glosses learn a lot more vocabulary than those without.

d$posttest <- c((rnorm(50, m = 50, sd = 7)+7), rnorm(50, m = 50, sd = 7))

Now let’s take a look at the posttest histograms:

d %>% ggplot(aes(x=posttest))+
  geom_histogram(binwidth = 5)+
  scale_x_continuous(limits = c(0, 100))+
  theme_bw()

d %>% ggplot(aes(x=posttest))+
  geom_histogram(binwidth = 5)+
  scale_x_continuous(limits = c(0, 100))+
  theme_bw()+
  facet_wrap(~group, nrow = 1)

d %>% group_by(group) %>% summarise(mean = mean(posttest),
                                    sd = sd(posttest),
                                    range = range(posttest))

`summarise()` has grouped output by 'group'. You can override using the `.groups` argument.

Let’s make this data tidy by putting it in long format:

d <- d %>% pivot_longer(pretest:posttest, names_to = "time", values_to = "vocab")

d <- d %>% mutate(time = factor(time, levels = c("pretest", "posttest")))

Take a look at the data - you should now see two rows per participant. This will open up some neat plotting options for us.

d %>% ggplot(aes(x=vocab))+
  geom_histogram(binwidth = 5)+
  scale_x_continuous(limits = c(0, 100))+
  theme_bw()+
  facet_wrap(~group+time, nrow = 2)

Or this…

d %>% ggplot(aes(x=vocab, fill = time))+
  geom_histogram(binwidth = 5)+
  scale_x_continuous(limits = c(0, 100))+
  theme_bw()+
  facet_wrap(~group, nrow = 2)

And to summarize numerically:

d %>% group_by(group, time) %>% summarise(mean = mean(vocab),
                                    sd = sd(vocab),
                                    min = min(vocab),
                                    max = max(vocab))

`summarise()` has grouped output by 'group'. You can override using the `.groups` argument.