Connect with me on Open Science Framework | Contact me via LinkedIn
Depending on your machine it might be necessary to use right-click -> open in new browser window.
R analysis script presenting the results for exercise 3 in consumer behavior. These are the questions regarding Topolinski, Lindner, & Freudenberg (2014)
Beware: 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 is already installed to 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 load all packages that we provide as arguments (e.g., pwr, which provides functions for power analysis).
if(!"pacman" %in% rownames(installed.packages())) install.packages("pacman")
library(pacman)
pacman::p_load(tidyverse, dplyr, pwr, tidyselect, compute.es)
Here is the R session info which gives you information on my machine, all loaded packages and their version:
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252
## [3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C
## [5] LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] compute.es_0.2-4 tidyselect_1.0.0 pwr_1.3-0 forcats_0.5.0
## [5] stringr_1.4.0 dplyr_0.8.5 purrr_0.3.3 readr_1.3.1
## [9] tidyr_1.0.2 tibble_2.1.3 ggplot2_3.3.0 tidyverse_1.3.0
## [13] pacman_0.5.1
##
## loaded via a namespace (and not attached):
## [1] xfun_0.12 haven_2.2.0 lattice_0.20-38 colorspace_1.4-1
## [5] vctrs_0.2.4 generics_0.0.2 htmltools_0.4.0 yaml_2.2.1
## [9] rlang_0.4.5 pillar_1.4.3 glue_1.3.2 withr_2.1.2
## [13] DBI_1.1.0 dbplyr_1.4.2 modelr_0.1.6 readxl_1.3.1
## [17] lifecycle_0.2.0 munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0
## [21] rvest_0.3.5 codetools_0.2-16 evaluate_0.14 knitr_1.28
## [25] fansi_0.4.1 broom_0.5.5 Rcpp_1.0.3 scales_1.1.0
## [29] backports_1.1.5 jsonlite_1.6.1 fs_1.3.2 hms_0.5.3
## [33] digest_0.6.25 stringi_1.4.3 grid_3.6.3 cli_2.0.2
## [37] tools_3.6.3 magrittr_1.5 crayon_1.3.4 pkgconfig_2.0.3
## [41] xml2_1.2.5 reprex_0.3.0 lubridate_1.7.4 assertthat_0.2.1
## [45] rmarkdown_2.1 httr_1.4.1 rstudioapi_0.11 R6_2.4.1
## [49] nlme_3.1-145 compiler_3.6.3
Estimate how many respondents are necessary to replicate the main effect of Study 1 with a power of 80% and a Type I error of 5%.
We have seen from our exercise that the program G*Power requires the effect size measure Cohen’s f (Cohen, 1988). The same hold for the pwr package, which we will again, use to conduct power analysis in R. Therefore, we first need to convert information provided by the JCP article into a sample estimate of f.
Information the article provides: “A 2 (oral interference: yes, no; between) × 2 (type of oral interference: popcorn, gum; between) × 2 (product: buying a lotion, donating for an organization; within) mixed ANOVA on the likelihood of choosing an advertised option found a significant main effect of oral interference, F(1, 184) = 26.85, p < .0001, \(\eta^{2}_{p}\) = .13. Thus, we collapsed over product and type […].”
Something to note: Remember what you have learned in statistics (Sarstedt & Mooi, 2019, p. 184). Generally, in a one-way ANOVA as applied in the article the denominator degrees of freedom (df=184 in our case) correspond to n-k, where n is the sample size and k is the number of groups that is compared. Furthermore, the nominator degrees of freedom (df=1 in our case) correspond to k-1, where k is the number of groups that is compared. We extract that the authors compared two groups (2-1=1) for a total of n=186 (186-2=184). You should notice: in the description of Study 1 the authors explained that the original sample contains n=188. Obviously, their description of the ANOVA therefore includes a typo (df=184 should be df=186). For this tutorial we assume that n=188 were used.
From this we are able to extract:
This is all that we need to first calculate the observed effect size Cohen’s f (Cohen, 1988). According the Cohen (1988, p. 281) \(\eta^{2}_{p}\) can easily be converted to f by a simple formula: f= \(\sqrt{\frac{\eta^{2}}{1-\eta^{2}}}\) Therefore, we simply need to solve \(\sqrt{\frac{0.13}{1-0.13}}\)
We achieve this by using the next code chunk. Within this code chunk we use the sqrt() function (place the cursor in the function and press ‘F1’ to see help), which calculates the square root of a number. We feed this function with the value of \(\eta^{2}_{p}\)=0.13. We assign the results to an object that we call ‘Cohen_f’.
In the next line we call the object and simultaneously round the results to 4 digits. This is obtained by using the round() function.
Cohen_f <- sqrt(0.13 / (1 - 0.13))
round(Cohen_f, digits = 4)
## [1] 0.3866
We can see that the calculated effect size 0.3866 perfectly mirrors the one we have seen in the exercise slides using G*Power. A common classification for Cohen’s f is: [0.1 | 0.25[ - small effect, [0.25 | 0.4[ - medium effect, and [0.4 | 1] - large effect.
In the question we are asked for \(\beta\)=20%, which corresponds to 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.anova.test() function from the pwr package. This function needs us to provide 5 arguments to fill its parameters. These are (see help by pressing ‘F1’):
The function, furthermore, assumes us to set one of these arguments to NULL. By doing so, we tell the function to use the remaining 4 parameters to search for the value of the fifth. In our case, we are searching for ‘n’, therefore, we set ‘n=NULL’.
We assign the results of our power analysis to a new object named ‘results’. Then we call for its content.
results <- pwr.anova.test(k = 2, n = NULL, f = Cohen_f, sig.level = 0.05, power = 0.8)
results
Balanced one-way analysis of variance power calculation
k = 2
n = 27.25799
f = 0.3865557
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
We can see that the calculated sample size per group 27.257993 only closely mirrors the one we have seen in G*Power 3 (Faul, Erdfelder, Lang, & Buchner, 2007) at n=28. We have to round upwards to end up with 28. For the total sample size, we just multiply by 2, which gives 56.
In a last step we can visualize the relationship between the expected statistical power and different sample sizes.
For this purpose we apply the plot() function with the ‘results’ object and a catchy label for the x-axes as arguments.
plot(results, xlab = "sample size")
From this plot we can alternatively extract the same information for sample size planning.
In the above-discussed example, we were planning a replication of the original findings from Topolinski, Lindner, & Freudenberg (2014). Here, like in Exercise 2, we assume that the observed effect size Cohen’s f = 0.3866 or \(\eta^{2}_{p}\) = 0.13 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.
Again, we can calculate confidence levels for the observed effect size in order to get an idea of a more conservative sample size for our replication. Like in Exercise 2, we will rely on the compute.es package to do so. However, this package does not provide Cohen’s f, nor \(\eta^{2}_{p}\). However, \(\eta^{2}_{p}\) is equivalent to the R² measure used in regression analysis (Sarstedt & Mooi, 2019, p. 184), which is the squared correlation r between observed and predicted observations. This provides us with a workaround. We first will estamate confidence intervals for the effect size r. Then, we will transform the lower limt to \(\eta^{2}_{p}\), before we finally convert to Cohen’s f that is used by the pwr package.
We can apply the fes() function to obtain confidence intervals for the observed effect size r. We assign the results to an object named ‘CI_effect’.
In the call below we feed the function with 5 arguments:
Hint: The numbers for the sample sizes in each group (n.1=92 and n.2=96) can be found in the section ‘Non-parametric analyses’ of the article.
CI_effect <- fes(f = 26.85, n.1 = 92, n.2 = 96, level = 95, dig = 10, verbose = TRUE)
## Mean Differences ES:
##
## d [ 95 %CI] = 0.7559996 [ 0.4580723 , 1.053927 ]
## var(d) = 0.02280627
## p-value(d) = 1.2837e-06
## U3(d) = 77.51753 %
## CLES(d) = 70.35272 %
## Cliff's Delta = 0.4070544
##
## g [ 95 %CI] = 0.7529471 [ 0.4562227 , 1.049671 ]
## var(g) = 0.02262248
## p-value(g) = 1.2837e-06
## U3(g) = 77.42591 %
## CLES(g) = 70.27803 %
##
## Correlation ES:
##
## r [ 95 %CI] = 0.3535122 [ 0.2207161 , 0.4734402 ]
## var(r) = 0.003818257
## p-value(r) = 1.1764e-06
##
## z [ 95 %CI] = 0.369452 [ 0.2244088 , 0.5144951 ]
## var(z) = 0.005405405
## p-value(z) = 1.1764e-06
##
## Odds Ratio ES:
##
## OR [ 95 %CI] = 3.9402 [ 2.295272 , 6.763982 ]
## p-value(OR) = 1.2837e-06
##
## Log OR [ 95 %CI] = 1.371232 [ 0.8308512 , 1.911612 ]
## var(lOR) = 0.07502963
## p-value(Log OR) = 1.2837e-06
##
## Other:
##
## NNT = 3.761045
## Total N = 188
Along the output of the fes()-function we see for the relevant effect size r that we are additionally provided with a test for significance for the effect size itself. This is an alternative to the classical ANOVA. If we reach at a significant result here, this means that the effect size obtained is significantly greater than 0. Put differently, the effect is assumed to exist in the underlying population.
We now recognize that the effect size in the underlying population is likely to be located somewhere in the interval ranging from 0.2207161 to 0.4734402. Put differently, in the worst case we have to plan for a replication of an effect that is, indeed, only of the size of r = 0.2207161.
Now, we are transforming the lower limit of r to \(\eta^{2}_{p}\) and assign the result to an object named ‘Effect_size_lower_limit’. Additionally, we print-out the result.
Effect_size_lower_limit <- CI_effect$l.r^2
Effect_size_lower_limit
## [1] 0.04871562
In R, the ‘^’ symbol is used to calculate the square of a numeric value.
Finally, we transform to Cohen’s f (see formula above).
Effect_size_lower_limit <- sqrt(Effect_size_lower_limit / (1 - Effect_size_lower_limit))
Effect_size_lower_limit
## [1] 0.2262971
We can briefly calculate how this will change our assumptions about the minimum sample size. We, again, use the pwr.anova.test() function (see above). This time, we hand in ‘Effect_size_lower_limit’ as the argument for the effect size parameter f.
results_worst_case <- pwr.anova.test(k = 2, n = NULL, f = Effect_size_lower_limit, sig.level = 0.05, power = 0.8)
results_worst_case
Balanced one-way analysis of variance power calculation
k = 2
n = 77.60566
f = 0.2262971
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
We can see that the more conservative sample size per group is 77.6056598, which makes a replication of the original findings harder. We should ensure a total sample size of at minimum n=156.
Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd ed.). Lawrence Erlbaum.
Faul, F., Erdfelder, E., Lang, A.-G., & Buchner, A. (2007). G* power 3: A flexible statistical power analysis program for the social, behavioral, and biomedical sciences. Behavior Research Methods, 39(2), 175–191. doi: 10.3758/BF03193146
Sarstedt, M., & Mooi, E. (2019). A concise guide to market research: The process, data, and methods using ibm spss statistics (3. ed.). Berlin - Heidelberg: Springer.
Topolinski, S., Lindner, S., & Freudenberg, A. (2014). Popcorn in the cinema: Oral interference sabotages advertising effects. Journal of Consumer Psychology, 24(2), 169–176. doi: 10.1016/j.jcps.2013.09.008