Validation of the German PC Scale

Overview

This notebook analyses the German PC scale validation dataset.

Variable	Description
SK01	Hand lowering (0 = Normal severity, 5 = Very severe)
SK02	Moving hands together (0 = No force, 5 = Strong force)
SK03	Mosquito hallucination (0 = No mosquito, 5 = Like a real mosquito)
SK04	Sweet taste (0 = No taste, 5 = Strong taste)
SK05	Sour taste (0 = No taste, 5 = Strong taste)
SK04SK05	Taste mean (calculated in the script)
SK06	Arm rigidity (0 = Normal / no stiffness, 5 = Very stiff)
SK07	Arm immobilization (0 = Normal / no heaviness, 5 = Very heavy)
SK08	Music hallucination (0 = No music, 5 = Clear music)
SK09	Negative visual hallucination (0 = Three balls seen, 5 = Two balls seen)
SK10	Amnesia (0 = Normal memory, 5 = No memory)
SK12	Urge to press space bar (0 = No urge, 5 = Clear urge)
SK13	Space bar amnesia (0 = Normal memory of instructions, 5 = No memory of instructions)
SK12SK13	Space bar mean (calculated in the script)
KK02_01	Items remembered before being asked again
KK03_01	Items remembered after being asked again
SO01	Gender (1 = female, 2 = male, 3 = diverse)
SO02_01	Age
SC02	Answered all questions seriously (0 = yes, 2 = no)
TIME005	Time needed to complete the experiment (seconds)

Import libraries

knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  fig.width = 8,
  fig.height = 5,
  fig.align = "center"
)

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(Hmisc)

## 
## Attache Paket: 'Hmisc'
## 
## Die folgenden Objekte sind maskiert von 'package:dplyr':
## 
##     src, summarize
## 
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     format.pval, units

library(psych)

## 
## Attache Paket: 'psych'
## 
## Das folgende Objekt ist maskiert 'package:Hmisc':
## 
##     describe
## 
## Die folgenden Objekte sind maskiert von 'package:ggplot2':
## 
##     %+%, alpha

library(ggdist)
library(patchwork)
library(ggthemr)

ggthemr("flat")

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the ggthemr package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## ℹ The deprecated feature was likely used in the ggthemr package.
##   Please report the issue to the authors.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Read data

dat <- read.csv("GitHub/german_pc_scale_validation_data.CSV",
                sep = ";")

dat$CASE |> unique() |> length()

## [1] 357

Exclusion Criteria

# Participant answered all questions meaningfully
dat |> count(SC02)

outlierMeaningful <- dat[dat$SC02 == 2 ,]$CASE


# Exclude participants how have taken less than 15min for the study (indicating that they skipped parts of the videos)
timeTaken <- dat$TIME005
timeTaken <- timeTaken[!is.na(timeTaken)]
lower_bound <- min(timeTaken, 800) - 1e-8
upper_bound <- max(timeTaken, 1500) + 1e-8


breaks <- c(lower_bound, seq(800, 1500, by = 100), upper_bound)

ggplot(data.frame(x = timeTaken), aes(x = x)) +
  geom_histogram(breaks = breaks,
                 closed = "left") +
  geom_vline(xintercept = 900, color = "red", size = 1) +
  coord_cartesian(xlim = c(800, 1500)) +
  labs(x = "Time taken",
       y = "Count")

outlierTime <- dat[dat$TIME005 < 900 ,]$CASE

# remove outlier from dataset
dat <- dat[! dat$CASE %in% outlierMeaningful, ]
dat <- dat[! dat$CASE %in% outlierTime, ]

outlierMeaningful |> unique() |> length()

## [1] 2

outlierTime |> unique() |> length()

## [1] 57

c(outlierMeaningful, outlierTime) |> unique() |> length()

## [1] 58

dat$CASE|> unique() |> length()

## [1] 299

Prepare dataframe

# Refactor columns to factors 
cols <- c("SK01","SK02","SK03","SK04","SK05", "SK06","SK07","SK08","SK09","SK10", "SK12","SK13") 
dat[cols] <- lapply(dat[cols], function(x) as.numeric(x) - 1)

# Label data
label(dat$SK01) <- "Hand lowering"
label(dat$SK02) <- "Moving hands together"
label(dat$SK03) <- "Mosquito hallucination"
label(dat$SK04) <- "Sweet taste"
label(dat$SK05) <- "Sour taste"
label(dat$SK06) <- "Arm rigidity"
label(dat$SK07) <- "Arm immobilization"
label(dat$SK08) <- "Music hallucination"
label(dat$SK09) <- "Negative visual hallucination"
label(dat$SK10) <- "Amnesia"
label(dat$SK12) <- "Urge to press space bar"
label(dat$SK13) <- "Space bar amnesia"

label(dat$KK02_01) <- "Items remembered before being asked again"
label(dat$KK03_01) <- "Items remembered after being asked again"

label(dat$SO01) <- "Gender"
label(dat$SO02_01) <- "Age"

label(dat$SC02) <- "Answered all questions seriously"

