pwrcode_cb.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "![Introduction to Power Analysis with R - AcqVA Aurora workshop](https://slcladal.github.io/images/acqvalab.png)\n",
                "\n",
                "\n",
                "# Introduction\n",
                "\n",
                "This tutorial introduces power analysis for determining sample size using R with the aim of informing about the theoretical underpinnings, showcasing how to perform power analysis, and enabling participants on how to adapt, change, and modify data sets so that they can be fed into a power analysis using the `simr` (Green and MacLeod 2016a). \n",
                "\n",
                "\n",
                "\n",
                "This tutorial is aimed at intermediate users of R. The aim is not to provide a fully-fledged guide but rather to show and exemplify how to perform a power analysis for a generalized linear mixed-effects model as this is the most common type of analysis for experimental studies.\n",
                "\n",
                "\n",
                "<div class=\"warning\" style='padding:0.1em; background-color:#f2f2f2; color:#51247a'>\n",
                "<span>\n",
                "<p style='margin-top:1em; text-align:center'>\n",
                "The entire R Notebook for the tutorial can be downloaded  via https://github.com/MartinSchweinberger/AcqVA_PowerR_WS/acqvapwrws.Rmd.  If you want to render the R Notebook on your machine, i.e. knitting the document to html or a pdf, you need to make sure that you have R and RStudio installed and you also need to download the bibliography file (see https://github.com/MartinSchweinberger/AcqVA_DataVisR_WS/bibliography.bib) and store it in the same folder where you store the Rproj file.<br><\/p>\n",
                "<p style='margin-left:1em;'>\n",
                "<\/p><\/span>\n",
                "<\/div>\n",
                "\n",
                "<br>\n",
                "\n",
                "\n",
                "\n",
                "## What to do before the workshop\n",
                "\n",
                "To get the most out of this workshop, you will need to have some (basic) R skills and (basic) knowledge of how to work with R, RStudio, R Projects, and R Notebooks. If you have no or little experience with this or if you need to refresh your skills, please carefully read (or optimally go through) these tutorials:\n",
                "\n",
                "* [Getting started with R](https://slcladal.github.io/intror.html)\n",
                "\n",
                "\n",
                "Before attending the workshop, you need to install the following packages in RStudio:\n",
                "\n",
                "* `here` (for easy pathing)\n",
                "\n",
                "* `tidyverse` (for data processing)\n",
                "\n",
                "* `pwr` (for simple power analyses)\n",
                "\n",
                "* `lme4` (for mixed-models)\n",
                "\n",
                "* `sjPlot` (for summaries of mixed-models)\n",
                "\n",
                "* `report` (for reporting mixed-models)\n",
                "\n",
                "* `simr` (for power analyses for glmms)\n",
                "\n",
                "* `knitr` (for knitting R Notebooks)\n",
                "\n",
                "* `markdown` (for rendering Rmd files)\n",
                "\n",
                "* `rmarkdown`  (for R Markdown formatting)\n",
                "\n",
                "* `installr` (for updating R)\n",
                "\n",
                "\n",
                "You can update R and install these packages by clicking on `packages` on the top of the lower right pane in RStudio or you can execute the following code in the lower left, `Console` pane of RStudio. \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# update R\n",
                "#install.packages(\"installr\")\n",
                "#library(installr)\n",
                "#updateR()\n",
                "# install required packages\n",
                "install.packages(\"devtools\")\n",
                "library(\"devtools\")\n",
                "install_version(\"nloptr\", version = \"1.2.1\", repos = \"http://cran.us.r-project.org\")\n",
                "install.packages(c(\"tidyverse\", \"pwr\", \"sjPlot\", \"lme4\", \"simr\"))\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**It is really important that you have knowledge of R and RStudio and that you have installed the packages before the workshop so that we do not have to deal with technical issues too much.**\n",
                "\n",
                "You can find a more detailed version of the content of this workshop in these two LADAL tutorials\n",
                "\n",
                "* [Power Analysis with R](https://slcladal.github.io/pwr.html)\n",
                "\n",
                "You can follow this workshop in different ways - you can sit back and watch it like a lecture or take a more active role - that said, the intention for this workshop is clearly to be practical so that I show something and then you do it on you computer and we have exercises where you can try out what you have just learned.\n",
                "\n",
                "So, in essence, there are the following three options for following this workshop:\n",
                "\n",
                "1. You can sit back and enjoy and focus on what you see and then go back home and try it by yourself later.\n",
                "\n",
                "2. You can follow this workshop in RStudio and execute code, see what it does, understand it, and adapt it. This requires some skills in R and RStudio - although I have trued to keep things simple.\n",
                "\n",
                "3. You can use the Jupyter Notebook on Google Colab to follow the workshop (and you are then ready to go - you only  have to install the packages which takes a couple of minutes). You can then copy that Notebook and save it in your own Google Drive and then execute code, modify it, and understand it. This does require less skills and will be easier for people who just want to focus on the code that produces certain visualizations without doing it in RStudio. Also, this get rid of most technical issues (hopefully) but the installation of packages will take a while and you will need to re-install the packages every time you want to re-use the Jupyter Notebook.\n",
                "\n",
                "Choose which option suits you best and then go with it. \n",
                "\n",
                "## Timeline\n",
                "\n",
                "Here is what we have planned to cover in this workshop:\n",
                "\n",
                "Tuesday, February 1, 10-12am\n",
                "\n",
                "* Introduction\n",
                "\n",
                "* What is statistical power?\n",
                "\n",
                "* Why perform power analyses?\n",
                "\n",
                "* Session preparation \n",
                "\n",
                "* Basic Power Analysis (`pwr`)\n",
                "\n",
                "\n",
                "Tuesday, February 1, 1-3pm\n",
                "\n",
                "* Advanced power analyses\n",
                "\n",
                "* Power analysis for piloted data (`simr`)\n",
                "\n",
                "* Linear Mixed-Effects Model\n",
                "\n",
                "* Linear Mixed-Effects Model\n",
                "\n",
                "* Post-hoc power analyses \n",
                "\n",
                "* Setting effect sizes\n",
                "\n",
                "* Wrap-up and Resources\n",
                "\n",
                "\n",
                "\n",
                "# Getting started\n",
                "\n",
                "If you choose option 2, you need to set up our R session and prepare our R project at the very beginning of the workshop. \n",
                "\n",
                "For everything to work, please do the following:\n",
                "\n",
                "* Create a folder for this workshop somewhere on your computer, e.g. called *AcqVA_DataVisR_WS*\n",
                "\n",
                "* In that folder, create two subfolders called *data* and *images*\n",
                "\n",
                "* Open RStudio, go to File > New Project > Existing Directory (Browse to project folder) > Create (and hit Enter)\n",
                "\n",
                "This will then create an R Project in the project folder.\n",
                "\n",
                "\n",
                "## Theoretical underpinning\n",
                "\n",
                "\n",
                "What is power and why do we need power analyses?\n",
                "\n",
                "Let's have a look at some distributions to understand what is going on...\n",
                "\n",
                "## Effect size\n",
                "\n",
                "And let's start with varying effect size.\n",
                "\n",
                "\n",
                "![2 groups from 1 distribution, N = 30](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist01.png)\n",
                "\n",
                "\n",
                "Now, lets have a look at at the distribution of two different groups (group has a weak effect, Cohen's *d* = .2)\n",
                "\n",
                "\n",
                "![2 groups from 2 distribution (weak effect, N = 30)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist02.png)\n",
                "\n",
                "Now, lets have a look at at the distribution of two different groups (group has a medium effect, Cohen's *d* = .5)\n",
                "\n",
                "![2 groups from 2 distribution (medium effect, N = 30)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist03.png)\n",
                "\n",
                "\n",
                "\n",
                "\n",
                "Now, lets have a look at at the distribution of two different groups (group has a strong effect, Cohen's *d* = .8)\n",
                "\n",
                "![2 groups from 2 distribution (strong effect, N = 30)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist04.png)\n",
                "\n",
                "\n",
                "\n",
                "Now, lets have a look at at the distribution of two different groups (group has a very strong effect, Cohen's *d* = 1.2)\n",
                "\n",
                "![2 groups from 2 distribution (very strong effect, N = 30)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist05.png)\n",
                "\n",
                "\n",
                "**If variability and sample size remain constant, larger effects are easier to detect than smaller effects!**\n",
                "\n",
                "## Sample size\n",
                "\n",
                "And let's now look at sample size.\n",
                "\n",
                "![2 groups from 2 distribution (medium effect, N = 30)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist06.png)\n",
                "\n",
                "\n",
                "\n",
                "\n",
                "Let us now increase the sample size to N = 50.\n",
                "\n",
                "![2 groups from 2 distribution (medium effect, N = 50)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist07.png)\n",
                "\n",
                "\n",
                "\n",
                "**If variability and effect size remain constant, effects are easier to detect with increasing sample size!**\n",
                "\n",
                "## Variability\n",
                "\n",
                "And let's now look at variability.\n",
                "\n",
                "![2 groups from 2 distribution (medium effect, N = 30, SD = 10)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist08.png)\n",
                "\n",
                "\n",
                "Let's decrease the variability to sd = 5.\n",
                "\n",
                "![2 groups from 2 distribution (medium effect, N = 30, SD = 5)](https://raw.githubusercontent.com/MartinSchweinberger/AcqVAPowerRWS/main/images/dist09.png)\n",
                "\n",
                "\n",
                "**If the sample and effect size remain constant, effects are easier to detect with decreasing variability!**\n",
                "\n",
                "**In summary**, there are three main factors that determine if a model finds an effect. The accuracy (i.e., the probability of finding an effect):\n",
                "\n",
                "* the size of the effect (bigger effects are easier to detect)\n",
                "* the variability of the effect (less variability makes it easier to detect an effect), and \n",
                "* the sample size (the bigger the sample size, the easier it is to detect an effect); \n",
                "  + number of subjects/participants\n",
                "  + number of items/questions\n",
                "  + number of observations per item within subjects/participants\n",
                "  \n",
                "Now, if a) we dealing with a very big effect, then we need only few participants and few items to accurately find this effect.\n",
                "\n",
                "Or b) if we dealing with an effect that has low variability (it is observable for all subjects with the same strength), then we need only few participants and few items to accurately find this effect.\n",
                "\n",
                "\n",
                "\n",
                "<div class=\"warning\" style='padding:0.1em; background-color:#f2f2f2; color:#51247a'>\n",
                "<span>\n",
                "<p style='margin-top:1em; text-align:center'>\n",
                "What power analysis allow us to do is to calculate at what sample size we would be able to detect an effect of a certain strength with a predefined certainty, given a certain amount of variability. <br><br>In other words, assuming the variability in the pilot data is the same as in the population, how many subjects do we need to detect a weak (or median or strong) effect with 80 percent certainty?<br><\/p>\n",
                "<p style='margin-left:1em;'>\n",
                "<\/p><\/span>\n",
                "<\/div>\n",
                "\n",
                "<br>\n",
                "\n",
                "# Today's Data\n",
                "\n",
                "\n",
                "We will use 1 data set that is based on a two data set that was send in by Isabel. \n",
                "\n",
                "\n",
                "Here are summaries of the data sets:\n",
                "\n",
                "1. **`AJT_V2`** (we will call this data `regdat`: contains 480 observations of the following 5 variables:\n",
                "  + ID: 20 levels (0% missing)\n",
                "  + Sentence: 24 levels\n",
                "  + Group: 2 levels, namely *English* (n = 240) and *NorwegianEnglish* (n = 240)\n",
                "  + WordOrder: 2 levels, namely *V2* (n = 240) and *V3* (n = 240)\n",
                "  + SentenceType: 2 levels, namely *NonSubjectInitial* (n = 240) and *SentenceAdverbials* (n = 240)\n",
                "\n",
                "\n",
                "**To upload the data into Google Colab, do the following**\n",
                "\n",
                "To be able to load data into Google Colab, you need to click on the folder symbol to the left of the screen:\n",
                "\n",
                "![Colab Folder Symbol](https://slcladal.github.io/images/ColabFolder.png)\n",
                "\n",
                "Then on the upload symbol. \n",
                "\n",
                "![Colab Upload Symbol](https://slcladal.github.io/images/ColabUpload.png)\n",
                "\n",
                "\n",
                "Next, upload the files you want to use. When you then execute the code (like to code chunk below, you will load the data into this notebook.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "regdat <- readRDS(\"regdat.rda\")\n",
                "# inspect\n",
                "regdat\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "To apply the code and functions below to your own data, you will need to modify the code chunks and replace the data we use here with your own data object. \n",
                "\n",
                "# Basic Power Analysis\n",
                "\n",
                "Let's start with a simple power analysis to see how power analyses work for simpler or basic statistical tests such as t-test, $\\chi$^2^-test, or linear regression.\n",
                "\n",
                "The `pwr` package  (Champely 2020) implements power analysis as outlined by Cohen (1988) and  allows to perform power analyses for the following tests (selection):\n",
                "\n",
                "* balanced one way ANOVA (`pwr.anova.test`) \n",
                "\n",
                "* chi-square test (`pwr.chisq.test`)\n",
                "\n",
                "* Correlation (`pwr.r.test`)\n",
                "\n",
                "* general linear model (`pwr.f2.test`)\n",
                " \t\n",
                "* paired (one and two sample) t-test (`pwr.t.test`)\n",
                "\n",
                "* two sample t-test with unequal N (`pwr.t2n.test`)\n",
                "\n",
                "For each of these functions, you enter three of the four parameters (effect size, sample size, significance level, power) and the fourth is calculated. \n",
                "\n",
                "So how does this work? Have a look at the code chunk below.\n",
                "\n",
                "## Power Analysis for One-Way ANOVAs\n",
                "\n",
                "Let check how to calculate the necessary sample size for each group for a one-way ANOVA that compares 5 groups and that has a power of 0.80 (80 percent), when the effect size is moderate (*f* = 0.25) and the significance level is 0.05 (5 percent). Regarding the effect size *f*, *f* has been categorized as small = .1, medium .25, and large $≥$ 0.4 (see Cohen 1988).\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# load package\n",
                "library(pwr)\n",
                "# calculate minimal sample size\n",
                "pwr.anova.test(k=5,            # 5 groups are compared\n",
                "               f=.25,          # moderate effect size\n",
                "               sig.level=.05,  # alpha/sig. level = .05\n",
                "               power=.8)       # confint./power = .8\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "In this case, the minimum number of participants in each group would be 40.\n",
                "\n",
                "Let's check how we could calculate the power if we had already collected data (with 30 participants in each group) and we want to report the power of our analysis (and et us assume that the effect size was medium).\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# calculate minimal sample size\n",
                "pwr.anova.test(k=5,            # 5 groups are compared\n",
                "               f=.25,          # moderate effect size\n",
                "               sig.level=.05,  # alpha/sig. level = .05\n",
                "               n=30)           # n of participants\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "In this case our analysis would only have a power or 66.8 percent. This means that we would only detect a medium sized effect in 66.7 percent of cases (which is considered insufficient).\n",
                "\n",
                "\n",
                "## Power Analysis for GLMs\n",
                "\n",
                "When determining the power of a generalized linear model, we need to provide\n",
                "\n",
                "* the degrees of freedom for numerator (´u´)\n",
                "\n",
                "* the degrees of freedom for denominator (`v`)\n",
                "\n",
                "* the effect size (the estimate of the intercept and the slope/estimates of the predictors)\n",
                "\n",
                "* the level of significance (i.e.,  the type I error probability)\n",
                "\n",
                "The values of the parameters in the example below are adapted from the fixed-effects regression example that was used to analyze different teaching styles (see [here](https://slcladal.github.io/regression.html#Example_2:_Teaching_Styles)).\n",
                "\n",
                "The effect size used here is $f^2$ that has be categorized as follows (see Cohen 1988): small $≥$ 0.02, medium $≥$ 0.15, and large $≥$ 0.35. So in order to determine if the data is sufficient to find a weak effect when comparing 2 groups with 30 participants in both groups (df_numerator_: 2-1; df_denominator_: (30-1) + (30-1)) and a significance level at $\\alpha$ = .05, we can use the following code.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# general linear model\n",
                "pwrglm <- pwr.f2.test(u = 1,\n",
                "                      v = 58,\n",
                "                      f2 = .02,\n",
                "                      sig.level = 0.05)\n",
                "# inspect results\n",
                "pwrglm\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "The results show that the regression analyses used to evaluate the effectiveness of different teaching styles only had a power of `r pwrglm$power`.\n",
                "\n",
                "\n",
                "## Power Analysis for t-tests\n",
                "\n",
                "For t-tests (both paired and 2-sample t-tests), the effect size measure is Cohen's $d$ that has be categorized as follows (see Cohen 1988): small $≥$ 0.2, medium $≥$ 0.5, and large $≥$ 0.8. \n",
                "\n",
                "**Paired t-test**\n",
                "\n",
                "So in order to determine if the data is sufficient to find a weak effect when comparing the pre- and post-test results of a group with 30 participants, evaluating an undirected hypothesis (thus the *two-tailed* approach), and a significance level at $\\alpha$ = .05, we can use the following code.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# paired t-test\n",
                "pwrpt <- pwr.t.test(d=0.2,\n",
                "                    n=30,\n",
                "                    sig.level=0.05,\n",
                "                    type=\"paired\",\n",
                "                    alternative=\"two.sided\")\n",
                "# inspect\n",
                "pwrpt\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Given the data, a weak effect in this design can only be detected with a certainty of `r pwrpt$power` percet. This means that we would need to substantively increase the sample size to detect a small effect with this design.\n",
                "\n",
                "**Two-sample (independent) t-test**\n",
                "\n",
                "The power in a similar scenario but with two different groups (with 25 and 35 subjects) can be determined as follows (in this case we test a directed hypothesis that checks if the intervention leads to an increase in the outcome - hence the *greater* in the `alternative` argument):\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# independent t-test\n",
                "pwr2t <- pwr.t2n.test(d=0.2,\n",
                "                      n1=35,\n",
                "                      n2 = 25,\n",
                "                      sig.level=0.05,\n",
                "                      alternative=\"greater\")\n",
                "# inspect\n",
                "pwr2t\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Given the data, a weak effect in this design can only be detected with a certainty of `r pwr2t$power` percent. This means that we would need to substantively increase the sample size to detect a small effect with this design.\n",
                "\n",
                "\n",
                "## Power Analysis for $\\chi$^2^-tests\n",
                "\n",
                "Let us now check the power of a $\\chi^2$^-test. For $\\chi^2$^-test, the effect size measure used in the power analysis is $w$ that has be categorized as follows (see Cohen 1988): small $≥$ 0.1, medium $≥$ 0.3, and large $≥$ 0.5. Also, keep in mind that for $\\chi^2$^-tests, at least 80 percent of cells need to have values $≥$ 5 and none of the cells should have expected values smaller than 1 (see Bortz, Lienert, and Boehnke 1990).. \n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# x2-test\n",
                "pwrx2 <- pwr.chisq.test(w=0.2,\n",
                "                        N = 25, # total number of observations\n",
                "                        df = 1,\n",
                "                        sig.level=0.05)\n",
                "# inspect\n",
                "pwrx2\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Given the data, a weak effect in this design can only be detected with a certainty of `r pwrx2$power` percent. This means that we would need to substantively increase the sample size to detect a small effect with this design.\n",
                "\n",
                "**EXERCISE**\n",
                "\n",
                "For the tests above (anova, glm, paired and independent t-test, and the $\\chi^2$-test), how many participants would you need to have a power of 80 percent?\n",
                "\n",
                "Here are the commands we used to help you:\n",
                "\n",
                "* ANOVA: `pwr.anova.test(k=5, f=.25, sig.level=.05, n=30) `\n",
                "\n",
                "* GLM: `pwr.f2.test(u = 1, v = 58, f2 = .02, sig.level = 0.05)`\n",
                "\n",
                "* paired t-test: `pwr.t.test(d=0.2, n=30, sig.level=0.05, type=\"paired\", alternative=\"two.sided\")`\n",
                "\n",
                "* independent t-test: `pwr.t2n.test(d=0.2,n1=35, n2 = 25, sig.level=0.05, alternative=\"greater\")`\n",
                "\n",
                "* $\\chi^2$-test: `pwr.chisq.test(w=0.2, N = 25,  df = 1, sig.level=0.05)`\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Advanced Power Analysis\n",
                "\n",
                "The basis for the present tutorial is Green and MacLeod (2016b) (which you can find [here](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504)). Green and MacLeod (2016b) is a highly recommendable and thorough tutorial on performing power analysis in R. Recommendable literature on this topic are, e.g. Arnold et al. (2011)  and Johnson et al. (2015)  and [this tutorial](https://www.journalofcognition.org/articles/10.5334/joc.10/). \n",
                "\n",
                "\n",
                "\n",
                "Before we conduct a study, we should figure out, what sample we need to detect a small/medium effect with medium variability so that our model is sufficient to detect this kind of effect. In order to do this, we would generate a data set that mirrors the kind of data that we expect to get (with the properties that we expect to get). We can then fit a model to this data and check if a model would be able to detect the expected effect. However, because a single model does not tell us that much (ift could simply be luck that it happened to find the effect), we run many different models on variations of the data and see how many of them find the effect. As a general rule of thumb, we want a data set that allows a model to find a medium sized effect with at least an accuracy of 80 percent (Field et al. 2007).\n",
                "\n",
                "In the following, we will go through how to determine what sample size we need for an example analysis.\n",
                "\n",
                "\n",
                "## Preparation and session set up\n",
                "\n",
                "As we we have installed the packages, we can activate them as shown below.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# load packages\n",
                "library(tidyverse)\n",
                "library(lme4)\n",
                "library(sjPlot)\n",
                "library(simr)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Once you have installed R and RStudio and initiated the session by executing the code shown above, you are good to go.\n",
                "\n",
                "## Using piloted data and a lmer\n",
                "\n",
                "Load data\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# load data\n",
                "regdat <- readRDS(\"regdat.rda\")\n",
                "# inspect\n",
                "head(regdat, 10);str(regdat)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Check how many levels we have.\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "table(regdat$Group, regdat$WordOrder); \n",
                "table(regdat$Group, regdat$SentenceType); \n",
                "table(regdat$ID)\n",
                "table(regdat$Sentence)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "We could also add a response variable (but not now ;) )\n",
                "\n",
                "\n",
                "### Generating the model\n",
                "\n",
                "Set effects and specify model parameters\n",
                "\n",
                "* fixed effects (and intercept)\n",
                "\n",
                "* random effects\n",
                "\n",
                "* variance\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Intercept + slopes for fixed effects \n",
                "# (Group, Sentence_type, Word_order, and an interaction between Group * Sentence_type)\n",
                "fixed <- c(.52, .52, .52, .52, .52)\n",
                "# Random intercepts for Sentence and ID\n",
                "#rand <- list(0.5, 0.1)\n",
                "rand <- list(0.5)\n",
                "# res. variance\n",
                "res <- 2\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Generate model\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "#model <- makeLmer(y ~ (1|Sentence) + (1|ID) + Group * Sentence_type + Word_order, fixef=fixed, VarCorr=rand, sigma=res, data=regdat)\n",
                "m1 <- makeLmer(y ~ (1|Sentence) + Group * SentenceType + WordOrder, \n",
                "               fixef=fixed, \n",
                "               VarCorr=rand, \n",
                "               sigma=res, \n",
                "               data=regdat)\n",
                "# inspect\n",
                "m1\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Inspect summary table\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "sjPlot::tab_model(m1)\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**The summary table shows that the effects are correct but none of them are reported as being significant!**\n",
                "\n",
                "### Power analysis\n",
                "\n",
                "Let us now check if the data is sufficient to detect the interaction between `Group` and `Sentence_type`.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "sim_wo <- simr::powerSim(m1, \n",
                "                        nsim=10,\n",
                "                        test = fcompare(y ~ WordOrder))\n",
                "# inspect results\n",
                "sim_wo\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**The data is sufficient to detect a weak effect of `WordOrder` with `r sim_wo$x` percent accuracy**\n",
                "\n",
                "**Why did I set the estimates to 0.52?**\n",
                "\n",
                "Let us now inspect the effect sizes and see why I used 0.52 as an effect size in the model.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# extract fixed effect estimates\n",
                "estimatesfixedeffects <- fixef(m1)\n",
                "# convert estimates into odds ratios\n",
                "exp(estimatesfixedeffects)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "We will now check the effect size which can be interpreted according to Chen, Cohen, and Chen (2010; see also Cohen 1988).\n",
                "\n",
                "* small effect (Cohen's d 0.2, OR = 1.68)\n",
                "\n",
                "* medium effect (Cohen's d 0.5, OR = 3.47)\n",
                "\n",
                "* strong effect  (Cohen's d 0.8, OR = 6.71)\n",
                "\n",
                "### Power analysis for interactions\n",
                "\n",
                "Power Analysis: let's check if the data set has enough power to detect a weak effect for `Group * SentenceType`.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "sim_gst <- simr::powerSim(m1, \n",
                "                          nsim=100, \n",
                "                          test = fcompare(y ~ Group * SentenceType))\n",
                "# inspect\n",
                "sim_gst\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**The data is not sufficient to detect a weak effect of `Group * SentenceType` with `r sim_gst$x` percent accuracy**\n",
                "\n",
                "\n",
                "## Extending Data\n",
                "\n",
                "We will now extend the data to see how many participants and sentences you need.\n",
                "\n",
                "\n",
                "### Adding sentences\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m1_as <- simr::extend(m1,\n",
                "                      along=\"Sentence\", \n",
                "                      n=50)\n",
                "# inspect model\n",
                "m1_as\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Check power when using 60 sentences\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "sim_m1_as_gst <- powerSim(m1_as, \n",
                "                   nsim=20, \n",
                "                   test = fcompare(y ~ Group * SentenceType))\n",
                "# inspect\n",
                "sim_m1_as_gst\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Plot power curve\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "pcurve_m1_as_gst <- simr::powerCurve(m1_as, \n",
                "                               test=fcompare(y ~ Group * SentenceType), \n",
                "                               along=\"Sentence\",\n",
                "                               nsim = 20, \n",
                "                               breaks=seq(10, 50, 10))\n",
                "# show plot\n",
                "plot(pcurve_m1_as_gst)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Checking participants\n",
                "\n",
                "What about increasing the number of participants?\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m2 <- makeLmer(y ~ (1|ID) + Group * SentenceType + WordOrder, \n",
                "               fixef=fixed, \n",
                "               VarCorr=rand, \n",
                "               sigma=res, \n",
                "               data=regdat)\n",
                "# inspect\n",
                "sjPlot::tab_model(m2)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Power analysis\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "sim_m2_wo <- powerSim(m2, \n",
                "                   nsim=10, \n",
                "                   test = fcompare(y ~ WordOrder))\n",
                "# inspect\n",
                "sim_m2_wo\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "An alternative way to perform the simulation\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# set seed for replicability\n",
                "set.seed(12345)\n",
                "# perform power analysis for present model\n",
                "rsim_m2_wo <- powerSim(m2, fixed(\"WordOrderV3\", \"z\"), nsim=10)\n",
                "# inspect\n",
                "rsim_m2_wo\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Plotting the power curve.\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "pcurve_m2_wo <- simr::powerCurve(m2, \n",
                "                               test=fcompare(y ~ WordOrder),\n",
                "                               along = \"ID\", \n",
                "                               nsim=10)\n",
                "# show plot\n",
                "plot(pcurve_m2_wo)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Generating data and a glmer\n",
                "\n",
                "In order to perform a  power analysis, we will start by loading the tidyverse package to process the data and by generating a data that we will use to determine the power of a regression model.\n",
                "\n",
                "This simulated data set has\n",
                "\n",
                "* 200 data points\n",
                "* 2 Conditions (Control, Test)\n",
                "* 10 Subjects\n",
                "* 10 Items\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# generate data\n",
                "simdat <- data.frame(\n",
                "  sub <- rep(paste0(\"Sub\", 1:10), each = 20),\n",
                "  cond <- rep(c(\n",
                "    rep(\"Control\", 10),\n",
                "    rep(\"Test\", 10))\n",
                "    , 10),\n",
                "  itm <- as.character(rep(1:10, 20))\n",
                ") %>%\n",
                "  dplyr::rename(Subject = 1,\n",
                "                Condition = 2,\n",
                "                Item = 3) %>%\n",
                "  dplyr::mutate_if(is.character, factor)\n",
                "# inspect\n",
                "head(simdat, 15)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**EXERCISE**\n",
                "\n",
                "Can you create a data set with 5 Subjects, 5 Items (each) for 4 Conditions?\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "**Solution**\n",
                "\n",
                "The first 20 rows of the data should look like shown below:\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# generate data\n",
                "Items <- rep(rep(paste0(\"Item\", 1:5), 4), 5)\n",
                "Condition <- rep(rep(paste0(\"Condition\", 1:4), each = 5), 5)\n",
                "Subjects <- rep(paste0(\"Subject\", 1:5), each = 20)\n",
                "# combine into data frame\n",
                "exdat <- data.frame(Subjects, Items, Condition)\n",
                "# inspect\n",
                "head(exdat, 20)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Generating the model\n",
                "\n",
                "Set effects and specify model parameters\n",
                "\n",
                "* fixed effects (and intercept)\n",
                "\n",
                "* random effects\n",
                "\n",
                "* variance\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Intercept + slopes for fixed effects \n",
                "# (Group, Sentence_type, Word_order, and an interaction between Group * Sentence_type)\n",
                "fixed <- c(.52, .52)\n",
                "# Random intercepts for Sentence and ID\n",
                "rand <- list(0.5, 0.1)\n",
                "#rand <- list(0.5)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Generate model\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m4 <- simr::makeGlmer(y ~ (1|Subject) + (1|Item) + Condition, \n",
                "                      family = \"binomial\", \n",
                "                      fixef=fixed, \n",
                "                      VarCorr=rand, \n",
                "                      data=simdat)\n",
                "# inspect\n",
                "sjPlot::tab_model(m4)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Extract power\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# set seed for replicability\n",
                "set.seed(12345)\n",
                "# perform power analysis for present model\n",
                "rsim_m4_c <- powerSim(m4, fixed(\"ConditionTest\", \"z\"), nsim=20)\n",
                "# inspect\n",
                "rsim_m4_c\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Extending the data\n",
                "\n",
                "### Adding Items\n",
                "\n",
                "* 2 Conditions\n",
                "\n",
                "* 10 Subjects\n",
                "\n",
                "* 40 Items (from 10)\n",
                "\n",
                "(each combination occurs once!)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m4_ai <- simr::extend(m4,\n",
                "                      along=\"Item\", \n",
                "                      n=40)\n",
                "# inspect model\n",
                "sjPlot::tab_model(m4_ai)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Plot power curve\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "pcurve_m4_ai_c <- powerCurve(m4_ai, \n",
                "                             fixed(\"ConditionTest\", \"z\"), \n",
                "                             along = \"Item\",\n",
                "                             nsim = 20,\n",
                "                             breaks=seq(10, 40, 5), \n",
                "                             nsim=10)\n",
                "# show plot\n",
                "plot(pcurve_m4_ai_c)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Adding subjects\n",
                "\n",
                "* 2 Conditions\n",
                "\n",
                "* 10 Items\n",
                "\n",
                "* 40 Subjects (from 10)\n",
                "\n",
                "(each combination occurs once!)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m4_as <- simr::extend(m4,\n",
                "                      along=\"Subject\", \n",
                "                      n=40)\n",
                "# inspect model\n",
                "sjPlot::tab_model(m4_as)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Plot power curve\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "pcurve_m4_as_c <- powerCurve(m4_as, \n",
                "                             fixed(\"ConditionTest\", \"z\"), \n",
                "                             along = \"Subject\",\n",
                "                             nsim = 20,\n",
                "                             breaks=seq(10, 40, 5), \n",
                "                             nsim=10)\n",
                "# show plot\n",
                "plot(pcurve_m4_as_c)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "How often does each *combination* occur?\n",
                "\n",
                "**Only once!**\n",
                "\n",
                "So, what if we increase the number of combinations? (**Repeated Measures Design!**)\n",
                "\n",
                "* 2 Conditions\n",
                "\n",
                "* 10 Items\n",
                "\n",
                "* 10 Subjects \n",
                "\n",
                "(But each combination of item and subject now occurs 10 times!)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m4_asi_c <- extend(m4, \n",
                "                   within=\"Subject+Item\", \n",
                "                   n=10)\n",
                "# perform power calculation\n",
                "pcurve_m4_asi_c <- powerCurve(m4_asi_c, \n",
                "                              fixed(\"ConditionTest\", \"z\"), \n",
                "                              within=\"Subject+Item\",\n",
                "                              nsim = 20,\n",
                "                              breaks=seq(2, 10, 2), \n",
                "                              nsim=10)\n",
                "# show plot\n",
                "plot(pcurve_m4_asi_c)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Adding subjects and items\n",
                "\n",
                "* 2 Conditions\n",
                "\n",
                "* 30 Items (from 10)\n",
                "\n",
                "* 30 Subjects (from 10)\n",
                "\n",
                "(each combination occurs once!)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "m4_as <- simr::extend(m4,\n",
                "                      along=\"Subject\", \n",
                "                      n=30)\n",
                "m4_asi <- simr::extend(m4_as,\n",
                "                      along=\"Item\", \n",
                "                      n=30)\n",
                "# inspect model\n",
                "sjPlot::tab_model(m4_asi)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Plot power curve\n",
                "\n",
                "*How many subjects do you need if you have 30 Items?*\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "pcurve_m4_asi_c <- powerCurve(m4_asi, \n",
                "                             fixed(\"ConditionTest\", \"z\"), \n",
                "                             along = \"Subject\", \n",
                "                             breaks=seq(5, 30, 5), \n",
                "                             nsim=10)\n",
                "# show plot\n",
                "plot(pcurve_m4_asi_c)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Post-Hoc Analyses\n",
                "\n",
                "We now turn to post-hoc power analyses.\n",
                "\n",
                "<div class=\"warning\" style='padding:0.1em; background-color:#f2f2f2; color:#51247a'>\n",
                "<span>\n",
                "<p style='margin-top:1em; text-align:center'>\n",
                "<b>NOTE<\/b><br>Power analysis have also been used post-hoc to test if the sample size of studies was sufficient to detect meaningful effects. However, such post-hoc power calculations where the target effect size comes from the data, give misleading results (Hoenig and Heisey 2001; Perugini, Gallucci, and Costantini 2018)  and should thus be treated with extreme care!<\/p>\n",
                "<p style='margin-left:1em;'>\n",
                "<\/p><\/span>\n",
                "<\/div>\n",
                "\n",
                "<br>\n",
                "\n",
                "We begin by adding a response variable to our data. In this case, we vary the response variable ti a higher likelihood of obtaining gazes in the *area of interests* (AOI) in the test condition.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "simdat2 <- simdat %>%\n",
                "  dplyr::arrange(Condition) %>%\n",
                "  dplyr::mutate(\n",
                "  AOI = c(sample(c(\"AOI\", \"NotAOI\"), 100, replace = T, prob = c(.5, .5)),\n",
                "           sample(c(\"AOI\", \"NotAOI\"), 100, replace = T, prob = c(.7, .3)))\n",
                "  ) %>%\n",
                "  dplyr::mutate_if(is.character, factor)\n",
                "# inspect\n",
                "head(simdat2, 20)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Now that we have generated some data, we will fit a model to it and perform a power analysis on the observed effects. \n",
                "\n",
                "\n",
                "We will fit a first model to the data. Thus, in a first step, we load the `lme4` package to create a model, set a seed (to save the results and so that the results can be replicated), and then create an initial mixed-effects model.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# set seed for replicability\n",
                "set.seed(12345)\n",
                "# fit model\n",
                "m5 <- glmer(AOI ~ (1|Subject) +(1|Item) + Condition, family=\"binomial\", data=simdat2)\n",
                "# inspect results\n",
                "summary(m5)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "We now check the effect sizes of the predictors in the model. We can do this by displaying the results of the model using the `tab_model` function from the `sjPlot` package.\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# tabulate results\n",
                "sjPlot::tab_model(m5)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Now, we perform power analysis on an observed effect. This analysis tells us how likely the model is to find an observed effect given the data.\n",
                "\n",
                "***\n",
                "\n",
                "<div class=\"warning\" style='padding:0.1em; background-color:#f2f2f2; color:#51247a'>\n",
                "<span>\n",
                "<p style='margin-top:1em; text-align:center'>\n",
                "<b>NOTE<\/b><br>We use a very low number of simulations (10) and we use the default z-test which is suboptimal for small samples (Bolker et al. 2009). In a proper study, you should increase the number of simulations (to at least 1000) and you should use a bootstrapping rather than a z-test (cf. Halekoh, Højsgaard, and others 2014). <\/p>\n",
                "<p style='margin-left:1em;'>\n",
                "<\/p><\/span>\n",
                "<\/div>\n",
                "\n",
                "***\n",
                "\n",
                "\n",
                "*What is the probability of finding the observed effect given the data?*\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# set seed for replicability\n",
                "set.seed(12345)\n",
                "# perform power analysis for present model\n",
                "rsim0 <- powerSim(m5, \n",
                "                  fixed(\"ConditionTest\", \"z\"), \n",
                "                  nsim=10)\n",
                "# inspect results\n",
                "rsim0\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "The results of the power analysis show that, given the data at hand, the model would have detected the effect of `Conidition:Test` with a probability of `r rsim0$x` percent. However, and as stated above, the results of such post-hoc power calculations (where the target effect size comes from the data) give misleading results (Hoenig and Heisey 2001)  and should thus be treated with extreme care!\n",
                "\n",
                "# Fixing effects\n",
                "\n",
                "We will set the effects that we obtained based on our \"observed\" data to check if, given the size of the data, we have enough power to detect a small effect of `Condition`. In a first step, we check the observed effects.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "estimatesfixedeffects <- fixef(m5)\n",
                "exp(estimatesfixedeffects)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "We can see that the effect is already rather small which makes it very hard to detect an effect. We will now change the size of the effect of ConditionTest to represent a truly small effect, i.e. on the brink of being noise but being just strong enough to be considered small. In other words, we will set the effect so that its odds ratio is exactly 1.68.\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# set seed for replicability\n",
                "set.seed(12345)\n",
                "# perform power analysis for small effect\n",
                "fixef(m5)[\"ConditionTest\"] <- 0.519\n",
                "estimatesfixedeffects <- fixef(m5)\n",
                "exp(estimatesfixedeffects)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "*What is the probability of finding a weak effect given the data?*\n",
                "\n",
                "Perform power analysis\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# set seed for replicability\n",
                "set.seed(12345)\n",
                "# perform power analysis for present model\n",
                "rsim0 <- powerSim(m5, \n",
                "                  fixed(\"ConditionTest\", \"z\"), \n",
                "                  nsim=20)\n",
                "# inspect results\n",
                "rsim0\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Power Analysis of Extended Data\n",
                "\n",
                "We will now extend the data to see what sample size is needed to get to the 80 percent accuracy threshold. We begin by increasing the number of items from 10 to 30 to see if this would lead to a sufficient sample size.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# increase sample size\n",
                "m5_ai <- extend(m5, \n",
                "             along=\"Item\", \n",
                "             n=30)\n",
                "# perform power simulation\n",
                "rsim2 <- powerSim(m5_ai, \n",
                "                  fixed(\"ConditionTest\", \"z\"), \n",
                "                  nsim=10)\n",
                "# show results\n",
                "rsim2\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Wrap-up and Resources\n",
                "\n",
                "If you want to know more, please have a look/go at the following resources:\n",
                "\n",
                "* [Power Analysis in R](https://slcladal.github.io/pwr.html)\n",
                "\n",
                "* [SIMR: an R package for power analysis of generalized linear mixed models by simulation](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504)\n",
                "\n",
                "* [Power analysis in Statistics with R](https://www.r-bloggers.com/2021/05/power-analysis-in-statistics-with-r/)\n",
                "\n",
                "* [Power Analysis with simr](https://humburg.github.io/Power-Analysis/simr_power_analysis.html)\n",
                "\n",
                "* [Parallelizing simr::powercurve() in R](https://www.r-bloggers.com/2021/07/parallelizing-simrpowercurve-in-r/)\n",
                "\n",
                "That's all folks!\n",
                "\n",
                "# Citation & Session Info \n",
                "\n",
                "Schweinberger, Martin. `r format(Sys.time(), '%Y')`. *Introduction to Power Analysis with R - AcqVA Aurora workshop*. Tromsø: The Arctic University of Norway, Tromsø. url: https://github.com/MartinSchweinberger/AcqVA_PowerR_WS/acqvapwrws.html.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "sessionInfo()\n",
                "\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# References \n",
                "\n",
                "Arnold, Benjamin F, Daniel R Hogan, John M Colford, and Alan E Hubbard. 2011. Simulation Methods to Estimate Design Power: An Overview for Applied Research. *British Medical Journal Medical Research Methodology* 11(1): 1–10.\n",
                "\n",
                "Bolker, Benjamin M, Mollie E Brooks, Connie J Clark, Shane W Geange, John R Poulsen, M Henry H Stevens, and Jada-Simone S White. 2009. Generalized Linear Mixed Models: A Practical Guide for Ecology and Evolution. *Trends in Ecology & Evolution* 24(3): 127–35.\n",
                "\n",
                "Bortz, J, GA Lienert, and K Boehnke. 1990. *Verteilungsfreie Methoden in der Biostatistik*. Berlin: Springer Verlag.\n",
                "\n",
                "Champely, Stephane. 2020. *Pwr: Basic Functions for Power Analysis*. https://CRAN.R-project.org/package=pwr.\n",
                "\n",
                "Chen, Henian, Patricia Cohen, and Sophie Chen. 2010. How Big Is a Big Odds Ratio? Interpreting the Magnitudes of Odds Ratios in Epidemiological Studies. *Communications in Statistics–Simulation and Computation* 39(4): 860–64.\n",
                "\n",
                "Cohen, Jacob. 1988. *Statistical Power Analysis for the Behavioral Sciences*. Hillsdale, NJ: Erlbaum.\n",
                "\n",
                "Field, Scott A, Patrick J O’Connor, Andrew J Tyre, and Hugh P Possingham. 2007. Making Monitoring Meaningful. *Australian Ecology* 32(5): 485–91.\n",
                "\n",
                "Green, Peter, and Catriona J. MacLeod. 2016a. Simr: An R Package for Power Analysis of Generalised Linear Mixed Models by Simulation. *Methods in Ecology and Evolution* 7(4): 493–98. https://doi.org/10.1111/2041-210X.12504.\n",
                "\n",
                "Green, Peter, and Catriona J MacLeod. 2016b. SIMR: An R Package for Power Analysis of Generalized Linear Mixed Models by Simulation. *Methods in Ecology and Evolution* 7(4): 493–98.\n",
                "\n",
                "Halekoh, Ulrich, Søren Højsgaard, and others. 2014. A Kenward-Roger Approximation and Parametric Bootstrap Methods for Tests in Linear Mixed Models–the R Package Pbkrtest. *Journal of Statistical Software* 59(9): 1–30.\n",
                "\n",
                "Hoenig, John M, and Dennis M Heisey. 2001. The Abuse of Power: The Pervasive Fallacy of Power Calculations for Data Analysis. *The American Statistician* 55(1): 19–24.\n",
                "\n",
                "Johnson, Paul CD, Sarah JE Barry, Heather M Ferguson, and Pie Müller. 2015. Power Analysis for Generalized Linear Mixed Models in Ecology and Evolution. *Methods in Ecology and Evolution* 6(2): 133–42.\n",
                "\n",
                "Perugini, Marco, Marcello Gallucci, and Giulio Costantini. 2018. A Practical Primer to Power Analysis for Simple Experimental Designs. *International Review of Social Psychology* 31(1): 1–23.\n"
            ]
        }
    ],
    "metadata": {
        "anaconda-cloud": "",
        "kernelspec": {
            "display_name": "R",
            "langauge": "R",
            "name": "ir"
        },
        "language_info": {
            "codemirror_mode": "r",
            "file_extension": ".r",
            "mimetype": "text/x-r-source",
            "name": "R",
            "pygments_lexer": "r",
            "version": "3.4.1"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 1
}