This script assesses and contrasts the stimulus-response characteristics of the SPARS, sensation rating scale (SRS), and pain NRS, through the use of descriptive plots.

The three scales measure the following ranges of somatic sensation:

  • pain NRS: 0 (no pain) to 100 (worst pain you can imagine)

  • SRS: 0 (no sensation) to 100 (pain)

  • SPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)

Because the the stimulus range was centred on the pre-determined pain threshold of each participant (compared to the fixed range of intensities used in Trial A), all analyses use the rank order of the nine stimulus intensities each participant was exposed to rather than the absolute intensities of the stimuli used.

The experimental design involved exposing each participant to four successive experimental blocks of 27 trials (laser stimulations) each for each of the three measurement scales. The sequence of stimulus intensities used within each block was pre-determined, and differed between blocks. The order of in which the measurement scales were assessed was randomized, but for convenience of reporting, the plots are always shown in the order: pain NRS_P, SRS, and SPARS.


Import and inspect data

# Import
data <- read_rds('./data-cleaned/SPARS_B.rds')

# Inspect
glimpse(data)
## Observations: 2,256
## Variables: 10
## $ PID               <chr> "ID01", "ID01", "ID01", "ID01", "ID01", "ID0...
## $ scale             <chr> "SPARS", "SPARS", "SPARS", "SPARS", "SPARS",...
## $ block_number      <int> 2, 2, 2, 4, 4, 4, 6, 6, 6, 8, 8, 8, 11, 11, ...
## $ trial_number      <int> 9, 15, 23, 7, 20, 25, 9, 18, 22, 3, 17, 23, ...
## $ intensity         <dbl> 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2....
## $ intensity_char    <chr> "2.25", "2.25", "2.25", "2.25", "2.25", "2.2...
## $ intensity_rank    <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ rating            <dbl> -31, -20, -48, -48, -21, -23, -48, -45, -47,...
## $ rating_positive   <dbl> 19, 30, 2, 2, 29, 27, 2, 5, 3, 0, 1, 3, 50, ...
## $ rating_equivalent <dbl> -31, -20, -48, -48, -21, -23, -48, -45, -47,...

Clean and transform data

############################################################
#                                                          #
#                          Clean                           #
#                                                          #
############################################################
data %<>%
  # Select required columns
  select(PID, scale, block_number, intensity_rank, rating, rating_equivalent) %>% 
  # Recode scales to CNRS_NP = NRS, NRS_NP = SRS
  mutate(scale = case_when(
      scale == 'NRS_NP' ~ 'SRS',
      scale == 'CNRS_P' ~ 'NRS',
      scale == 'SPARS' ~ 'SPARS'
  ))

############################################################
#                                                          #
#                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                    #
#                                                          #
############################################################
# Add sequential block numbers, tagged by scale
data %<>%
    # Group by ID and scale
    group_by(PID, scale) %>% 
    # Arrange blocks
    arrange(block_number) %>%
    # Add sequential 
    mutate(block_sequential = dense_rank(block_number),
           block_sequential = paste0(scale, ' - ', block_sequential)) %>%
    ungroup()

# Calculate the participant trimean for each scale (using 'rating_equivalent')
data_tm <- data %>% 
  group_by(PID, scale, intensity_rank) %>%
  summarise(tri_mean = tri.mean(rating_equivalent)) %>%
  ungroup()

Exploratory plots

Raw intensity ratings (participant-level)

Data are shown on their original scales: pain NRS: 0 to 100, SRS: 0 to 100, and SPARS: -50 to +50.

# Generate plots
plot_raw <- data %>%
  group_by(PID) %>%
  nest() %>%
  mutate(plots = map2(.x = data,
                      .y = PID,
                      ~ ggplot(data = .x) +
                          aes(x = intensity_rank,
                              y = rating,
                              colour = scale,
                              fill = scale) +
                          geom_smooth(se = FALSE) +
                          geom_point() +
                          scale_y_continuous(limits = c(-50, 100)) +
                          scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
                          scale_fill_manual(name = 'Scale',
                                            values = grey_pal) +
                          scale_colour_manual(name = 'Scale',
                                              values = grey_pal) +
                          labs(title = paste0(.y, ': Raw participant-level stimulus-response ratings on the pain NRS, SRS and SPARS'),
                               subtitle = 'Plots are conditioned on measurement scale and experimental block\npain NRS: 0 (no pain) to 100 (worst pain you can imagine)\nSRS: 0 (no sensation) to 100 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
                               x = 'Stimulus intensity (rank)',
                               y = 'Intensity rating') +
                          geom_hline(yintercept = 0, linetype = 2) +
                          facet_wrap(~ block_sequential)))

# Print plots
walk(plot_raw$plots, ~ print(.x))