label(dat$TIME005) <- "Completion time (experiment duration)"

# calculate item 4 (arithmetic mean of taste sweet and sour) and 10 (geometric mean of Urge to press space bar and space bar amnesia)
dat <- dat |> mutate(SK04SK05 = (SK04 + SK05) / 2)
label(dat$SK04SK05) <- "Taste Halluctionation"

dat <- dat |> mutate(SK12SK13 = sqrt(SK12 * SK13))
label(dat$SK12SK13) <- "Post-session suggestion"

Single item descriptives

#function to get bootstrapped drop item cronbach's alpha

alpha_drop_boot_ci <- function(df, B = 2000, conf = 0.95, seed = NULL, warnings = FALSE, ...) {
  if (!is.null(seed)) {
    set.seed(seed)
  }

  if (!requireNamespace("psych", quietly = TRUE)) {
    stop("Package 'psych' must be installed.")
  }

  if (!is.data.frame(df) && !is.matrix(df)) {
    stop("'df' must be a data frame or matrix containing only the scale items.")
  }

  df <- as.data.frame(df)

  if (ncol(df) < 2) {
    stop("At least 2 items are required.")
  }

  if (!all(vapply(df, is.numeric, logical(1)))) {
    stop("All columns in 'df' must be numeric.")
  }

  if (!is.numeric(B) || length(B) != 1 || B < 1) {
    stop("'B' must be a positive integer.")
  }
  B <- as.integer(B)

  if (!is.numeric(conf) || length(conf) != 1 || conf <= 0 || conf >= 1) {
    stop("'conf' must be a single number between 0 and 1.")
  }

  alpha_level <- 1 - conf
  probs <- c(alpha_level / 2, 1 - alpha_level / 2)

  # Original estimates
  a0 <- psych::alpha(df, warnings = warnings, ...)
  item_names <- rownames(a0$alpha.drop)

  drop_est <- a0$alpha.drop$raw_alpha
  names(drop_est) <- item_names

  # Bootstrap matrix: rows = resamples, cols = items
  boot_drop <- matrix(NA_real_, nrow = B, ncol = length(drop_est))
  colnames(boot_drop) <- item_names

  for (b in seq_len(B)) {
    idx <- sample.int(nrow(df), size = nrow(df), replace = TRUE)

    boot_drop[b, ] <- tryCatch({
      ab <- psych::alpha(df[idx, , drop = FALSE], warnings = warnings, ...)
      vals <- ab$alpha.drop$raw_alpha
      names(vals) <- rownames(ab$alpha.drop)
      vals <- vals[item_names]
      as.numeric(vals)
    }, error = function(e) {
      rep(NA_real_, length(drop_est))
    })
  }

  # Percentile bootstrap CIs
  cis <- t(apply(boot_drop, 2, function(x) {
    if (all(is.na(x))) {
      c(NA_real_, NA_real_)
    } else {
      quantile(x, probs = probs, na.rm = TRUE, names = FALSE)
    }
  }))

  out <- data.frame(
    item_dropped = item_names,
    alpha_if_dropped = as.numeric(drop_est[item_names]),
    ci_lower = cis[, 1],
    ci_upper = cis[, 2],
    row.names = NULL
  )

  out
}

cit_boot_ci <- function(df, B = 2000, conf = 0.95, seed = NULL) {
  if (!is.null(seed)) set.seed(seed)

  if (!is.data.frame(df) && !is.matrix(df)) {
    stop("'df' must be a data frame or matrix containing only the scale items.")
  }

  df <- as.data.frame(df)

  if (ncol(df) < 2) stop("At least 2 items are required.")

  if (!all(vapply(df, is.numeric, logical(1)))) {
    stop("All columns in 'df' must be numeric.")
  }

  if (!is.numeric(B) || length(B) != 1 || B < 1) {
    stop("'B' must be a positive integer.")
  }

  B <- as.integer(B)

  if (!is.numeric(conf) || length(conf) != 1 || conf <= 0 || conf >= 1) {
    stop("'conf' must be a single number between 0 and 1.")
  }

  items <- names(df)

  cit_fun <- function(dat) {
    sapply(items, function(item) {
      other_items <- setdiff(items, item)

      total_other_mean <- rowMeans(
        dat[, other_items, drop = FALSE],
        na.rm = TRUE
      )

      cor(
        dat[[item]],
        total_other_mean,
        use = "pairwise.complete.obs"
      )
    })
  }

  # Original estimates
  cit0 <- cit_fun(df)

  alpha_level <- 1 - conf
  probs <- c(alpha_level / 2, 1 - alpha_level / 2)

  # Bootstrap estimates
  boot_cit <- matrix(NA_real_, nrow = B, ncol = length(items))
  colnames(boot_cit) <- items

  for (b in seq_len(B)) {
    idx <- sample.int(nrow(df), size = nrow(df), replace = TRUE)

    boot_cit[b, ] <- tryCatch({
      cit_fun(df[idx, , drop = FALSE])
    }, error = function(e) {
      rep(NA_real_, length(items))
    })
  }

  # Percentile bootstrap CIs
  cis <- t(apply(boot_cit, 2, function(x) {
    if (all(is.na(x))) {
      c(NA_real_, NA_real_)
    } else {
      quantile(x, probs = probs, na.rm = TRUE, names = FALSE)
    }
  }))

  data.frame(
    item = items,
    corrected_item_total = as.numeric(cit0[items]),
    ci_lower = cis[, 1],
    ci_upper = cis[, 2],
    row.names = NULL
  )
}

