In this script, we combined Figure 4 from (suppl_05_4A-stimulus-response-1.Rmd) and Figure 5 from (suppl_06_4A-stimulus-response-2.Rmd) into a single figure, Figure 4, which shows the individual participants’ Tukey trimeans by stimulus intensity (black dots), the group means by intensity (white circles), the (best fit) cubic model (black line) and its 95% CI (grey band).
All figures downstream of Figure 5 have been re-labelled accordingly.
# Import
data <- read_rds('./data-cleaned/SPARS_A.rds')
# Inspect
glimpse(data)
## Observations: 1,927
## Variables: 19
## $ PID <chr> "ID01", "ID01", "ID01", "ID01", "ID01", "ID0...
## $ block <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A",...
## $ block_order <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,...
## $ trial_number <dbl> 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, ...
## $ intensity <dbl> 3.00, 2.25, 4.00, 3.25, 2.75, 2.25, 2.75, 4....
## $ intensity_char <chr> "3.00", "2.25", "4.00", "3.25", "2.75", "2.2...
## $ rating <dbl> -40, -25, 10, 2, -10, -25, -20, 10, -25, -50...
## $ rating_positive <dbl> 10, 25, 60, 52, 40, 25, 30, 60, 25, 0, 25, 3...
## $ EDA <dbl> 75270.55, 43838.67, 35967.67, 26720.61, 1931...
## $ age <dbl> 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, ...
## $ sex <dbl> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
## $ panas_positive <dbl> 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, ...
## $ panas_negative <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...
## $ dass42_depression <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ dass42_anxiety <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ dass42_stress <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ pcs_magnification <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,...
## $ pcs_rumination <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, ...
## $ pcs_helplessness <dbl> 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, ...
We performed a basic clean-up of the data, and then calculated Tukey trimean at each stimulus intensity for each participant (participant average), and finally the median of the trimeans at each stimulus intensity across participants (group average).
############################################################
# #
# Clean #
# #
############################################################
data %<>%
# Select required columns
select(PID, block, block_order, trial_number, intensity, intensity_char, rating)
############################################################
# #
# Calculate 'Tukey trimean' #
# #
############################################################
# Define tri.mean function
tri.mean <- function(x) {
# Calculate quantiles
q1 <- quantile(x, probs = 0.25, na.rm = TRUE)[[1]]
q2 <- median(x, na.rm = TRUE)
q3 <- quantile(x, probs = 0.75, na.rm = TRUE)[[1]]
# Calculate trimean
tm <- (q2 + ((q1 + q3) / 2)) / 2
# Convert to integer
tm <- as.integer(round(tm))
return(tm)
}
############################################################
# #
# Generate core data #
# #
############################################################
# Calculate the participant average
data_tm <- data %>%
group_by(PID, intensity) %>%
summarise(tri_mean = tri.mean(rating)) %>%
ungroup()
# Calculate the group average
data_group <- data_tm %>%
group_by(intensity) %>%
summarise(median = median(tri_mean)) %>%
ungroup()
Generate the 3rd-order (cubic) polynomial regression model (random slope and intercept), and get predicted values for plotting.
# Generate model
lmm3b <- lmer(tri_mean ~ poly(intensity, 3) + (intensity | PID),
data = data_tm)
# Generate predicted values
predicted <- ggeffect(model = lmm3b,
terms = 'intensity',
ci.lvl = 0.95)
p <- ggplot() +
geom_ribbon(data = predicted,
aes(x = x,
ymin = conf.low,
ymax = conf.high),
fill = '#CCCCCC') +
geom_point(data = data_tm,
aes(x = intensity,
y = tri_mean),
position = position_jitter(width = 0.05)) +
geom_line(data = predicted,
aes(x = x,
y = predicted),
size = 0.8) +
geom_point(data = data_group,
aes(x = intensity,
y = median),
shape = 21,
size = 5,
stroke = 1,
fill = '#FFFFFF') +
geom_segment(aes(x = 0.8, xend = 0.8,
y = -50.25, yend = 50.25),
size = 1.2) +
geom_segment(aes(x = 0.995, xend = 4.006,
y = -55, yend = -55),
size = 1.2) +
geom_segment(aes(x = 0.8, xend = 4,
y = 0, yend = 0),
size = 0.6,
linetype = 2) +
labs(x = 'Stimulus intensity (J)',
y = 'SPARS rating (-50 to 50)') +
scale_y_continuous(limits = c(-55, 50.25),
expand = c(0, 0),
breaks = c(-50, -25, 0, 25, 50)) +
scale_x_continuous(limits = c(0.8, 4.2),
expand = c(0, 0),
breaks = seq(from = 1, to = 4, by = 0.5)) +
theme_bw() +
theme(panel.border = element_blank(),
panel.grid = element_blank(),
axis.text = element_text(size = 16,
colour = '#000000'),
axis.title = element_text(size = 16,
colour = '#000000'))
# Print plot
print(p)
# Save plot
ggsave(filename = 'figures/figure_4.pdf',
plot = p,
width = 6,
height = 5)
```r
sessionInfo()
```
```
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
##
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-r0.2.19.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=C
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] bindrcpp_0.2.2 ggeffects_0.5.0 lmerTest_3.0-1 lme4_1.1-18-1
## [5] Matrix_1.2-14 forcats_0.3.0 stringr_1.3.1 dplyr_0.7.6
## [9] purrr_0.2.5 readr_1.1.1 tidyr_0.8.1 tibble_1.4.2
## [13] ggplot2_3.0.0 tidyverse_1.2.1 magrittr_1.5
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-137 lubridate_1.7.4 httr_1.3.1
## [4] rprojroot_1.3-2 numDeriv_2016.8-1 tools_3.5.1
## [7] TMB_1.7.14 backports_1.1.2 R6_2.2.2
## [10] sjlabelled_1.0.13 lazyeval_0.2.1 colorspace_1.3-2
## [13] nnet_7.3-12 withr_2.1.2 tidyselect_0.2.4
## [16] emmeans_1.2.3 compiler_3.5.1 cli_1.0.0
## [19] rvest_0.3.2 xml2_1.2.0 sandwich_2.5-0
## [22] effects_4.0-3 scales_1.0.0 mvtnorm_1.0-8
## [25] ggridges_0.5.0 digest_0.6.16 minqa_1.2.4
## [28] rmarkdown_1.10 stringdist_0.9.5.1 pkgconfig_2.0.2
## [31] htmltools_0.3.6 pwr_1.2-2 rlang_0.2.2
## [34] readxl_1.1.0 rstudioapi_0.7 bindr_0.1.1
## [37] zoo_1.8-3 jsonlite_1.5 modeltools_0.2-22
## [40] bayesplot_1.6.0 Rcpp_0.12.18 munsell_0.5.0
## [43] prediction_0.3.6 stringi_1.2.4 multcomp_1.4-8
## [46] yaml_2.2.0 snakecase_0.9.2 carData_3.0-1
## [49] MASS_7.3-50 plyr_1.8.4 grid_3.5.1
## [52] parallel_3.5.1 sjmisc_2.7.4 crayon_1.3.4
## [55] lattice_0.20-35 haven_1.1.2 splines_3.5.1
## [58] sjstats_0.17.0 hms_0.4.2 knitr_1.20
## [61] pillar_1.3.0 estimability_1.3 codetools_0.2-15
## [64] stats4_3.5.1 glue_1.3.0 evaluate_0.11
## [67] data.table_1.11.4 modelr_0.1.2 nloptr_1.0.4
## [70] cellranger_1.1.0 gtable_0.2.0 assertthat_0.2.0
## [73] coin_1.2-2 xtable_1.8-3 broom_0.5.0
## [76] survey_3.33-2 survival_2.42-3 glmmTMB_0.2.2.0
## [79] TH.data_1.0-9
```