Exercise 15: Power analyses

Exercise 15: Power analyses#

This assignment is designed to give you practice with Monte Carlo methods to conduct power analyses via simulation. You won’t need to load in any data for this homework. We will, however, be using parts of the homework from last week.

1. Simulating data (1 points)#

Pull your simulate_data() function from your last homework and add it below.

As a reminder, this function simulates the relationship between age, word reading experience, and reading comprehension skill.

c is reading comprehension, and x is word reading experience.

sample_size = 100 # How many children in data set? 
age_lo = 80     # minimum age, in months
age_hi = 200    # maximum age, in months
beta_xa = 0.5   # amount by which experience changes for increase of one month in age
beta_x0 = -5    # amount of experience when age = 0 (not interpretable, since minimum age for this data is 80 months)
sd_x = 50       # standard dev of gaussian noise term, epsilon_x
beta_ca = 0.8   # amount that comprehension score improves for every increase of one unit in age
beta_cx = 3     # amount that comprehension score improves for every increase of one unit in reading experience
beta_c0 = 10    # comprehension score when reading experience is 0. 
sd_c = 85      # standard dev of gaussian noise term, epsilon_c

simulate_data <- function(sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c) {
      # WRITE YOUR CODE HERE
     
      return(data.frame(age=age,x=x,c=c)) # it's actually bad form to have a variable named "c" in R, my bad...
}

dat <- simulate_data(sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c)
head(dat)

2. `run_analysis()` function (2 pts)#

Last week, we looked at whether word reading experience(x) mediated the relation between age and reading comprehension (c).

Now we’re going to use our simulate_data() function to conduct a power analysis. The goal is to determine how many participants we would need in order to detect both the mediated and the direct effects in this data.

Note: We’re going to pretend for the sake of simplicity that we don’t have any control over the ages of the children we get (so ages are generated using runif(sample_size, age_lo, age_hi), although of course this would be an unusual situation in reality.

First, write a function, run_analysis(), that takes in simulated data, runs your mediation from last week, and returns a vector containing the ACME and ADE estimates and p-values (these are the d0, d0.p, z0, and z0.p features of the mediated model object, e.g., fitMed$d0.p). Print this function’s output for the data we simulated previously.

# WRITE YOUR CODE HERE
library(mediation)

3. `repeat_analysis()` function (3 pts)#

Next fill in the function repeat_analysis() below so that it simulates and analyzes data num_simulations times. Store the outputs from each simulation in the simouts matrix. Calculate and return the coverage across all the simulations run for both ACME and ADE.

repeat_analysis <- function(num_simulations, alpha, sample_size, age_lo, age_hi, 
        beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c) {
    # Initialize simouts matrix for storing each output from run_analysis()
    simouts <- matrix(rep(NA, num_simulations*4), nrow=num_simulations, ncol=4)
    
    # Start simulating
    for (i in 1:num_simulations) {
      # WRITE YOUR CODE HERE
     
    }
    
    # Calculate coverage for both ACME and ADE estimates using p-values in simouts
    ACME_cov = # FILL THIS IN 
    ADE_cov =  # FILL THIS IN

    return(list(ACME_cov = ACME_cov, ADE_cov = ADE_cov))
}

Now run the repeat_analysis() function using the same parameter settings as above, for 10 simulations, with an alpha criterion of 0.01.

# WRITE YOUR CODE HERE

4. Testing different sample sizes (2 pts)#

Finally, do the same thing (10 simulations, alpha criterion of 0.01) but for 5 different sample sizes: 50, 75, 100, 125, 150. You can do this using map (as in the tutorial), or a simple for loop, or by calculating each individually. Up to you! This should take around 3 minutes to run.

# WRITE YOUR CODE HERE

Print your results.

# WRITE YOUR CODE HERE

5. Reflection (2 pts)#

If this were a real power analysis, we’d want to run more simulations per sample size (to get a more precise estimate of power) and we may also want to test out some other values of the parameters we used to simulate our data. However, what would you conclude just based on the results above?

Write your response here.

Given how we generated the data, why was the direct effect harder to detect than the mediated effect?

Write your response here.

DUE: 5pm EST, April 3, 2024

IMPORTANT Did you collaborate with anyone on this assignment? If so, list their names here.

Someone’s Name