# Items to include
items <- c("SK01","SK02","SK03","SK04SK05","SK06","SK07","SK08","SK09","SK10","SK12SK13")
  # c("SK01","SK02","SK03","SK04","SK05",
  #          "SK06","SK07","SK08","SK09","SK10",
  #          "SK12","SK13","SK04SK05","SK12SK13")

# Subset data
dat_items <- dat[, items]

# Cronbach's alpha
alpha_res <- psych::alpha(dat_items, n.iter=5000)

# 95% C for # Cronbach's alpha
alpha_res$boot.ci

##      2.5%       50%     97.5% 
## 0.6697629 0.7246929 0.7682909

# Extract stats
means <- colMeans(dat_items, na.rm = TRUE)
sds   <- apply(dat_items, 2, sd, na.rm = TRUE)

# Skewness and Kurtosis
skew  <- apply(dat_items, 2, psych::skew, na.rm = TRUE)
kurt  <- apply(dat_items, 2, psych::kurtosi, na.rm = TRUE)


# Alpha if item removed
alpha_drop <- alpha_res$alpha.drop

# Bootstrap 95% CI for alpha-if-item-dropped
drop_boot <- alpha_drop_boot_ci(dat_items, B = 5000, seed = 123)
drop_boot

# Calculate PC score per person
PC <- dat |>
  mutate(
    PC = rowMeans(
      across(c(SK01, SK02, SK03,
               SK06, SK07, SK08, SK09, SK10,
               SK04SK05, SK12SK13))
    )
  ) |>
  select(CASE, PC)

# Item-total correlation: each item correlated with PC
item_total <- sapply(items, function(x) {
  cor(dat[[x]], PC$PC, use = "pairwise.complete.obs" , method = "pearson")
})

result <- data.frame(
  item_code  = items,
  Item       = sapply(items, function(x) Hmisc::label(dat[[x]])),
  Mean       = round(means, 2),
  SD         = round(sds, 2),
  Skew       = round(skew, 2),
  Kurt       = round(kurt, 2),
  Alpha      = round(alpha_drop$raw_alpha, 2),
  CI_Low     = round(drop_boot$ci_lower, 2),
  CI_Up      = round(drop_boot$ci_upper, 2),
  ItemTotal  = round(item_total, 2)
)

result 

corrected_item_total <- sapply(items, function(item) {
    other_items <- setdiff(items, item)
  total_other_mean <- rowMeans(dat_items[, other_items], na.rm = TRUE)
  cor(dat_items[[item]], total_other_mean, use = "pairwise.complete.obs")
})

corrected_ItemTotal_CIs <- cit_boot_ci(dat_items, B = 2000, conf = 0.95, seed = 123)

result |> 
  left_join(
    data.frame(
      item_code = items,
      corrected_ItemTotal = round(corrected_item_total, 2),
      row.names = NULL
    )
  ) %>% 
  left_join(corrected_ItemTotal_CIs %>% rename(item_code = item) %>% select(-c(corrected_item_total)) %>% mutate(ci_lower = round(ci_lower, 2), ci_upper = round(ci_upper, 2)))

a0 <- psych::alpha(dat_items, warnings = TRUE)
a0

## 
## Reliability analysis   
## Call: psych::alpha(x = dat_items, warnings = TRUE)
## 
##   raw_alpha std.alpha G6(smc) average_r S/N   ase mean   sd median_r
##       0.73      0.71    0.73       0.2 2.5 0.023  1.8 0.73     0.16
## 
##     95% confidence boundaries 
##          lower alpha upper
## Feldt     0.68  0.73  0.77
## Duhachek  0.68  0.73  0.77
## 
##  Reliability if an item is dropped:
##          raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## SK01          0.68      0.67    0.68      0.18 2.0    0.027 0.017  0.15
## SK02          0.70      0.69    0.70      0.20 2.2    0.025 0.019  0.17
## SK03          0.70      0.68    0.71      0.19 2.2    0.025 0.027  0.15
## SK04SK05      0.71      0.69    0.71      0.20 2.2    0.024 0.028  0.15
## SK06          0.66      0.65    0.67      0.17 1.9    0.029 0.017  0.15
## SK07          0.67      0.66    0.68      0.18 2.0    0.028 0.020  0.15
## SK08          0.72      0.71    0.73      0.22 2.5    0.023 0.026  0.18
## SK09          0.73      0.72    0.74      0.22 2.6    0.023 0.026  0.18
## SK10          0.71      0.69    0.71      0.20 2.3    0.024 0.029  0.14
## SK12SK13      0.74      0.73    0.74      0.23 2.7    0.022 0.024  0.18
## 
##  Item statistics 
##            n raw.r std.r r.cor r.drop mean  sd
## SK01     299  0.65  0.64  0.62   0.53 3.56 1.3
## SK02     299  0.57  0.55  0.50   0.42 3.10 1.5
## SK03     299  0.57  0.56  0.48   0.41 1.54 1.4
## SK04SK05 299  0.54  0.53  0.43   0.37 1.87 1.5
## SK06     299  0.73  0.70  0.70   0.60 2.92 1.6
## SK07     299  0.69  0.66  0.64   0.55 2.52 1.5
## SK08     299  0.38  0.41  0.29   0.23 0.58 1.2
## SK09     299  0.32  0.38  0.23   0.19 0.34 1.0
## SK10     299  0.50  0.52  0.42   0.35 1.23 1.3
## SK12SK13 299  0.32  0.34  0.19   0.14 0.81 1.3