Equivalent units ratings (participant-level)

Raw scores for the pain NRS and SRS scales (both rated on 0 to 100 scales) were converted to SPARS equivalent ranges (-50 to +50). The scaling of the pain NRS and SRS were as follows:

  • pain NRS: raw 0 to 100 scores were converted to a 0 to 50 range by dividing 2, such that after the scaling, 0 = no pain and 50 = worst pain you can imagine. This is equivalent to the positive range of the SPARS.

  • SRS: raw 0 to 100 scores were converted to a -50 to 0 range by subtracting 100, and then dividing 2, such that after the scaling, -50 = no sensation and 0 = pain. This is equivalent to the negative range of the SPARS.

# Generate plots
plot_equi <- data %>%
  group_by(PID) %>%
  nest() %>%
  mutate(plots = map2(.x = data,
                      .y = PID,
                      ~ ggplot(data = .x) +
                          aes(x = intensity_rank,
                              y = rating_equivalent,
                              colour = scale,
                              fill = scale) +
                          geom_smooth(se = FALSE) +
                          geom_point() +
                          scale_y_continuous(limits = c(-50, 50)) +
                          scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
                          scale_fill_manual(name = 'Scale',
                                            values = grey_pal) +
                          scale_colour_manual(name = 'Scale',
                                              values = grey_pal) +
                          labs(title = paste0(.y, ': Scaled participant-level stimulus-response ratings on the pain NRS, SRS and SPARS'),
                               subtitle = 'Plots are conditioned on measurement scale and experimental block\npain NRS: 0 (no pain) to 50 (worst pain you can imagine)\nSRS: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
                               x = 'Stimulus intensity (rank)',
                               y = 'Scaled intensity rating (-50 to 50)') +
                          geom_hline(yintercept = 0, linetype = 2) +
                          facet_wrap(~ block_sequential)))

# Print plots
walk(plot_equi$plots, ~ print(.x))

Tukey trimean plots (participant-level)

For each participant, we calculated the Tukey trimean of the intensity rating at each stimulus intensity and for each of the measurement scales. The scaled versions of the pain NRS and SRS were used.

# Generate plot
plot_tm <- data_tm %>%
    # Nest
    group_by(PID) %>%
    nest() %>%
    # Plot
    mutate(plot = map2(.x = data,
                       .y = PID,
                       ~ ggplot(data = ..1) +
                           aes(x = intensity_rank,
                               y = tri_mean) +
                           geom_smooth(method = 'loess',
                                       se = FALSE,
                                       colour = '#656565') +
                           geom_point(shape = 21,
                                      size = 3,
                                      fill = '#CCCCCC') +
                           scale_fill_manual(name = 'Scale',
                                            values = grey_pal) +
                           scale_colour_manual(name = 'Scale',
                                              values = grey_pal) +
                           scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
                           scale_y_continuous(limits = c(-50, 50)) +
                           labs(title = paste0(.y, ': Scaled participant-level stimulus-response rating Tukey trimeans on the pain NRS, SRS and SPARS'),
                                subtitle = 'White points: Tukey trimeans | Grey curve: loess curve\nPlots are conditioned on the three scales\npain NRS: 0 (no pain) to 50 (worst pain you can imagine)\nSRS: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
                                x = 'Stimulus intensity (rank)',
                                y = 'Scaled intensity rating (-50 to 50)') +
                           geom_hline(yintercept = 0, linetype = 2) +
                           facet_wrap(~ scale, ncol = 3)))

# Print plots
walk(plot_tm$plot, ~ print(.x))

Summary plots (group-level)

# Generate plot
data_tm %>%
  ggplot(data = .) +
    aes(x = intensity_rank,
        y = tri_mean,
        colour = scale,
        fill = scale) +
    geom_smooth(method = 'loess',
                se = FALSE) +
    geom_point(size = 2,
               position = position_dodge(width = 0.2)) +
    scale_fill_manual(name = 'Measurement scale',
                      values = grey_pal) +
    scale_colour_manual(name = 'Measurement scale',
                        values = grey_pal) +
    scale_y_continuous(limits = c(-50, 50)) +
    scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
    labs(title = 'Scaled group-level stimulus-response ratings on the pain NRS, SRS and SPARS',
         subtitle = 'Points (dodged for clarity): Tukey trimeans | Curves: loess lines\npain NRS: 0 (no pain) to 50 (worst pain you can imagine)\nSRS: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
         x = 'Stimulus intensity (rank)',
         y = 'Scaled intensity rating (-50 to 50)') +
    geom_hline(yintercept = 0, linetype = 2) 

