This script is part 5 of our analysis of the stimulus-response characteristics of the SPARS. In this analysis we examined whether there is a linear relationship between change in stimulus intensity and change in rating magnitude between two successive stimuli, irrespective of the direction of the change?
Descriptive plots of the data are provided in “outputs/suppl_05_4A-stimulus-response-1.html”, modelling of the stimulus-response relationship is described in “outputs/suppl_06_4A-stimulus-response-2.html”, the diagnostics on the final linear mixed model are described in “outputs/suppl_07_4A-stimulus-response-3.html”, the stability of the model is described in “outputs/suppl_08_4A-stimulus-response-4.html”, and the variance in ratings at each stimulus intensity is described in “outputs/suppl_10_4A-stimulus-reponse-6.html”.
############################################################
# #
# Import #
# #
############################################################
data <- read_rds('./data-cleaned/SPARS_A.rds')
############################################################
# #
# Clean #
# #
############################################################
data %<>%
# Select required columns
select(PID, block, block_order, trial_number, intensity, intensity_char, rating)
############################################################
# #
# Define Wald CI function for robust LMM #
# #
############################################################
# Adapted from code provided Ben Bolker on StackExchange: https://stats.stackexchange.com/questions/233800/how-can-i-get-confidence-intervals-for-fixed-effects-using-the-rlmer-function-r
confint.rlmerMod <- function(object, level = 0.95) {
# Extract beta coefficients
beta <- fixef(object)
# Extract names of coefficients
parm <- names(beta)
# Extract standard errors for the coefficients
se <- sqrt(diag(vcov(object)))
# Set level of confidence interval
z <- qnorm((1 + level) / 2)
# Calculate CI
ctab <- cbind(beta - (z * se),
beta + (z * se))
# label column names
colnames(ctab) <- c(paste(100 * ((1 - level) / 2), '%'),
paste(100 * ((1 + level) / 2), '%'))
# Output
return(as.data.frame(ctab[parm, ]))
}
############################################################
# #
# Generate core data #
# #
############################################################
# Calculate lagged data
data_scale <- data %>%
# Get lag-1 stimulus intensity and FEST rating by PID
group_by(PID) %>%
mutate(intensity_lag = lag(intensity),
rating_lag = lag(rating)) %>%
# Ungroup and remove incomplete cases greated by lag function
ungroup() %>%
filter(complete.cases(.))
# Calculate change in stimulus intensity and SPARS rating between successive stimuli
data_scale %<>%
mutate(intensity_delta = intensity - intensity_lag,
rating_delta = rating - rating_lag) %>%
# Determine the direction of the change in stimulus intensity
mutate(change_direction = case_when(
intensity_delta > 0 ~ 'up',
intensity_delta < 0 ~ 'down',
intensity_delta == 0 ~ 'no change'
))
# Filter out 'no change' change_direction
data_reduced <- data_scale %>%
filter(change_direction != 'no change')
# Generate the plots for each individual
data_scale %>%
# Filter out 'no change'
filter(change_direction != 'no change') %>%
ggplot(data = .) +
aes(y = rating_delta,
x = intensity_delta,
colour = change_direction) +
geom_point() +
geom_smooth(method = 'lm') +
scale_color_manual(name = 'Direction of intensity change: ',
labels = c('Down', 'Up'),
values = c("#999999", "#323232")) +
scale_x_continuous(limits = c(-3.5, 3.5),
breaks = seq(from = -3, to = 3, by = 1),
expand = c(0,0)) +
scale_y_continuous(limits = c(-70, 70),
breaks = seq(from = -60, to = 60, by = 20),
expand = c(0,0)) +
labs(title = 'Group-level plot: Change in SPARS rating vs Change in stimulus intensity \nbetween two successive stimuli',
subtitle = 'Grey points: Individual responses for all trials | Coloured line: Linear trend for context\nDark grey points/line: When a stimulus was of less intensity than the preceding stimulus | \nLight grey points/line: When a stimulus was of greater intensity than the preceding stimulus.\nData have been omitted when successive stimuli were of the same intensity (zero change).',
x = expression(Delta~stimulus~intensity~(J)),
y = expression(Delta~SPARS~rating~(predicted))) +
theme(legend.position = 'top')
## Publication plot
p <- data_scale %>%
# Filter out 'no change'
filter(change_direction != 'no change') %>%
ggplot(data = .) +
aes(y = rating_delta,
x = intensity_delta,
colour = change_direction) +
geom_point() +
geom_smooth(method = 'lm',
se = FALSE) +
geom_segment(x = -3.3, xend = 3.05,
y = 0, yend = 0,
size = 0.6,
linetype = 2,
colour = '#000000') +
geom_segment(x = -3.3, xend = -3.3,
y = -90.23, yend = 90.23,
size = 1.2,
colour = '#000000') +
geom_segment(x = -3.005, xend = 3.005,
y = -98, yend = -98,
size = 1.2,
colour = '#000000') +
labs(x = 'Change in stimulus intensity (J)',
y = 'Change in SPARS rating') +
scale_colour_manual(name = 'Direction of change: ',
values = c("#999999", "#323232"),
labels = c('Decreasing intensity', 'Increasing intensity')) +
scale_y_continuous(limits = c(-98, 90.25),
expand = c(0, 0),
breaks = c(-90, -60, -30, 0, 30, 60, 90)) +
scale_x_continuous(limits = c(-3.3, 3.2),
expand = c(0, 0),
breaks = seq(from = -3, to = 3, by = 1)) +
theme_bw() +
theme(legend.title = element_text(size = 14,
face = 'bold'),
legend.text = element_text(size = 12),
legend.background = element_rect(colour = '#000000'),
legend.position = c(0.23, 0.895),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.text = element_text(size = 16,
colour = '#000000'),
axis.title = element_text(size = 16,
colour = '#000000'))
ggsave(filename = 'figures/figure_7.pdf',
plot = p,
width = 6,
height = 5)
Participant-level plot: Change in SPARS rating vs change in stimulus intensity between two successive stimuli
Coloured points: Per trial responses | Coloured line: Linear trend for the group
Green points/line: When a stimulus was less than the preceding stimulus intensity | Orange points/line: When a stimulus was greater than the preceding stimulus intensity.
Data when successive stimuli were of the same intensity (zero change) have been omitted.
# Generate the plots for each individual
scale_plots <- data_scale %>%
# Filter out 'no change'
filter(change_direction != 'no change') %>%
# Nest data by PID
group_by(PID) %>%
nest() %>%
# Plot data
mutate(plot = map2(.x = data,
.y = unique(PID),
~ ggplot(data = .x) +
aes(y = rating_delta,
x = intensity_delta,
colour = change_direction) +
geom_point() +
geom_smooth(method = 'lm') +
scale_color_manual(name = 'Direction of\nintensity change',
labels = c('Down', 'Up'),
values = c('#999999', '#323232')) +
scale_x_continuous(limits = c(-3.5, 3.5),
breaks = seq(from = -3, to = 3, by = 1),
expand = c(0,0)) +
scale_y_continuous(limits = c(-70, 70),
breaks = seq(from = -60, to = 60, by = 20),
expand = c(0,0)) +
labs(subtitle = paste0('[', .y, ']'),
x = expression(Delta~stimulus~intensity~(J)),
y = expression(Delta~SPARS~rating)) +
theme(legend.position = 'none')))
# Print plots
wrap_plots(scale_plots$plot, ncol = 4)
To examine whether the relationship between change in rating intensity and change in stimulus intensity was affected by the direction of the change in stimulus intensity (i.e., whether the preceding stimulus was greater than or less than the current stimulus), we fit a linear mixed model regression (with random intercept) that included an interaction term between change in stimulus intensity and the direction of the change in stimulus intensity (‘up’ or ‘down’):
\[\Delta~SPARS~rating \thicksim \Delta~stimulus~intenisty~\ast~direction~of~change~+~(1~|~participant)\]
Preliminary analysis indicated that influence points may be an issue, and so we modelled the data using robust methods.
We generated two models. The first model included an interaction term between change in stimulus intensity and direction of change, and the second model did not. Comparing the models allows one you check whether the slope of the SPARS rating vs change in stimulus intensity relationship was dependent on the direction of the change in stimulus intensity.
# Linear mixed model with interaction
scale_lmm <- rlmer(rating_delta ~ intensity_delta * change_direction + (1 | PID),
method = 'DASvar',
data = data_reduced)
# Linear mixed model without interaction
scale_lmm2 <- rlmer(rating_delta ~ intensity_delta + change_direction + (1 | PID),
method = 'DASvar',
data = data_reduced)
# Inspect models
## With interaction
summary(scale_lmm)
## Robust linear mixed model fit by DASvar
## Formula: rating_delta ~ intensity_delta * change_direction + (1 | PID)
## Data: data_reduced
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2631 -0.6422 0.0132 0.6522 3.8558
##
## Random effects:
## Groups Name Variance Std.Dev.
## PID (Intercept) 0.0 0.00
## Residual 232.8 15.26
## Number of obs: 1789, groups: PID, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.8809 0.9712 1.937
## intensity_delta 14.1475 0.7171 19.729
## change_directionup -1.6309 1.4269 -1.143
## intensity_delta:change_directionup -1.4151 1.0304 -1.373
##
## Correlation of Fixed Effects:
## (Intr) intns_ chng_d
## intnsty_dlt 0.848
## chng_drctnp -0.681 -0.577
## intnsty_d:_ -0.590 -0.696 -0.051
##
## Robustness weights for the residuals:
## 1416 weights are ~= 1. The remaining 373 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.315 0.585 0.739 0.739 0.901 0.999
##
## Robustness weights for the random effects:
## All 19 weights are ~= 1.
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal II (k = 1.345, s = 10)
## Random Effects, variance component 1 (PID):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal II (k = 1.345, s = 10)
ci_interaction <- confint.rlmerMod(object = scale_lmm)
knitr::kable(ci_interaction,
caption = 'Interaction model: Wald 95% confidence intervals of fixed effects')
2.5 % | 97.5 % | |
---|---|---|
(Intercept) | -0.0226197 | 3.7843721 |
intensity_delta | 12.7420347 | 15.5529405 |
change_directionup | -4.4276617 | 1.1658267 |
intensity_delta:change_directionup | -3.4345943 | 0.6043071 |
## Without interaction
summary(scale_lmm2)
## Robust linear mixed model fit by DASvar
## Formula: rating_delta ~ intensity_delta + change_direction + (1 | PID)
## Data: data_reduced
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3212 -0.6286 0.0315 0.6683 3.8525
##
## Random effects:
## Groups Name Variance Std.Dev.
## PID (Intercept) 0.0 0.00
## Residual 232.7 15.25
## Number of obs: 1789, groups: PID, 19
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.1289 0.7839 1.44
## intensity_delta 13.4273 0.5148 26.08
## change_directionup -1.6806 1.4247 -1.18
##
## Correlation of Fixed Effects:
## (Intr) intns_
## intnsty_dlt 0.754
## chng_drctnp -0.882 -0.854
##
## Robustness weights for the residuals:
## 1407 weights are ~= 1. The remaining 382 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.311 0.580 0.757 0.744 0.912 0.999
##
## Robustness weights for the random effects:
## All 19 weights are ~= 1.
##
## Rho functions used for fitting:
## Residuals:
## eff: smoothed Huber (k = 1.345, s = 10)
## sig: smoothed Huber, Proposal II (k = 1.345, s = 10)
## Random Effects, variance component 1 (PID):
## eff: smoothed Huber (k = 1.345, s = 10)
## vcp: smoothed Huber, Proposal II (k = 1.345, s = 10)
ci_main <- confint.rlmerMod(object = scale_lmm2)
knitr::kable(ci_main,
caption = 'Main effects model: Wald 95% confidence intervals of fixed effects')
2.5 % | 97.5 % | |
---|---|---|
(Intercept) | -0.4075402 | 2.665404 |
intensity_delta | 12.4183089 | 14.436217 |
change_directionup | -4.4728418 | 1.111722 |
## Compare the models
### Prepare for plotting
ci_interaction %<>%
rownames_to_column(var = 'Fixed effects') %>%
mutate(model = 'interaction') %>%
filter(`Fixed effects` != 'intensity_delta:change_directionup')
ci_main %<>%
rownames_to_column(var = 'Fixed effects') %>%
mutate(model = 'main effects')
ci <- bind_rows(ci_main, ci_interaction)
ci %>%
mutate(`Fixed effects` = fct_relevel(`Fixed effects`,
'(Intercept)',
'intensity_delta',
'change_directionup')) %>%
ggplot(data = .) +
aes(x = `Fixed effects`,
y = (`2.5 %` + `97.5 %`) / 2,
ymin = `2.5 %`,
ymax = `97.5 %`,
colour = model,
fill = model) +
geom_crossbar(position = position_dodge(width = 1),
alpha = 0.7) +
scale_fill_manual(name = 'Model',
labels = c('Interaction', 'Main effects only'),
values = c('#999999', '#323232')) +
scale_colour_manual(name = 'Model',
labels = c('Interaction', 'Main effects only'),
values = c('#999999', '#323232')) +
labs(title = 'Comparison of the Wald 95% confidence intervals of the fixed effects (intercept and main effects)',
subtitle = 'The middle line in each bar shows the mid-point between the interval edges',
y = 'Units')
Using the confidence intervals of the fixed effect estimates to assess significance (intervals that excluded zero were considered significant at the 5% level), we assessed the interaction term to be not significant. When comparing the models, the fixed effect estimates for the two models were similar (see the plot above), and therefore we conclude that the direction of the change in stimulus intensity between two successive stimuli does not influence the slope of the relationship between change in stimulus intensity and change in SPARS rating.
A basic diagnostic assessment did not identify any violations of the model assumptions.
# Generate plots
p_list <- list(plot(scale_lmm2, which = 1),
plot(scale_lmm2, which = 2),
plot(scale_lmm2, which = 3))
x <- plot(scale_lmm2, which = 1) +
scale_colour_viridis_c()
x
## NULL
str(x)
## NULL
class(x)
## [1] "NULL"
# Print plots
walk(p_list, ~print(.x))
# Predicted values from the main effects model
## Generate 'newdata'
new_data <- data_reduced %>%
select(intensity_delta, change_direction, PID)
## Get predictions
yhat <- data_frame(y = predict(object = scale_lmm2,
newdata = new_data))
## Bind to 'new_data/extra_data'
predicted <- bind_cols(yhat, new_data) %>%
rename(x = intensity_delta,
group = change_direction)
## Plot
ggplot(data = predicted) +
aes(x = x,
y = y,
colour = group) +
geom_line() +
geom_point() +
scale_color_manual(name = 'Direction of\nintensity change',
labels = c('Down', 'Up'),
values = c('#999999', '#323232')) +
scale_x_continuous(limits = c(-3, 3),
breaks = seq(from = -3, to = 3, by = 1)) +
scale_y_continuous(limits = c(-60, 60),
breaks = seq(from = -60, to = 60, by = 20)) +
labs(title = 'Plot of predicted change in SPARS ratings vs change stimulus intensity',
x = expression(Delta~stimulus~intensity~(J)),
y = expression(Delta~SPARS~rating~(predicted)))
There is linear scaling between change in stimulus intensity and change in SPARS rating, and the sensitivity of this response (the slope of the curve) is consistent across increasing and decreasing changes in stimulus intensities. The estimate for change in the change in stimulus intensity (intensity_delta) main effect was 13.3 (95% CI: 12.3 to 14.4), indicating that a single unit change (1 J) in stimulus intensity results in 12.3 to 14.4 unit change in SPARS rating.
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 influence.ME_0.9-9 sjPlot_2.6.0
## [4] robustlmm_2.2-1 lmerTest_3.0-1 lme4_1.1-18-1
## [7] Matrix_1.2-14 patchwork_0.0.1 forcats_0.3.0
## [10] stringr_1.3.1 dplyr_0.7.6 purrr_0.2.5
## [13] readr_1.1.1 tidyr_0.8.1 tibble_1.4.2
## [16] 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 TMB_1.7.14
## [7] tools_3.5.1 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 mnormt_1.5-5
## [16] tidyselect_0.2.4 emmeans_1.2.3 compiler_3.5.1
## [19] cli_1.0.0 rvest_0.3.2 xml2_1.2.0
## [22] sandwich_2.5-0 labeling_0.3 effects_4.0-3
## [25] scales_1.0.0 psych_1.8.4 DEoptimR_1.0-8
## [28] mvtnorm_1.0-8 robustbase_0.93-2 ggridges_0.5.0
## [31] digest_0.6.16 foreign_0.8-70 minqa_1.2.4
## [34] rmarkdown_1.10 stringdist_0.9.5.1 pkgconfig_2.0.2
## [37] htmltools_0.3.6 highr_0.7 pwr_1.2-2
## [40] rlang_0.2.2 readxl_1.1.0 rstudioapi_0.7
## [43] bindr_0.1.1 zoo_1.8-3 jsonlite_1.5
## [46] modeltools_0.2-22 bayesplot_1.6.0 Rcpp_0.12.18
## [49] munsell_0.5.0 prediction_0.3.6 fastGHQuad_0.2
## [52] stringi_1.2.4 multcomp_1.4-8 yaml_2.2.0
## [55] snakecase_0.9.2 carData_3.0-1 MASS_7.3-50
## [58] plyr_1.8.4 grid_3.5.1 parallel_3.5.1
## [61] sjmisc_2.7.4 crayon_1.3.4 lattice_0.20-35
## [64] ggeffects_0.5.0 haven_1.1.2 splines_3.5.1
## [67] sjstats_0.17.0 hms_0.4.2 knitr_1.20
## [70] pillar_1.3.0 estimability_1.3 codetools_0.2-15
## [73] stats4_3.5.1 glue_1.3.0 evaluate_0.11
## [76] data.table_1.11.4 modelr_0.1.2 nloptr_1.0.4
## [79] cellranger_1.1.0 gtable_0.2.0 assertthat_0.2.0
## [82] coin_1.2-2 xtable_1.8-3 broom_0.5.0
## [85] survey_3.33-2 viridisLite_0.3.0 survival_2.42-3
## [88] glmmTMB_0.2.2.0 TH.data_1.0-9