a0$alpha.drop

a0$total

# Calculacte PC Score per person
PC <- dat |> mutate(PC = rowMeans(across(c(SK01,SK02,SK03,
           SK06,SK07,SK08,SK09,SK10,
           SK04SK05,SK12SK13)))) |> 
  select(CASE, PC)

PC

percentile_rank <- function(sample, value) {
  100 * sum(sample <= value) / length(sample)
}

percentile_rank(PC$PC, 2.8)

## [1] 92.97659

Demographics

# age 
mean(dat$SO02_01)

## [1] 31.65217

range(dat$SO02_01)

## [1] 18 72

sd(dat$SO02_01)

## [1] 12.35083

# Gender (1=female, 2=male, 3=diverse)
dat |> count(SO01)

# Time taken
median(dat$TIME005)/60

## Completion time (experiment duration) 
## [1] 16

min(dat$TIME005)/60

## [1] 15.51667

max(dat$TIME005)/60

## [1] 65.35

Plots

# Histogram of PC Scores with Density
ggplot(PC, aes(x = PC)) +
  geom_histogram(aes(y = after_stat(density)),
                 color = "black",
                 fill = "gray",
                 binwidth = 0.5,
                 boundary = 0,
                 alpha = 0.7) +
  geom_density(linewidth = 1.2) +
  labs(x = "PC-Ger Score",
       y = "Density") +
  theme(
    axis.title = element_text(face = "bold", size = 14)
  )

# raincloud plot for PC
my_col <- "#4C72B0"

