if (!"pacman" %in% rownames(installed.packages())) install.packages("pacman")
library(pacman)
pacman::p_load(tidyverse, SensoMineR, purrr, sensR)
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.
If this grabs your attention
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!
R analysis script presenting the examples studied in Chapter 6 of the Sensory Marketing lecture
In this exercise students should practice analyzing data with R and RStudio. These skills will serve them well in the progress of their study time and in the future. The internet provides countless further examples. Besides, there exist very good text books (Chapman & Feit, 2019; Lê & Worch, 2018; Næs, Brockhoff, & Tomic, 2010).
1 Loading packages that we will use
Note
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 to your machine. If yes, the corresponding package will be loaded. If not, R will install the package (as long as the package is still maintained on the Comprehensive R Archive Network CRAN.
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., tidyverse (Wickham et al., 2019), sensR (Christensen & Brockhoff, 2023), SensoMineR (Husson, Le, & Cadoret, 2020)). This only works if all packages are still available on CRAN.
Here is the R session info which gives you information on my machine, all loaded packages and their version:
R version 4.3.2 (2023-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22631)
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] sensR_1.5-3 SensoMineR_1.26 FactoMineR_2.9 compute.es_0.2-5
[5] pwr_1.3-0 htmltools_0.5.7 quarto_1.3 ggpubr_0.6.0
[9] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.3
[13] purrr_1.0.2 readr_2.1.4 tidyr_1.3.0 tibble_3.2.1
[17] ggplot2_3.4.4 tidyverse_2.0.0 labelled_2.12.0 knitr_1.45
[21] kableExtra_1.3.4 haven_2.5.3 pacman_0.5.1
loaded via a namespace (and not attached):
[1] tidyselect_1.2.0 viridisLite_0.4.2 fastmap_1.1.1
[4] TH.data_1.1-2 digest_0.6.33 timechange_0.2.0
[7] estimability_1.4.1 lifecycle_1.0.4 cluster_2.1.4
[10] multcompView_0.1-9 survival_3.5-7 processx_3.8.2
[13] magrittr_2.0.3 compiler_4.3.2 rlang_1.1.2
[16] tools_4.3.2 utf8_1.2.4 yaml_2.3.7
[19] ggsignif_0.6.4 htmlwidgets_1.6.2 scatterplot3d_0.3-44
[22] plyr_1.8.9 xml2_1.3.5 KernSmooth_2.23-22
[25] abind_1.4-5 multcomp_1.4-25 numDeriv_2016.8-1.1
[28] withr_2.5.2 grid_4.3.2 fansi_1.0.5
[31] AlgDesign_1.2.1 xtable_1.8-4 colorspace_2.1-0
[34] gtools_3.9.4 emmeans_1.8.9 scales_1.2.1
[37] MASS_7.3-60 flashClust_1.01-2 cli_3.6.1
[40] mvtnorm_1.2-3 rmarkdown_2.25 generics_0.1.3
[43] rstudioapi_0.15.0 reshape2_1.4.4 httr_1.4.7
[46] tzdb_0.4.0 splines_4.3.2 rvest_1.0.3
[49] vctrs_0.6.4 webshot_0.5.5 Matrix_1.6-3
[52] sandwich_3.0-2 jsonlite_1.8.7 carData_3.0-5
[55] car_3.1-2 hms_1.1.3 rstatix_0.7.2
[58] ggrepel_0.9.4 systemfonts_1.0.5 glue_1.6.2
[61] codetools_0.2-19 ps_1.7.5 DT_0.30
[64] stringi_1.8.1 gtable_0.3.4 later_1.3.1
[67] munsell_0.5.0 pillar_1.9.0 R6_2.5.1
[70] evaluate_0.23 lattice_0.22-5 backports_1.4.1
[73] leaps_3.1 broom_1.0.5 Rcpp_1.0.11
[76] svglite_2.1.2 xfun_0.41 zoo_1.8-12
[79] pkgconfig_2.0.3 2 Example for Discrimination test - Triangle test.
Triangle tests are very popular in sensory product research as tests for discrimination, and help to evaluate whether consumers can significantly distinguish between product alternatives (see, e.g., Lawless & Heymann, 2010, Chapter 4.4). In our lecture example on Lasagna variants, we observe the following results on the basis of n=100 consumers:
| Correctly identified odd sample | Incorrect response | |
|---|---|---|
| Observed | 43 (O1) | 57 (O2) |
| Expected | 33.\(\bar{3}\) (E1) | 66.\(\bar{6}\) (E2) |
In our example, analysis draws on a (Pearson) \(\chi^{2}\)-test (without continuity correction). Here, the test statistic is calculated as follows:
\(\chi^{2}_{(1)}=\sum\frac{(observed-expected)^2}{expected}=[\frac{(O_{1}-E_{1})^2}{E_{1}}]+[\frac{(O_{2}-E_{2})^2}{E_{2}}]=[\frac{(43-33.\bar{3})^2}{33.\bar{3}}]+\frac{(57-66.\bar{6})^2}{66.\bar{6}}]=2.803+1.402=4.205\)
We know from lectures in statistics that the df in a \(\chi^{2}\)-test equal: (number of rows -1) * (number of columns -1). Thus, df=(2-1)*(2-1)=1. The degrees of freedom are 1.
The package sensR is providing functions to handle nearly all types of discrimination tests, including Triangle tests. We will use the function discrim(). This function allows us to specify many arguments (see help by pressing ‘F1’). We will do so for
- correct = the number of correct answers
- total = total number of responses
- conf.level = significance level (Type I error probability \(\alpha\))
- method = the discrimination protocol applied (see lecture slides for the various options)
- statistic = The type of analysis strategy to apply for the data (whether one wants to use normal distribution or Chi² etc.)
- test = Whether a test for similarity or difference should be obtained.
Below, we assign the results of our Triangle Discrimination test to a new object named ‘triangle_test’. Then we call for its content.
triangle_test <- discrim(
correct= 43,
total= 100,
conf.level = 0.95,
method = "triangle",
statistic = "score",
test = "difference")
triangle_test
Estimates for the triangle discrimination protocol with 43 correct
answers in 100 trials. One-sided p-value and 95 % two-sided confidence
intervals are based on the Pearson and score statistics.
Estimate Std. Error Lower Upper
pc 0.430 0.04951 0.3333 0.5328
pd 0.145 0.07426 0.0000 0.2992
d-prime 1.075 0.30276 0.0000 1.6355
Result of difference test:
Pearson X-square statistic = 4.205, df = 1, p-value: 0.02015
Alternative hypothesis: d-prime is greater than 0
Interpretation
As the results of the \(\chi^{2}\)-test is significant (0.02), we have to reject H0. Thus, consumers can reliably distinguish between the new and the old lasagna.
- In the output above, pc means the percentage of correct discriminators (43/100).
- pd stands for the estimated proportion of individuals in the (relevant) population that would detect the product difference. Which is in the case of a Triangle test: \(pd=\frac{pc-1/3}{2/3}\). Note that the confidence interval for pd has an upper limit of 29.92%. Thus, in the worst case, the true proportion of discriminators is almost as high as 30%. Schlich (1993) proposed the following limits for pd as small, medium, large differences: 0.25, 0.375, and 0.50.
- The Thurstonian approach of transforming the number of correct answers into an estimate, called d-prime of the underlying (relative) sensory difference is an attempt to overcome the concepts of pc and pd, since these are dependent on the concrete test protocol one has applied (e.g., Triangle vs. 3-AFC tests). The higher the value the higher is the degree of difference between the two products that were tested.
3 Example of a Post-hoc Power Analyis for a Triangle Test
To get an idea about the statistical Power of a Similarity or Dissimilarity test protocol, the sensR package also provides functions that facilitate sample size planning for discrimination tests (i.e., statistical power analysis). To provide an example, we will use the d-prime estimate from the Triangle test on lasagna above (1.07) as input the determine the statistical power of the Triangle test with n=100.
In the code chunk below, we use the obtained d-prime as input for the function d.primePwr(). Furthermore, we specify all other aspects of our exemplary Triangle test.
power_triangle <- d.primePwr(
d.primeA = 1.075,
sample.size = 100,
alpha = 0.05,
method = "triangle",
test = "difference", statistic = "exact"
)
power_triangle[1] 0.617578
Interpretation
With the observed effect size for the difference, our Triangle test protocol only has a power of 61.76%. Thus, if the observed effect is the true effect in the underlying population, then we would have this chance to find a significant difference in a Triangle test with n=100 and \(\alpha\)=5%, which is insufficient.
As a final remark, we can estimate how many consumers should be tested in order to obtain a sufficient statistical power based on the observed effect size for the difference, a Triangle test, and \(\alpha\)=5%. For this purpose, we first establish a list of varying sample sizes below. Afterwards, we use this list as input for the sample.size argument of the d.primePwr() function. Usually, this function does not accept a list of values as arguments. To handle this problem, we additionally apply the map() function provided by the purrr package (Wickham & Henry, 2023). In the corresponding call of d.primePwr() we set sample.size = ., to tell the map function to use each element of the input list as input for the sample size argument in repeated calls of d.primePwr().
sample_sizes <- list(n100 = 100,
n200 = 200,
n300 = 300,
n400 = 400,
n500 = 500)
sample_sizes %>%
map(~ d.primePwr(
d.primeA = 1.075,
sample.size = .,
alpha = 0.05,
method = "triangle",
test = "difference", statistic = "exact"
))$n100
[1] 0.617578
$n200
[1] 0.8583079
$n300
[1] 0.9552637
$n400
[1] 0.9888657
$n500
[1] 0.9972308
Interpretation
We can see that with a sample size of n=200 and above, we realize a power of over 85%.
References
Chapman, C., & Feit, E. M. (2019). R for marketing research and analytics. doi: 10.1007/978-3-030-14316-9
Christensen, R. H. B., & Brockhoff, P. B. (2023). sensR: Thurstonian models for sensory discrimination. Retrieved from https://CRAN.R-project.org/package=sensR
Husson, F., Le, S., & Cadoret, M. (2020). SensoMineR: Sensory data analysis. Retrieved from https://CRAN.R-project.org/package=SensoMineR
Lawless, H. T., & Heymann, H. (2010). Sensory evaluation of food: Principles and practices (2nd ed.). New York, Dordrecht, Heidelberg, London: Springer.
Lê, S., & Worch, T. (2018). Analyzing sensory data with r. doi: 10.1201/9781315373416
Næs, T., Brockhoff, P. B., & Tomic, O. (2010). Statistics for sensory and consumer science. doi: 10.1002/9780470669181
Rinker, T. W., & Kurkiewicz, D. (2018). Pacman: Package management for r. Retrieved from http://github.com/trinker/pacman
Schlich, P. (1993). Risk tables for discrimination tests. Food Quality and Preference, 4(3), 141–151. doi: 10.1016/0950-3293(93)90157-2
Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., … Yutani, H. (2019). Welcome to the tidyverse. 4, 1686. doi: 10.21105/joss.01686
Wickham, H., & Henry, L. (2023). Purrr: Functional programming tools. Retrieved from https://CRAN.R-project.org/package=purrr
Citation
For attribution, please cite this work as:
Prof. Dr. Lichters, M. (2023, November 16). Sensory Marketing &
Product Innovation: Chapter 6 - Sensory Product Research - Examples.
Retrieved from https://rpubs.com/M_Lichters/smpiexercise3