if (!"pacman" %in% rownames(installed.packages())) install.packages("pacman")
library(pacman)
pacman::p_load(tidyverse, pwr, compute.es)Connect with me on Open Science Framework | Contact me via LinkedIn
It might be necessary to use right-click -> open in a new browser window depending on your machine.
R analysis script presenting part of the solutions for exercise #1 in Sensory Marketing and Product Innovation regarding Li, Heuvinck, & Pandelaere (2021). In one of the questions, students are asked to estimate the appropriate sample size to replicate Study 1’s findings (under \(\alpha\)=5%, Power=0.8).
The purpose of this script does not solely lay in answering the exercise question. Moreover, studying these scripts should students make become familiar with some aspects of working in R.
If this exercise grabs your attention, please check out our master study programs at the Otto-von-Guericke-University in Magdeburg (Germany) by clicking on the logo:
1 Loading packages
R is a context-sensitive language. Thus, ‘data’ will be interpreted not in the same way as ‘Data’ will.
In R most functionality is provided by additional packages.
Most of the packages are well-documented, see: https://cran.r-project.org/
The code chunk below first evaluates if the package pacman (Rinker & Kurkiewicz, 2018) is already installed on your machine. If yes, the corresponding package will be loaded. If not, R will install the package.
Alternatively, you can do this manually first by executing install.packages(“pacman”) and then library(pacman).
The second line then loads the package pacman.
The third line uses the function p_load() from the pacman package to install (if necessary) and loads all packages that we provide as arguments (e.g., pwr (Champely, 2020), which provides functions for statistical power calculations).
In all code chunks throughout this script, you can receive additional help on each used function by clicking on its name (or via right-click and then opening in a new browser tab). Alternatively, when coding, we can see which arguments a function understands by pressing ‘F1’ while setting the cursor to the function’s name.
Here is the R session info which gives you information on my machine, all loaded packages, and their version:
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] compute.es_0.2-5 pwr_1.3-0 htmltools_0.5.6 quarto_1.3
[5] ggpubr_0.6.0 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.0
[9] dplyr_1.1.3 purrr_1.0.2 readr_2.1.4 tidyr_1.3.0
[13] tibble_3.2.1 ggplot2_3.4.3 tidyverse_2.0.0 labelled_2.12.0
[17] knitr_1.44 kableExtra_1.3.4 haven_2.5.3 pacman_0.5.1
loaded via a namespace (and not attached):
[1] gtable_0.3.4 xfun_0.40 processx_3.8.2 rstatix_0.7.2
[5] tzdb_0.4.0 vctrs_0.6.3 tools_4.3.1 ps_1.7.5
[9] generics_0.1.3 fansi_1.0.4 pkgconfig_2.0.3 webshot_0.5.5
[13] lifecycle_1.0.3 compiler_4.3.1 munsell_0.5.0 carData_3.0-5
[17] yaml_2.3.7 pillar_1.9.0 later_1.3.1 car_3.1-2
[21] abind_1.4-5 tidyselect_1.2.0 rvest_1.0.3 digest_0.6.33
[25] stringi_1.7.12 fastmap_1.1.1 grid_4.3.1 colorspace_2.1-0
[29] cli_3.6.1 magrittr_2.0.3 utf8_1.2.3 broom_1.0.5
[33] withr_2.5.1 scales_1.2.1 backports_1.4.1 timechange_0.2.0
[37] rmarkdown_2.25 httr_1.4.7 ggsignif_0.6.4 hms_1.1.3
[41] evaluate_0.22 viridisLite_0.4.2 rlang_1.1.1 Rcpp_1.0.11
[45] glue_1.6.2 xml2_1.3.5 svglite_2.1.1 rstudioapi_0.15.0
[49] jsonlite_1.8.7 R6_2.5.1 systemfonts_1.0.4
2 Finding an answer to question #4
For your master thesis, it could be an idea to replicate the IAT findings of Study 1 (D-score > 0, t-test). Based on a power analysis (Assumptions: \(\alpha\)=5%, Power=0.8), how many subjects should you recruit for such a study?
2.1 Extract information on the effect size Cohen’s d and the statistical analysis procedure
Information the article provides on page 3:
The D-score was significantly greater than 0 (D = 0.45, SD = 0.47, t(133) = 10.98, p < .001, 95% CI = [0.37, 0.53], Cohen’s d = 0.95), suggesting that participants responded significantly faster in assigning the stimuli when healthy and light were grouped together.
From this, we are able to extract:
- The effect size measure Cohen’s d equals 0.95.
- For analysis, the authors applied a one-sample t-test, which tests against the constant value of 0.
2.2 Conduct an a priori power analysis to search for minimum n
In the question, we are asked for a statistical power of 80% (a commonly used threshold in social sciences).
To estimate a minimum sample size n for the replication, we apply the pwr.t.test() function from the pwr package. This function needs us to provide 6 arguments to fill its parameters. These are (see help by pressing ‘F1’ when setting the cursor into the function’s name):
- d = Effect size Cohen’s d (see, e.g., Cohen (2013))
- n = Total number of observations
- sig.level = Significance level (Type I error probability \(\alpha\))
- power = Power of test (1 minus Type II error probability)
- type = Type of t-test: one- (“one.sample”) two- (“two.sample”) or paired-samples (“paired”)
- alternative = parameter indicating the test strategy: one of “two.sided” (default), “greater” or “less”
The function, furthermore, assumes us to set one of the arguments (d, n, sig.level, or power) to NULL. By doing so, we tell the function to use the remaining 5 parameters to search for the value of the sixth. In our case, we are searching for ‘n’, therefore, we set ‘n=NULL’. Furthermore, we are using a one sample t-test and the original study expects the d score to be greater than 0. Therefore, we set alternative to “greater”.
We assign the results of our power analysis to a new object named ‘results’. Then we call for its content.
results <- pwr.t.test(d = 0.95, n = NULL, type = "one.sample", sig.level = 0.05, power = 0.8, alternative = "greater")
results One-sample t test power calculation
n = 8.379928
d = 0.95
sig.level = 0.05
power = 0.8
alternative = greater
We can see that the (conservatively rounded) calculated sample size 9 perfectly mirrors the one we have seen in G*Power 3 (Faul, Erdfelder, Lang, & Buchner, 2007). As discussed in our exercise, the authors of the study applied standard exclusion criteria for IAT procedures. Therefore, 1.49 times the sample size estimate (1.49 times 9 = 14) would be necessary to end at n = 9 in the net sample.
Keep in mind: if an original study comes with a very large effect size, then a replication with a power of 80% often needs fewer subjects as compared to the original study. However, studies reporting only small effects usually need many more participants for a successful replication.
3 Additional reading (not relevant for the exam)
In the above-discussed example, we were planning a replication of the original findings. In doing so, we have made certain assumptions:
We assume that the observed effect size Cohen’s d = 0.95 is the true effect size to find. Strictly speaking, this is wrong. Our calculated value is only a sample estimate of the true effect size in the underlying population. Under the assumption of knowing the true effect size, we calculated that we need at least n=14 participants for a replication of the original findings with a power of 80%. This means that we will have approximately the chance of 80% to find the effect to be significant at \(\alpha\)=0.05. Because, the true effect to find might, however, be smaller, this calculation may present an overly optimistic approach. Therefore, other researchers prefer to focus on the confidence intervals for the estimate of effect size measures (Thompson, 2002).
To estimate the lower limit of the observed Cohen’s d, we can apply a formula proposed by Hedges & Olkin (1985):
\(d \pm se * z_{crit}\)
Where se is the standard error given by:
\(se = \sqrt{\frac{1}{n}+\frac{d^{2}}{2n}}\)
and \(z_{crit}\) is given by the standard normal distribution at p = 0.95 (\(\alpha\) = 5%).
Let us first calculate se. Table 4 in the article informs that n = 134 participants were used for analysis.
se <- sqrt((1/134)+((0.95*0.95)/(2*134)))
se[1] 0.1040684
Next, let us determine the critical z-value (\(\alpha\) = 5%). For this purpose, we use the qnorm() function. We only have to provide the probability (95%) as well as the mean and standard deviation.
z_crit <- qnorm(p = 0.95, mean = 0, sd = 1, lower.tail = TRUE)
z_crit[1] 1.644854
Finally, we calculate the lower limit of the 95% confidence interval of Cohen’s d:
d_lower95 <- 0.95 -se*z_crit
d_lower95[1] 0.7788228
The results tell us which lower Cohen’s d to expect in the underlying population (0.78). This is an alternative to the classical t-test. If we reach a significant result here (i.e., the lower limit is not crossing 0), this means that the effect size obtained is significantly greater than 0. Put differently, the effect is assumed to exist in the underlying population. Researchers in the context effect domain are starting to report tests on effect sizes, especially, when their articles span over multiple studies (Evangelidis, Levav, & Simonson, 2018).
We can briefly calculate how this will change our assumptions about the minimum sample size. We, again, use the pwr.t.test() function (see above). This time, we hand in d_lower95 as the argument for the effect size parameter d.
results_worst_case <- pwr.t.test(d = d_lower95, n = NULL, type = "one.sample", sig.level = 0.05, power = 0.8, alternative = "greater")
results_worst_case One-sample t test power calculation
n = 11.67058
d = 0.7788228
sig.level = 0.05
power = 0.8
alternative = greater
We can see that the more conservative sample size is 12. Combined with the exclusion factor discussed above, we end at a (worst-case) gross sample n of 18.