ggplot(PC, aes(x = PC, y = "")) +
  stat_slab(
    fill = my_col,
    alpha = 0.7,     
    adjust = .5,
    side = "top",
    color = NA
  ) +
  stat_dots(
    fill = my_col,
    alpha = 0.7,     
    side = "bottom",
    justification = 1,
    binwidth = 0.1,
    dotsize = .7,
    stackratio = 0.8
  ) +
  geom_boxplot(
    fill = my_col,
    alpha = 1,       
    width = 0.12,
    outlier.shape = NA,
    color = "black"
  ) +
  labs(
    x = "PC-Ger Score",
    y = "Density"
  ) +
  theme(
    axis.title = element_text(face = "bold", size = 14),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

items <- c("SK01","SK02","SK03","SK04SK05","SK06",
           "SK07","SK08","SK09","SK10","SK12SK13")

var_labels <- c(
  SK01 = "Hand\nlowering",
  SK02 = "Moving\nhands\ntogether",
  SK03 = "Mosquito\nhallucination",
  SK04 = "Sweet\ntaste",
  SK05 = "Sour\ntaste",
  SK06 = "Arm\nrigidity",
  SK07 = "Arm\nimmobilisation",
  SK08 = "Music\nhallucination",
  SK09 = "Negative\nvisual\nhallucination",
  SK10 = "Amnesia",
  SK12 = "Urge to\npress\nspace bar",
  SK13 = "Space bar\namnesia",
  SK04SK05 = "Taste\nhallucination",
  SK12SK13 = "Post-session\nsuggestion"
)

label_numbers <- c(
  SK01 = "1",
  SK02 = "2",
  SK03 = "3",
  SK04SK05 = "4",
  SK06 = "5",
  SK07 = "6",
  SK08 = "7",
  SK09 = "8",
  SK10 = "9",
  SK12SK13 = "10"
)

var_labels_num <- setNames(
  paste0(label_numbers[names(var_labels)], ". ", var_labels),
  names(var_labels)
)

vars_left  <- c("SK01","SK02","SK03","SK04SK05", "SK06")
vars_right <- c("SK07","SK08","SK09","SK10","SK12SK13")

singleItemPlot_data <- dat  |> 
      select(CASE, all_of(items)) %>%
      mutate(across(everything(), ~ as.numeric(.))) %>%
      pivot_longer(cols = -CASE,
               names_to = "Item",
               values_to = "value")

make_plot <- function(data, vars) {
  data_filtered <- data %>%
    filter(Item %in% vars) %>%
    mutate(Item = factor(Item, levels = rev(vars)))

  ggplot() +
        # Layer 2: SK12SK13 with different parameters
    geom_dots(
      data = data_filtered %>% filter(Item == "SK12SK13"),
      aes(x = value, y = Item, fill = Item, color = Item),
      smooth = smooth_discrete(width = 3.5, size = .2), 
      binwidth = 0.1,
      dotsize = .65
    ) +
    
    # Layer 1: all items except SK12SK13
    geom_dots(
      data = data_filtered %>% filter(Item != "SK12SK13"),
      aes(x = value, y = Item, fill = Item, color = Item),
      smooth = smooth_discrete(width = .8, size = .1),
      binwidth = 0.1,
      dotsize = .65
    ) +
    

    scale_y_discrete(labels = var_labels_num, expand = expansion(mult = c(0, 0))) +
    scale_x_continuous(breaks = 0:5) +

    scale_fill_manual(values = c(
      SK01 = "#2ECC71FF",
      SK02 = "#F1C40FFF",
      SK03 = "#E74C3CFF",
      SK04SK05 = "#9B59B6FF",
      SK06 = "#1ABC9CFF",
      SK07 = "#db00b6",
      SK08 = "#D35400FF",
      SK09 = "#65f7f4",
      SK10 = "#1f78b4",
      SK12SK13 = "#666666"
    )) +
    scale_color_manual(values = c(
      SK01 = "#2ECC71FF",
      SK02 = "#F1C40FFF",
      SK03 = "#E74C3CFF",
      SK04SK05 = "#9B59B6FF",
      SK06 = "#1ABC9CFF",
      SK07 = "#db00b6",
      SK08 = "#D35400FF",
      SK09 = "#65f7f4",
      SK10 = "#1f78b4",
      SK12SK13 = "#666666"
    )) +

    labs(
      x = "PC-Ger Score",
      y = NULL
    ) +
    theme_minimal() +
    theme(
      axis.title = element_text(face = "bold", size = 24),
      axis.text = element_text(face = "bold", size = 24),
      axis.text.y = element_text(
        face = "bold",
        size = 24,
        hjust = 0,
        margin = margin(r = 10)
      ),
      legend.position = "none"
    )
}

p1 <- make_plot(singleItemPlot_data, vars_left)
p2 <- make_plot(singleItemPlot_data, vars_right)

p1 + p2

Comparison with Lush et al. (2021), First test data

Lush_dat <- read.csv("PeteOSF/PCS_norms_first_test_data.csv", sep = ",")

Lush_dat <- Lush_dat %>%  
  filter(exclude == 0) %>%             
  filter(Condition == 0) 


rename_map <- c(
  ParticipantNumber                   = "CASE",
  HandLoweringSubjectiveRating        = "SK01",
  HandsTogetherSubjectiveRating       = "SK02",
  MosquitoSubjectiveRating            = "SK03",
  SweetSubRating                      = "SK04",
  SourSubRating                       = "SK05",
  ArmRigiditySubjectiveRating         = "SK06",
  ArmImmobilisationSubjectiveRating   = "SK07",
  MusicSubjectiveRating               = "SK08",
  VisualHallucinationSubjectiveRating = "SK09",
  AmnesiaSubjectiveRating             = "SK10",
  PostHypnoticSub1                    = "SK12",
  PostHypnoticSub2                    = "SK13",
  TasteSubjectiveRating               = "SK04SK05",
  PostHypnoticSubjectiveRating        = "SK12SK13"
)

label_map <- c(
  SK01   = "Hand lowering",
  SK02   = "Moving hands together",
  SK03   = "Mosquito hallucination",
  SK04   = "Sweet taste",
  SK05   = "Sour taste",
  SK06   = "Arm rigidity",
  SK07   = "Arm immobilization",
  SK08   = "Music hallucination",
  SK09   = "Negative visual hallucination",
  SK10   = "Amnesia",
  SK12   = "Urge to press space bar",
  SK13   = "Space bar amnesia",
  SK0405 = "Taste Halluctionation",
  SK1213 = "Post-session suggestion"
)


matched <- intersect(names(rename_map), names(Lush_dat))
names(Lush_dat)[match(matched, names(Lush_dat))] <- rename_map[matched]

# Items to include
items <- c("SK01","SK02","SK03","SK04SK05","SK06","SK07","SK08","SK09","SK10","SK12SK13")
itemsPC <- c("SK01","SK02","SK03","SK04SK05","SK06","SK07","SK08","SK09","SK10","SK12SK13", "PC")

Lush_dat <- Lush_dat %>% 
            select(CASE, all_of(items)) %>% 
            mutate(PC = rowMeans(across(all_of(items), ~ as.numeric(.)))) %>% 
            mutate(group = "Original") %>% 
            mutate(across(all_of(items), ~ as.numeric(.)))

Ger_dat <- dat %>% 
           select(CASE, all_of(items)) %>% 
           mutate(PC = rowMeans(across(all_of(items), ~ as.numeric(.)))) %>% 
           mutate(group = "German") %>% 
           mutate(across(all_of(items), ~ as.numeric(.)))


combined_dat <- bind_rows(Lush_dat, Ger_dat) %>% 
  select(CASE, group, SK01, SK02, SK03, SK04SK05, SK06, SK07, SK08, SK09, SK10, SK12SK13, PC)

for (nm in intersect(names(label_map), names(combined_dat))) {
  Hmisc::label(combined_dat[[nm]]) <- label_map[[nm]]
}


combined_dat_long <- combined_dat %>%
  mutate(across(all_of(itemsPC), ~ as.numeric(haven::zap_labels(.x)))) %>%
  pivot_longer(
    cols = all_of(itemsPC),
    names_to = "Item",
    values_to = "value"
  )

Item-Level Descriptives

item_labels <- sapply(items, function(x) Hmisc::label(combined_dat[[x]]))

item_desc <- combined_dat %>%
  group_by(group) %>%
  group_modify(~{
    
    dat_items <- .x[, items, drop = FALSE]
    
    # Cronbach's alpha
    alpha_res <- psych::alpha(dat_items, n.iter = 5000)
    alpha_drop <- alpha_res$alpha.drop
    
    # Bootstrap CI for alpha if item dropped
    drop_boot <- alpha_drop_boot_ci(dat_items, B = 5000, seed = 123)
    
    # Item descriptives
    means <- colMeans(dat_items, na.rm = TRUE)
    sds   <- apply(dat_items, 2, sd, na.rm = TRUE)
    skew  <- apply(dat_items, 2, psych::skew, na.rm = TRUE)
    kurt  <- apply(dat_items, 2, psych::kurtosi, na.rm = TRUE)
    
    # Person composite score within group
    PC <- .x %>%
      mutate(
        PC = rowMeans(across(all_of(items)), na.rm = TRUE)
      ) %>%
      pull(PC)
    
    # Item-total correlation within group
    item_total <- sapply(items, function(x) {
      cor(.x[[x]], PC, use = "pairwise.complete.obs", method = "pearson")
    })
    
    data.frame(
      item_code = items,
      Item      = unname(item_labels),
      Mean      = round(means, 2),
      SD        = round(sds, 2),
      Skew      = round(skew, 2),
      Kurt      = round(kurt, 2),
      Alpha     = round(alpha_drop$raw_alpha, 2),
      CI_Low    = round(drop_boot$ci_lower, 2),
      CI_Up     = round(drop_boot$ci_upper, 2),
      ItemTotal = round(item_total, 2)
    )
  }) %>%
  ungroup()



corrected_by_group <- combined_dat %>%
  group_by(group) %>%
  group_modify(~ {
    
    dat <- .x
    
    vals <- sapply(items, function(item) {
      
      other_items <- setdiff(items, item)
      total_other_mean <- rowMeans(dat[, other_items], na.rm = TRUE)
      
      cor(dat[[item]],
          total_other_mean,
          use = "pairwise.complete.obs")
    })
    
    data.frame(
      item_code = items,
      corrected_item_total = round(vals, 2)
    )
  }) %>%
  ungroup()

corrected_by_group_CIs <- combined_dat %>% 
  group_by(group) %>% 
  group_modify(~{
        dat_items <- .x[, items, drop = FALSE]

        CIs <- cit_boot_ci(dat_items, B = 2000, conf = 0.95, seed = 123)
  }) %>%
  ungroup()


item_desc |> left_join(
  corrected_by_group
) %>% 
  left_join(corrected_by_group_CIs %>% rename(item_code = item) %>% select(-c(corrected_item_total)) %>% mutate(ci_lower = round(ci_lower, 2), ci_upper = round(ci_upper, 2)))

ks_results <- combined_dat_long |>
  group_by(Item) |>
  summarise(
    ks_stat   = ks.test(
      value[group == "Original"],
      value[group == "German"]
    )$statistic,
    p_value   = ks.test(
      value[group == "Original"],
      value[group == "German"]
    )$p.value,
    .groups = "drop"
  ) |>
  mutate(significant = p_value < .05)

ks_results

label_map <- c(
  SK01      = "1) Hand lowering",
  SK02      = "2) Moving hands together",
  SK03      = "3) Mosquito hallucination",
  SK04SK05  = "4) Taste Hallucination",
  SK06      = "5) Arm rigidity",
  SK07      = "6) Arm immobilization",
  SK08      = "7) Music hallucination",
  SK09      = "8) Negative visual hallucination",
  SK10      = "9) Amnesia",
  SK12SK13  = "10) Post-session suggestion"
)

create_hist_plot <- function(var_name) {
  
  bins <- if (var_name %in% c("SK04SK05", "SK12SK13")) 11 else 6
  
  ggplot(
    combined_dat %>%
      filter(group %in% c("German", "Original")),
    aes(
      x = .data[[var_name]],
      color = group,
      fill = group
    )
  ) +
    geom_histogram(
      aes(y = after_stat(density)),
      alpha = 0.2,
      position = "identity",
      bins = bins
    ) +
    labs(
      title = label_map[[var_name]],
      x = NULL,
      y = "Density",
      color = "Dataset",
      fill = "Dataset"
    ) +
    scale_fill_manual(values = c(
      German = "orange",
      Original = "royalblue"
    )) +
    scale_color_manual(values = c(
      German = "orange",
      Original = "royalblue"
    )) +
    theme_minimal() +
    theme(
      plot.title = element_text(face = "bold", size = 20, hjust = 0.5),
      axis.title = element_text(face = "bold", size = 24),
      axis.text = element_text(face = "bold", size = 20),
      axis.text.y = element_text(
        face = "bold",
        size = 20,
        hjust = 0,
        margin = margin(r = 10)
      ),
      legend.title = element_text(size = 20),
      legend.text = element_text(size = 18)
    )
}

plot_order <- c(
  "SK01", "SK07",
  "SK02", "SK08",
  "SK03", "SK09",
  "SK04SK05", "SK10",
  "SK06", "SK12SK13"
)

hist_plots <- map(plot_order, create_hist_plot)
names(hist_plots) <- plot_order

combined_hist_plot <- wrap_plots(hist_plots, ncol = 2, nrow = 5, guides = "collect") &
  theme(legend.position = "bottom")

combined_hist_plot

Scale-Level (PC score) Descriptives

inner_join(combined_dat %>%
  group_by(group) %>%
  summarise(
    mean = mean(PC, na.rm = TRUE),
    sd   = sd(PC, na.rm = TRUE),
    min  = min(PC, na.rm = TRUE),
    max  = max(PC, na.rm = TRUE),
    skew = psych::skew(PC, na.rm = TRUE),
    kurt = psych::kurtosi(PC, na.rm = TRUE),
    p10  = quantile(PC, probs = 0.10, na.rm = TRUE),
    p25  = quantile(PC, probs = 0.25, na.rm = TRUE),
    p50  = quantile(PC, probs = 0.50, na.rm = TRUE),
    p75  = quantile(PC, probs = 0.75, na.rm = TRUE),
    p90  = quantile(PC, probs = 0.90, na.rm = TRUE),
  ),
item_desc %>% 
  group_by(group) %>%
  summarise(meanItemTotal = mean(ItemTotal))
) |> select(group, p10, p25, p50, p75, p90)

res <- ks.test(
  combined_dat_long |>
    filter(group == "Original", Item == "PC") |>
    pull(value),
  combined_dat_long |>
    filter(group == "German", Item == "PC") |>
    pull(value)
)

res

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  pull(filter(combined_dat_long, group == "Original", Item == "PC"), value) and pull(filter(combined_dat_long, group == "German", Item == "PC"), value)
## D = 0.089533, p-value = 0.2317
## alternative hypothesis: two-sided

library(transport)

wasserstein1d(
    combined_dat_long |>
    filter(group == "Original", Item == "PC") |>
    pull(value),
  combined_dat_long |>
    filter(group == "German", Item == "PC") |>
    pull(value)

)

## [1] 0.1008361

group_colors <- c(
  German = "orange",
  Original = "royalblue"
)


percentile_points <- combined_dat |>
  dplyr::group_by(group) |>
  dplyr::summarise(
    `10th percentile` = quantile(PC, 0.10, na.rm = TRUE),
    `90th percentile` = quantile(PC, 0.90, na.rm = TRUE),
    .groups = "drop"
  ) |>
  tidyr::pivot_longer(
    cols = c(`10th percentile`, `90th percentile`),
    names_to = "percentile",
    values_to = "PC"
  )

boxplot_top <- ggplot(
  combined_dat,
  aes(x = PC, y = group, color = group, fill = group)
) +
  geom_boxplot(
    alpha = 0.25,
    linewidth = 1,
    width = 0.6,
    outlier.shape = NA
  ) +
  geom_point(
    data = percentile_points,
    aes(
      x = PC,
      y = group,
      color = group,
      shape = percentile
    ),
    inherit.aes = FALSE,
    size = 6,
    stroke = 1.2
  ) +
  scale_fill_manual(values = group_colors) +
  scale_color_manual(values = group_colors) +
  scale_shape_manual(
    values = c(
      "10th percentile" = 17,
      "90th percentile" = 16
    )
  ) +
  labs(x = NULL, y = NULL) +
  theme_minimal() +
  theme(
    axis.title = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    panel.grid = element_blank(),
    legend.position = "none",
    plot.margin = margin(5, 5, 0, 5)
  )

histo <- ggplot(
  combined_dat,
  aes(x = PC, color = group, fill = group)
) +
  geom_histogram(
    aes(y = after_stat(density)),
    alpha = 0.2,
    position = "identity",
    bins = 30
  ) +
  geom_density(
    linewidth = 1.2,
    alpha = 0.2
  ) +
  labs(
    x = "PC Score",
    y = "Density",
    color = "Dataset",
    fill = "Dataset"
  ) +
  scale_fill_manual(values = group_colors) +
  scale_color_manual(values = group_colors) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold", size = 24),
    axis.text = element_text(face = "bold", size = 24),
    axis.text.y = element_text(
      face = "bold",
      size = 24,
      hjust = 0,
      margin = margin(r = 10)
    ),
    legend.title = element_text(size = 24),
    legend.text = element_text(size = 24),
    plot.margin = margin(0, 5, 5, 5)
  )

