Exercise15.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/choiitea/Choi_DSPN_S24/blob/main/Exercise15.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-2W919d2ZXp7"
      },
      "source": [
        "# Exercise 15: Power analyses"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L4nOzVhyZXqK"
      },
      "source": [
        "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."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9DsyBTB6ZXqN"
      },
      "source": [
        "---\n",
        "## 1. Simulating data (1 points)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J1f6MywYhyzN"
      },
      "source": [
        "Pull your `simulate_data()` function from your last homework and add it below.\n",
        "\n",
        "As a reminder, this function simulates the relationship between age, word reading experience, and reading comprehension skill.\n",
        "\n",
        "`c` is reading comprehension, and `x` is word reading experience."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "UnBVazYfZXqP",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 286
        },
        "outputId": "53b51b58-c99c-4f43-be7c-49021321bdb7"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<table class=\"dataframe\">\n",
              "<caption>A data.frame: 6 × 3</caption>\n",
              "<thead>\n",
              "\t<tr><th></th><th scope=col>age</th><th scope=col>word_reading_experience</th><th scope=col>reading_comprehension</th></tr>\n",
              "\t<tr><th></th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
              "</thead>\n",
              "<tbody>\n",
              "\t<tr><th scope=row>1</th><td>116.2587</td><td>  7.788229</td><td>171.1395</td></tr>\n",
              "\t<tr><th scope=row>2</th><td>172.7556</td><td>124.824699</td><td>513.9073</td></tr>\n",
              "\t<tr><th scope=row>3</th><td>176.6749</td><td> 51.140622</td><td>157.8951</td></tr>\n",
              "\t<tr><th scope=row>4</th><td>135.2671</td><td> 43.255888</td><td>234.3591</td></tr>\n",
              "\t<tr><th scope=row>5</th><td>100.8476</td><td> 88.163574</td><td>334.0828</td></tr>\n",
              "\t<tr><th scope=row>6</th><td>177.5615</td><td> 75.574170</td><td>425.3733</td></tr>\n",
              "</tbody>\n",
              "</table>\n"
            ],
            "text/markdown": "\nA data.frame: 6 × 3\n\n| <!--/--> | age &lt;dbl&gt; | word_reading_experience &lt;dbl&gt; | reading_comprehension &lt;dbl&gt; |\n|---|---|---|---|\n| 1 | 116.2587 |   7.788229 | 171.1395 |\n| 2 | 172.7556 | 124.824699 | 513.9073 |\n| 3 | 176.6749 |  51.140622 | 157.8951 |\n| 4 | 135.2671 |  43.255888 | 234.3591 |\n| 5 | 100.8476 |  88.163574 | 334.0828 |\n| 6 | 177.5615 |  75.574170 | 425.3733 |\n\n",
            "text/latex": "A data.frame: 6 × 3\n\\begin{tabular}{r|lll}\n  & age & word\\_reading\\_experience & reading\\_comprehension\\\\\n  & <dbl> & <dbl> & <dbl>\\\\\n\\hline\n\t1 & 116.2587 &   7.788229 & 171.1395\\\\\n\t2 & 172.7556 & 124.824699 & 513.9073\\\\\n\t3 & 176.6749 &  51.140622 & 157.8951\\\\\n\t4 & 135.2671 &  43.255888 & 234.3591\\\\\n\t5 & 100.8476 &  88.163574 & 334.0828\\\\\n\t6 & 177.5615 &  75.574170 & 425.3733\\\\\n\\end{tabular}\n",
            "text/plain": [
              "  age      word_reading_experience reading_comprehension\n",
              "1 116.2587   7.788229              171.1395             \n",
              "2 172.7556 124.824699              513.9073             \n",
              "3 176.6749  51.140622              157.8951             \n",
              "4 135.2671  43.255888              234.3591             \n",
              "5 100.8476  88.163574              334.0828             \n",
              "6 177.5615  75.574170              425.3733             "
            ]
          },
          "metadata": {}
        }
      ],
      "source": [
        "sample_size = 100 # How many children in data set?\n",
        "age_lo = 80     # minimum age, in months\n",
        "age_hi = 200    # maximum age, in months\n",
        "beta_xa = 0.5   # amount by which experience changes for increase of one month in age\n",
        "beta_x0 = -5    # amount of experience when age = 0 (not interpretable, since minimum age for this data is 80 months)\n",
        "sd_x = 50       # standard dev of gaussian noise term, epsilon_x\n",
        "beta_ca = 0.8   # amount that comprehension score improves for every increase of one unit in age\n",
        "beta_cx = 3     # amount that comprehension score improves for every increase of one unit in reading experience\n",
        "beta_c0 = 10    # comprehension score when reading experience is 0.\n",
        "sd_c = 85      # standard dev of gaussian noise term, epsilon_c\n",
        "\n",
        "simulate_data <- function(sample_size, age_lo, age_hi, beta_xa,\n",
        "                          beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c) {\n",
        "                            # Generate ages distributed between age_lo and age_hi\n",
        "                            ages <- runif(sample_size, min = age_lo, max = age_hi)\n",
        "                            # Calculate word reading experience (x) based on age and add Gaussian noise\n",
        "                            x <- beta_xa * ages + beta_x0 + rnorm(sample_size, mean = 0, sd = sd_x)\n",
        "                            # Calculate reading comprehension (c) based on age, word reading experience, and add Gaussian noise\n",
        "                            c <- beta_ca * ages + beta_cx * x + beta_c0 + rnorm(sample_size, mean = 0, sd = sd_c)\n",
        "                            # Return a data frame containing age, word reading experience, and reading comprehension\n",
        "                            return(data.frame(age = ages, word_reading_experience = x, reading_comprehension = c))\n",
        "                          }\n",
        "\n",
        "dat <- simulate_data(sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c)\n",
        "head(dat)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T3b3_KsHk-xD"
      },
      "source": [
        "---\n",
        "## 2. `run_analysis()` function (2 pts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TJyPIaTO8zCS"
      },
      "source": [
        "Last week, we looked at whether word reading experience(`x`) mediated the relation between `age` and reading comprehension (`c`).\n",
        "\n",
        "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.\n",
        "\n",
        "*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.*\n",
        "\n",
        "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."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "AidH_KidrX9L",
        "scrolled": true,
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "170e518d-dc5a-4bf7-8c11-89a92c91cb5c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1] 1.7805904 0.0000000 0.8975315 0.0040000\n"
          ]
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "#install.packages(\"mediation\")\n",
        "library(mediation)\n",
        "\n",
        "run_analysis <- function(simulated_data) {\n",
        "  fitM <- lm(word_reading_experience ~ age, data = simulated_data)\n",
        "\n",
        "  fitY <- lm(reading_comprehension ~ age + word_reading_experience, data = simulated_data)\n",
        "\n",
        "  # mediation\n",
        "  fitMed <- mediate(fitM, fitY, treat = \"age\", mediator = \"word_reading_experience\")\n",
        "  # vector containing ACME estimate, ACME pvalue, ADE estimate and ADE p-value\n",
        "  acme_ade_pvalues <- c(fitMed$d0, fitMed$d0.p, fitMed$z0, fitMed$z0.p)\n",
        "\n",
        "  return(acme_ade_pvalues)\n",
        "  }\n",
        "\n",
        "analysis_results <- run_analysis(dat)\n",
        "print(analysis_results)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bHxz5W2r5s79"
      },
      "source": [
        "---\n",
        "## 3. `repeat_analysis()` function (3 pts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G_jvc1XQBCpH"
      },
      "source": [
        "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."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "id": "NheOpBDsBBo9",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "repeat_analysis <- function(num_simulations, alpha, sample_size, age_lo, age_hi,\n",
        "        beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c) {\n",
        "    # Initialize simouts matrix for storing each output from run_analysis()\n",
        "    simouts <- matrix(rep(NA, num_simulations*4), nrow=num_simulations, ncol=4)\n",
        "\n",
        "    # Start simulating\n",
        "    for (i in 1:num_simulations) {\n",
        "      # WRITE YOUR CODE HERE\n",
        "      # simulate data\n",
        "      simulated_data <- simulate_data(sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c)\n",
        "      # Run analysis\n",
        "      analysis_results <- run_analysis(simulated_data)\n",
        "      # Store the outputs from each simulation in the simouts matrix\n",
        "      simouts[i, ] <- analysis_results\n",
        "      }\n",
        "      # Calculate coverage for both ACME and ADE estimates using p-values in simouts\n",
        "      ACME_cvg = mean(simouts[, 2] <= alpha)\n",
        "      ADE_cvg =  mean(simouts[, 4] <= alpha)\n",
        "\n",
        "    return(list(ACME_cvg = ACME_cvg, ADE_cvg = ADE_cvg))\n",
        "}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bM5J5pYyQoRL"
      },
      "source": [
        "Now run the `repeat_analysis()` function using the same parameter settings as above, for 10 simulations, with an alpha criterion of 0.01."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "id": "RpURGbZQQlEK",
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4b47e1cf-258b-4c80-a8f0-eef44b1565ee"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "$ACME_cvg\n",
            "[1] 0.6\n",
            "\n",
            "$ADE_cvg\n",
            "[1] 0.5\n",
            "\n"
          ]
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "num_simulations <- 10\n",
        "alpha <- 0.01\n",
        "\n",
        "sim10_results <- repeat_analysis(num_simulations, alpha, sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c)\n",
        "\n",
        "print(sim10_results)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "imIoN4w65s7_"
      },
      "source": [
        "---\n",
        "## 4. Testing different sample sizes (2 pts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CzeTaoY7SaXv"
      },
      "source": [
        "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."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "id": "mp-a19XsS8Md",
        "vscode": {
          "languageId": "r"
        }
      },
      "outputs": [],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "library(tidyverse)\n",
        "# parameters for testing sample sizes + 10 sim + alpha .01\n",
        "sample_sizes <- c(50, 75, 100, 125, 150)\n",
        "num_simulations <- 10\n",
        "alpha <- 0.01\n",
        "\n",
        "sample_size_sim10_results <- expand.grid(sample_size = sample_sizes, alpha = alpha)\n",
        "sample_size_sim10_results$id <- 1:nrow(sample_size_sim10_results)\n",
        "\n",
        "results <- sample_size_sim10_results %>%\n",
        "  nest(parameters := c(sample_size, alpha)) %>%\n",
        "    mutate(ACME_power = map(parameters, ~ repeat_analysis(num_simulations, .$alpha, .$sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c)$ACME_cvg),\n",
        "         ADE_power = map(parameters, ~ repeat_analysis(num_simulations, .$alpha, .$sample_size, age_lo, age_hi, beta_xa, beta_x0, sd_x, beta_ca, beta_cx, beta_c0, sd_c)$ADE_cvg)) %>%\n",
        "      unnest(c(parameters, ACME_power, ADE_power))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "akJwKCB7dkDw"
      },
      "source": [
        "Print your results."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "id": "gIkkWTG2aR4d",
        "scrolled": true,
        "vscode": {
          "languageId": "r"
        },
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "fe13ad02-2cea-469f-b17f-4bc317d6074f"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\u001b[90m# A tibble: 5 × 5\u001b[39m\n",
            "     id sample_size alpha ACME_power ADE_power\n",
            "  \u001b[3m\u001b[90m<int>\u001b[39m\u001b[23m       \u001b[3m\u001b[90m<dbl>\u001b[39m\u001b[23m \u001b[3m\u001b[90m<dbl>\u001b[39m\u001b[23m      \u001b[3m\u001b[90m<dbl>\u001b[39m\u001b[23m     \u001b[3m\u001b[90m<dbl>\u001b[39m\u001b[23m\n",
            "\u001b[90m1\u001b[39m     1          50  0.01        0.5       0.4\n",
            "\u001b[90m2\u001b[39m     2          75  0.01        0.8       0.6\n",
            "\u001b[90m3\u001b[39m     3         100  0.01        1         0.7\n",
            "\u001b[90m4\u001b[39m     4         125  0.01        0.9       0.8\n",
            "\u001b[90m5\u001b[39m     5         150  0.01        1         1  \n"
          ]
        }
      ],
      "source": [
        "# WRITE YOUR CODE HERE\n",
        "print(results)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q3Oo6npB5s8B"
      },
      "source": [
        "## 5. Reflection (2 pts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tpbGkul2dr-5"
      },
      "source": [
        "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?\n",
        "\n",
        "> based on the results, we can generally see that increasing sample size improves the power to detect mediation effect and power to detect the direct effect (of course). In particular, sample size of 150 seems to be the best sample size to select (though it is depndent on the specific parameters we have selected). Selecting the exact power we want is somewhat arbitrary - for some 80% power may be sufficient in light of realistic constraints such as sample availability. But the output alone seems to suggest that 150 is sweetspot to approximately get 100% power with an alpha of .01.\n",
        ">"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RcCo3QOifUi5"
      },
      "source": [
        "**Given** how we generated the data, why was the direct effect harder to detect than the mediated effect?\n",
        "> The direct effects are harder to detect than the mediated effect because when analyzing the direct effect, part of the variance it explains in reading comprehension is already accounted for by the mediator (word reading experience)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C4MPECMmZXqe"
      },
      "source": [
        "**DUE:** 5pm EST, April 3, 2024"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r9GUofXN4BVy"
      },
      "source": [
        "**IMPORTANT** Did you collaborate with anyone on this assignment? If so, list their names here.\n",
        "> *Someone's Name*"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "R",
      "language": "R",
      "name": "ir"
    },
    "language_info": {
      "codemirror_mode": "r",
      "file_extension": ".r",
      "mimetype": "text/x-r-source",
      "name": "R",
      "pygments_lexer": "r",
      "version": "4.2.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}