# Publication plot
p <- data_tm %>%
    ggplot(data = .) +
    aes(x = intensity_rank,
        y = tri_mean,
        colour = scale) +
    geom_segment(x = 0.8, xend = 9, 
                 y = 0, yend = 0, 
                 size = 0.6,
                 linetype = 2,
                 colour = '#000000') +
    geom_point(size = 2,
               position = position_dodge(width = 0.4)) +
    geom_smooth(method = 'loess',
                se = FALSE,
                size = 1) +
    geom_segment(x = 0.7, xend = 0.7, 
                 y = -50.15, yend = 50.15, 
                 size = 1.2,
                 colour = '#000000') +
    geom_segment(x = 0.995, xend = 9.006, 
                 y = -55, yend = -55, 
                 size = 1.2,
                 colour = '#000000') +
    labs(x = 'Relative stimulus intensity',
         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.7, 9.2), 
                       expand = c(0, 0),
                       breaks = seq(from = 1, to = 9, by = 1)) +
    scale_fill_manual(name = 'Scale',
                      values = grey_pal) +
    scale_colour_manual(name = 'Scale',
                        values = grey_pal) +
    theme_bw() +
    theme(panel.border = element_blank(),
          panel.grid = element_blank(),
          legend.text = element_text(size = 12),
          legend.background = element_rect(colour = '#000000'),
          legend.position = c(0.11, 0.9),
          legend.title = element_blank(),
          axis.text = element_text(size = 16,
                                    colour = '#000000'),
          axis.title = element_text(size = 16,
                                    colour = '#000000'))

ggsave(filename = 'figures/figure_11.pdf',
       plot = p,
       width = 6,
       height = 5)
# Generate plot
data %>%
  ggplot(data = .) +
    aes(x = factor(intensity_rank),
        y = rating_equivalent,
        colour = scale,
        fill = scale) +
    geom_boxplot(alpha = 0.6) +
    scale_fill_manual(name = 'Measurement scale',
                      values = grey_pal) +
    scale_colour_manual(name = 'Measurement scale',
                        values = grey_pal) +
    scale_y_continuous(limits = c(-50, 50)) +
    labs(title = 'Scaled group-level stimulus-response ratings on the pain NRS, SRS and SPARS',
         subtitle = 'pain NRS: 0 (no pain) to 50 (worst pain you can imagine)\nSRS: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
         x = 'Stimulus intensity (rank)',
         y = 'Scaled intensity rating (-50 to 50)') +
    geom_hline(yintercept = 0, linetype = 2)


Summary

In general, the ratings provided on the SPARS were less than those provided on the pain NRS (even at higher stimulus intensities), and greater than those provided on the SRS (even at lower stimulus intensities). This indicates that the extra dimensionality of the SPARS compared to the other two scales (the SPARS allowing the ratings of stimulus intensities from no sensation to worst pain imaginable) leads to a compression of the numeric value of the ratings into a narrower band around 0 (pain threshold) compared to the two polar scales (measuring either noxious or non-noxious stimulus ranges).

However, the trend lines indicate that the SPARS scale is reporting the same information as each of the other two scales in their respective domains (sensation of noxious and non-noxious stimuli), but the sensitivity of the SPARS, in each of the two domains, may be lower than they achieved from the two domain-specific scales.

(Please also see: outputs/6B-scale-agreement.html)


Session information

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  forcats_0.3.0   stringr_1.3.1   dplyr_0.7.6    
##  [5] purrr_0.2.5     readr_1.1.1     tidyr_0.8.1     tibble_1.4.2   
##  [9] ggplot2_3.0.0   tidyverse_1.2.1 magrittr_1.5   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.18     cellranger_1.1.0 pillar_1.3.0     compiler_3.5.1  
##  [5] plyr_1.8.4       bindr_0.1.1      tools_3.5.1      digest_0.6.16   
##  [9] lubridate_1.7.4  jsonlite_1.5     evaluate_0.11    nlme_3.1-137    
## [13] gtable_0.2.0     lattice_0.20-35  pkgconfig_2.0.2  rlang_0.2.2     
## [17] cli_1.0.0        rstudioapi_0.7   yaml_2.2.0       haven_1.1.2     
## [21] withr_2.1.2      xml2_1.2.0       httr_1.3.1       knitr_1.20      
## [25] hms_0.4.2        rprojroot_1.3-2  grid_3.5.1       tidyselect_0.2.4
## [29] glue_1.3.0       R6_2.2.2         readxl_1.1.0     rmarkdown_1.10  
## [33] modelr_0.1.2     backports_1.1.2  scales_1.0.0     htmltools_0.3.6 
## [37] rvest_0.3.2      assertthat_0.2.0 colorspace_1.3-2 labeling_0.3    
## [41] stringi_1.2.4    lazyeval_0.2.1   munsell_0.5.0    broom_0.5.0     
## [45] crayon_1.3.4