# Split the two groups
german <- combined_dat %>%
  filter(group == "German") %>%
  pull(PC)

original <- combined_dat %>%
  filter(group == "Original") %>%
  pull(PC)

# Create common probability sequence
n <- min(length(german), length(original))
probs <- seq(0, 1, length.out = n)

# Compute quantiles
qq_data <- data.frame(
  German = quantile(german, probs = probs, na.rm = TRUE),
  Original = quantile(original, probs = probs, na.rm = TRUE)
)

# Plot
qq1 <- ggplot(qq_data, aes(x = Original, y = German)) +
  geom_point() +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
  labs(
    x = "Original PC quantiles",
    y = "German PC quantiles"
  ) +
  coord_flip() +
  theme_minimal() +
    theme(
      axis.title = element_text(face = "bold", size = 24),
      axis.text = element_text(face = "bold", size = 24),
      axis.text.y = element_text(
        face = "bold",
        size = 24,
        hjust = 0,
        margin = margin(r = 10)
      ),
      legend.position = "none"
    )


# Q-Q plot against normal distribution
qq2 <- ggplot(data.frame(PC = german), aes(sample = PC)) +
  stat_qq() +
  stat_qq_line(linetype = "dashed", color = "red") +
  labs(
    x = "Normal distribution quantiles",
    y = "German PC quantiles"
  ) +
  coord_flip() +
  theme_minimal() +
    theme(
      axis.title = element_text(face = "bold", size = 24),
      axis.text = element_text(face = "bold", size = 24),
      axis.text.y = element_text(
        face = "bold",
        size = 24,
        hjust = 0,
        margin = margin(r = 10)
      ),
      legend.position = "none"
    )

qq1

qq2

Compare Cronbach’s Alpha

alpha_lush <- psych::alpha(Lush_dat[items], n.iter=5000)
alpha_ger  <- psych::alpha(Ger_dat[items], n.iter=5000)

  

data.frame(
  Dataset = c("Original", "German"),
  Alpha   = c(alpha_lush$total$raw_alpha,
              alpha_ger$total$raw_alpha)
)

alpha_lush$boot.ci

##      2.5%       50%     97.5% 
## 0.5981377 0.6795979 0.7396602

alpha_ger$boot.ci

##      2.5%       50%     97.5% 
## 0.6685916 0.7250248 0.7690652

Exploratory Factor Analysis

efa_vars <- combined_dat %>%
  select(SK01:SK12SK13) %>%
  names()

eigenvalues <- combined_dat %>%
  group_by(group) %>%
  group_modify(~{
    cor_mat <- .x %>%
      select(all_of(efa_vars)) %>%
      cor(use = "pairwise.complete.obs")

    tibble(
      Factor = seq_along(eigen(cor_mat)$values),
      Eigenvalue = eigen(cor_mat)$values
    )
  })

scree <- eigenvalues %>%
  ggplot(
  aes(x = Factor, y = Eigenvalue, color = group)
) +
  geom_line() +
  geom_point(size = 2) +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  labs(
    x = "Factor",
    y = "Eigenvalue",
    color = "Dataset",
  ) +
  scale_color_manual(values = group_colors) +
  theme_minimal() +
  theme(
    axis.title = element_text(face = "bold", size = 24),
    axis.text = element_text(face = "bold", size = 24),
    axis.text.y = element_text(
      face = "bold",
      size = 24,
      hjust = 0,
      margin = margin(r = 10)
    ),
    legend.title = element_text(size = 24),
    legend.text = element_text(size = 24),
    plot.margin = margin(0, 5, 5, 5)
  )




efa_result_german <- fa(
  r = combined_dat %>%
    filter(group == "German")|> 
    select(all_of(efa_vars)),
  nfactors = 2,
  rotate = "oblimin")

efa_result_original <- fa(
  r = combined_dat %>%
    filter(group == "Original")|> 
    select(all_of(efa_vars)),
  nfactors = 2,
  rotate = "oblimin")

factor.congruence(
  efa_result_german,
  efa_result_original
)

##      MR1  MR2
## MR1 0.93 0.15
## MR2 0.28 0.86

efa_result_german$Vaccounted

##                             MR1       MR2
## SS loadings           2.1720848 1.0842129
## Proportion Var        0.2172085 0.1084213
## Cumulative Var        0.2172085 0.3256298
## Proportion Explained  0.6670412 0.3329588
## Cumulative Proportion 0.6670412 1.0000000

efa_result_original$Vaccounted

##                             MR1       MR2
## SS loadings           1.6488453 1.0728513
## Proportion Var        0.1648845 0.1072851
## Cumulative Var        0.1648845 0.2721697
## Proportion Explained  0.6058152 0.3941848
## Cumulative Proportion 0.6058152 1.0000000

ggsave(
  filename = "EFA_scree.png",
  plot = scree,
  width = 20,
  height = 10,
  units = "in",
  dpi = 300
)