Study 1

Participants

library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)

df <- haven::read_sav("../../data/Murphy2020/Study 1.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other")))) 
  
paste0(
  "Data from the [study 1](https://osf.io/3m5nh/?view_only=a68051df4abe4ecb992f22dc8c17f769) (Murphy et al., 2020), downloaded from OSF, included ",
  report::report_participants(df, age = "Age", sex = NA, gender = "Gender"),
  "."
)

[1] “Data from the study 1 (Murphy et al., 2020), downloaded from OSF, included 451 participants (Mean age = 25.8, SD = 8.4, range: [18, 69]; Gender: 69.4% women, 29.5% men, 1.11% non-binary).”

Descriptive

data <- select(df, -Age, -Gender)

dens1 <- estimate_density(data, method="kernSmooth")

plot(dens1) 


data <- standardize(data)

EFA

n <- parameters::n_factors(data, n_max = 15)

n
## # Method Agreement Procedure:
## 
## The choice of 1 dimensions is supported by 5 (31.25%) methods out of 16 (Bentler, Acceleration factor, Scree (R2), VSS complexity 1, Velicer's MAP).

knitr::kable(as.data.frame(n)[n$n_Factors %in% c(1, 4), ])
n_Factors Method Family
1 1 Bentler Bentler
2 1 Acceleration factor Scree
3 1 Scree (R2) Scree_SE
4 1 VSS complexity 1 VSS
5 1 Velicer’s MAP Velicers_MAP
10 4 beta Multiple_regression
11 4 Optimal coordinates Scree
12 4 Parallel analysis Scree
13 4 Kaiser criterion Scree 
plot(n)

efa1 <- parameters::factor_analysis(data, n=1, sort = TRUE)
efa1
## # Loadings from Factor Analysis (no rotation)
## 
## Variable        | MR1  | Complexity | Uniqueness
## ------------------------------------------------
## Pain            | 0.60 |    1.00    |    0.64   
## Sneeze          | 0.59 |    1.00    |    0.65   
## Temp            | 0.57 |    1.00    |    0.67   
## Wind            | 0.57 |    1.00    |    0.68   
## Affective_touch | 0.56 |    1.00    |    0.69   
## Muscles         | 0.55 |    1.00    |    0.69   
## Tickle          | 0.54 |    1.00    |    0.70   
## Itch            | 0.54 |    1.00    |    0.71   
## Defecate        | 0.54 |    1.00    |    0.71   
## Vomit           | 0.54 |    1.00    |    0.71   
## Taste           | 0.53 |    1.00    |    0.72   
## Cough           | 0.53 |    1.00    |    0.72   
## Bruise          | 0.52 |    1.00    |    0.73   
## Burp            | 0.52 |    1.00    |    0.73   
## Urinate         | 0.50 |    1.00    |    0.75   
## Breathing       | 0.50 |    1.00    |    0.75   
## Thirsty         | 0.44 |    1.00    |    0.81   
## Sex_arousal     | 0.41 |    1.00    |    0.83   
## Hungry          | 0.41 |    1.00    |    0.84   
## Heart           | 0.39 |    1.00    |    0.85   
## Blood_Sugar     | 0.38 |    1.00    |    0.86   
## 
## The unique latent factor accounted for 26.51% of the total variance of the original data.

efa4 <- parameters::factor_analysis(data, n=4, rotation = "oblimin", sort = TRUE)
efa4
## # Rotated loadings from Factor Analysis (oblimin-rotation)
## 
## Variable        |    MR2    |  MR3  |    MR1    |    MR4    | Complexity | Uniqueness
## -------------------------------------------------------------------------------------
## Tickle          |   0.72    | 0.03  |   -0.10   |   0.11    |    1.09    |    0.47   
## Itch            |   0.67    | 0.14  | -9.44e-03 |   -0.10   |    1.13    |    0.46   
## Affective_touch |   0.49    | -0.13 |   0.29    |   0.13    |    1.97    |    0.59   
## Blood_Sugar     |   0.29    | 0.21  |   0.18    |   -0.23   |    3.51    |    0.75   
## Burp            |   -0.03   | 0.73  |   -0.03   | -1.21e-03 |    1.00    |    0.50   
## Cough           | 7.56e-03  | 0.69  | 7.00e-03  |   -0.01   |    1.00    |    0.52   
## Sneeze          |   0.07    | 0.54  | 7.09e-03  |   0.20    |    1.30    |    0.55   
## Wind            |   0.08    | 0.46  |   0.07    |   0.16    |    1.37    |    0.62   
## Hungry          |   -0.16   | 0.04  |   0.65    |   0.03    |    1.13    |    0.60   
## Thirsty         | -4.84e-03 | 0.01  |   0.49    |   0.11    |    1.11    |    0.70   
## Breathing       |   0.14    | 0.05  |   0.44    |   0.03    |    1.24    |    0.70   
## Heart           |   0.11    | -0.07 |   0.40    |   0.08    |    1.32    |    0.78   
## Bruise          |   0.31    | 0.25  |   0.34    |   -0.25   |    3.74    |    0.59   
## Pain            |   0.29    | 0.01  |   0.32    |   0.20    |    2.67    |    0.61   
## Temp            |   0.19    | 0.13  |   0.29    |   0.16    |    2.91    |    0.67   
## Muscles         |   0.20    | 0.13  |   0.26    |   0.16    |    3.24    |    0.69   
## Sex_arousal     |   0.18    | 0.02  |   0.21    |   0.15    |    2.81    |    0.82   
## Defecate        |   0.10    | 0.11  | 2.77e-03  |   0.66    |    1.10    |    0.46   
## Urinate         |   -0.11   | 0.13  |   0.24    |   0.52    |    1.65    |    0.54   
## Taste           |   0.23    | 0.06  |   0.17    |   0.30    |    2.60    |    0.69   
## Vomit           |   0.21    | 0.18  |   0.11    |   0.25    |    3.26    |    0.71   
## 
## The 4 latent factors (oblimin rotation) accounted for 37.94% of the total variance of the original data (MR2 = 10.46%, MR3 = 10.17%, MR1 = 10.08%, MR4 = 7.23%).

plot(efa4)

CFA

cfa1 <-  parameters::efa_to_cfa(efa1, names="IAS") |> 
  lavaan::cfa(data=data)

cfa4 <- parameters::efa_to_cfa(efa4, names=c("Skin", "Expulsion", "Interoception", "Elimination"), threshold = "max") |> 
  lavaan::cfa(data=data)

anova(cfa1, cfa4)
## 
## Chi-Squared Difference Test
## 
##       Df   AIC   BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)    
## cfa4 183 24748 24945   545                                        
## cfa1 189 24968 25141   777        232 0.289       6     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cfa4mod <- "
Skin =~ Affective_touch + Tickle + Itch
Expulsion =~ Sneeze + Cough + Wind + Burp
Interoception =~ Heart + Hungry + Breathing + Thirsty + Temp + Sex_arousal + Muscles + Bruise + Pain
Elimination =~ Urinate + Defecate + Vomit
" |> 
lavaan::cfa(data=data)

cfa5 <- "
Skin =~ Affective_touch + Tickle + Itch
Expulsion =~ Sneeze + Cough + Wind + Burp
Interoception =~ Heart + Hungry + Breathing + Thirsty + Temp + Sex_arousal
Nociception =~ Muscles + Bruise + Pain
Elimination =~ Urinate + Defecate + Vomit
" |> 
lavaan::cfa(data=data)

anova(cfa4mod, cfa5)
## 
## Chi-Squared Difference Test
## 
##          Df   AIC   BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)    
## cfa5    142 22362 22559   441                                        
## cfa4mod 146 22383 22563   469       28.8 0.117       4    8.4e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cfa5mod <- "
Skin =~ Affective_touch + Tickle + Itch
Expulsion =~ Sneeze + Wind + Burp
Interoception =~ Heart + Hungry + Breathing + Thirsty + Temp + Sex_arousal
Nociception =~ Muscles + Bruise + Pain
Elimination =~ Urinate + Defecate + Vomit
" |> 
lavaan::cfa(data=data)

anova(cfa5, cfa5mod)
## 
## Chi-Squared Difference Test
## 
##          Df   AIC   BIC Chisq Chisq diff  RMSEA Df diff Pr(>Chisq)    
## cfa5mod 125 21256 21446   379                                         
## cfa5    142 22362 22559   441       61.4 0.0761      17    6.1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cfa5h <- "
Skin =~ Affective_touch + Tickle + Itch
Expulsion =~ Sneeze + Wind + Burp
Nociception =~ Muscles + Bruise + Pain
Elimination =~ Urinate + Defecate + Vomit
Interoception =~ Heart + Hungry + Breathing + Thirsty + Temp + Sex_arousal + Skin + Expulsion + Nociception + Elimination
" |> 
lavaan::cfa(data=data)

anova(cfa5mod, cfa5h)
## 
## Chi-Squared Difference Test
## 
##          Df   AIC   BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)    
## cfa5mod 125 21256 21446   379                                        
## cfa5h   131 21288 21452   422       43.2 0.117       6      1e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
check_perf <- function(mod) {
  perf <- model_performance(mod, metrics = c("RMSEA", "CFI", "SRMR", "AGFI", "NNFI"))
  interpret(perf) |> 
    data_rename(pattern = c("Name", "Value"), replacement = c("Index", "Score")) |> 
    arrange(Index)
}

knitr::kable(check_perf(cfa5mod))
Index Score Threshold Interpretation
AGFI 0.882 0.90 poor
CFI 0.878 0.90 poor
NNFI 0.851 0.90 poor
RMSEA 0.067 0.05 poor
SRMR 0.057 0.08 satisfactory

New sample

df2 <- haven::read_sav("../../data/Murphy2020/Study 6 IAS.sav") |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(Gender == 1, "Male", ifelse(Gender == 2, "Female", "Other")))) |> 
  select(-ESL)
  
paste0(
  "Data from the [study 6](https://osf.io/3m5nh/?view_only=a68051df4abe4ecb992f22dc8c17f769) (Murphy et al., 2020), downloaded from OSF, included ",
  report::report_participants(df2, age = "Age", sex = NA, gender = "Gender"),
  "."
)

[1] “Data from the study 6 (Murphy et al., 2020), downloaded from OSF, included 375 participants (Mean age = 35.3, SD = 16.9, range: [18, 91]; Gender: 70.1% women, 28.5% men, 1.33% non-binary).”

Descriptive

data <- select(df2, -Age, -Gender)

dens2 <- estimate_density(data, method="kernSmooth")

plot(dens2) 


data <- standardize(data)

CFA

efa4b <- parameters::factor_analysis(data, n=4, rotation = "oblimin", sort = TRUE)
cfa1_refit <- "IAS =~ Pain + Sneeze + Temp + Wind + Affective_touch + Muscles + Tickle + Itch + Defecate + Vomit + Taste + Cough + Bruise + Burp + Urinate + Breathing + Thirsty + Sex_arousal + Hungry + Heart + Blood_Sugar" |> lavaan::cfa(data=data)
cfa4mod_refit <- update(cfa4mod, data=data)
cfa5_refit <- update(cfa5, data=data)
cfa5mod_refit <- update(cfa5mod, data=data)
cfa5h_refit <- update(cfa5h, data=data)

table <- anova(cfa5mod_refit, cfa5h_refit, cfa5_refit, cfa4mod_refit, cfa1_refit)
table
## 
## Chi-Squared Difference Test
## 
##                Df   AIC   BIC Chisq Chisq diff  RMSEA Df diff Pr(>Chisq)    
## cfa5mod_refit 125 17135 17316   339                                         
## cfa5h_refit   131 17166 17323   382       42.9 0.1281       6    1.2e-07 ***
## cfa5_refit    142 18068 18257   412       29.8 0.0676      11     0.0017 ** 
## cfa4mod_refit 146 18092 18265   444       32.1 0.1369       4    1.8e-06 ***
## cfa1_refit    189 20256 20421   757      312.5 0.1293      43    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(check_perf(cfa5mod_refit))
Index Score Threshold Interpretation
AGFI 0.880 0.90 poor
CFI 0.906 0.90 satisfactory
NNFI 0.885 0.90 poor
RMSEA 0.068 0.05 poor
SRMR 0.063 0.08 satisfactory

Study 2

Participants

df3 <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  select(Participant, Sex, Age, starts_with("Item_IAS")) |> 
  filter(!is.na(Item_IAS_Elimination_2))

names(df3) <- str_remove_all(names(df3), "Item_IAS_")
names(df3) <- case_when(
  names(df3) == "Interoception_1" ~ "Heart",
  names(df3) == "Interoception_2" ~ "Hungry",
  names(df3) == "Interoception_3" ~ "Breathing",
  names(df3) == "Interoception_4" ~ "Thirsty",
  names(df3) == "Interoception_5" ~ "Temp",
  names(df3) == "Interoception_6" ~ "Sex_arousal",
  names(df3) == "Elimination_1" ~ "Urinate",
  names(df3) == "Elimination_2" ~ "Defecate",
  names(df3) == "Elimination_3" ~ "Vomit",
  names(df3) == "Expulsion_1" ~ "Wind",
  names(df3) == "Expulsion_2" ~ "Burp",
  names(df3) == "Expulsion_3" ~ "Sneeze",
  names(df3) == "Nociception_1" ~ "Muscles",
  names(df3) == "Nociception_2" ~ "Bruise",
  names(df3) == "Nociception_3" ~ "Pain",
  names(df3) == "Skin_1" ~ "Affective_touch",
  names(df3) == "Skin_2" ~ "Tickle",
  names(df3) == "Skin_3" ~ "Itch",
  TRUE ~ names(df3)
  )

report::report_participants(df3, age = "Age", sex = "Sex")

[1] “485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Sex: 50.3% females, 49.7% males, 0.0% other)”

Descriptive

data <- select(df3, -Age, -Sex, -Participant)

dens3 <- estimate_density(data, method="kernSmooth")

plot(dens3) +
  theme_modern()



data[] <- standardize(data) |> 
  sapply(as.numeric)

EFA

n <- parameters::n_factors(data, n_max = 15)

n
## # Method Agreement Procedure:
## 
## The choice of 4 dimensions is supported by 6 (31.58%) methods out of 19 (beta, Optimal coordinates, Parallel analysis, Kaiser criterion, BIC, BIC (adjusted)).
plot(n)

efa4c <- parameters::factor_analysis(data, n=4, rotation = "oblimin", sort = TRUE, fm="ml")
efa4c
## # Rotated loadings from Factor Analysis (oblimin-rotation)
## 
## Variable        |   ML1    |    ML3    |    ML4    |   ML2    | Complexity | Uniqueness
## ---------------------------------------------------------------------------------------
## Hungry          |   0.64   |   -0.05   |   0.14    |   0.05   |    1.12    |    0.50   
## Thirsty         |   0.57   |   0.09    |   0.06    |   0.03   |    1.08    |    0.56   
## Sex_arousal     |   0.56   |   0.08    |   0.02    |   0.03   |    1.04    |    0.62   
## Pain            |   0.55   |   0.05    |   0.09    |  -0.01   |    1.07    |    0.61   
## Temp            |   0.36   |   0.11    |   0.19    |   0.07   |    1.83    |    0.65   
## Affective_touch |   0.26   |   -0.03   |   0.14    |   0.25   |    2.59    |    0.76   
## Defecate        |   0.06   |   0.67    | -5.79e-03 |   0.03   |    1.02    |    0.49   
## Sneeze          |  -0.18   |   0.53    |   0.29    |   0.05   |    1.84    |    0.59   
## Vomit           |   0.13   |   0.53    |   -0.05   |   0.08   |    1.17    |    0.61   
## Urinate         |   0.38   |   0.48    |   -0.11   |  -0.02   |    2.04    |    0.53   
## Burp            | 7.02e-03 |   0.45    |   0.17    |   0.09   |    1.36    |    0.64   
## Breathing       |   0.02   |   0.02    |   0.72    | 8.85e-03 |    1.00    |    0.45   
## Heart           |   0.17   | -1.25e-03 |   0.60    |   0.01   |    1.17    |    0.50   
## Muscles         |   0.15   |   0.07    |   0.47    |   0.07   |    1.31    |    0.60   
## Itch            | 1.25e-03 |   -0.03   |   -0.03   |   0.78   |    1.00    |    0.42   
## Tickle          |  -0.02   |   0.03    |   -0.01   |   0.73   |    1.00    |    0.47   
## Bruise          | 8.31e-03 |   0.04    |   0.07    |   0.37   |    1.10    |    0.81   
## Wind            |   0.09   |   0.08    |   0.11    |   0.25   |    1.92    |    0.83   
## 
## The 4 latent factors (oblimin rotation) accounted for 40.86% of the total variance of the original data (ML1 = 12.30%, ML3 = 10.18%, ML4 = 9.32%, ML2 = 9.06%).

plot(efa4c)

##CFA

cfa1 <- paste("Interoception =~ ", paste(names(data), collapse = " + ")) |> 
  lavaan::cfa(data=data)

cfa5 <- "
Skin =~ Affective_touch + Tickle + Itch
Expulsion =~ Wind + Burp + Sneeze
Interoception =~ Heart + Hungry + Breathing + Thirsty + Temp + Sex_arousal
Nociception =~ Muscles + Bruise + Pain
Elimination =~ Urinate + Defecate + Vomit
" |> 
  lavaan::cfa(data=data)

cfa4 <- "
Homeostasis =~ Hungry + Thirsty + Sex_arousal + Pain + Temp + Affective_touch + Urinate
Skin =~ Affective_touch + Tickle + Itch + Bruise + Wind
Elimination =~ Defecate + Sneeze + Vomit + Urinate + Burp
Muscle =~ Breathing + Heart + Muscles
" |> 
  lavaan::cfa(data=data)

anova(cfa1, cfa5, cfa4)
## 
## Chi-Squared Difference Test
## 
##       Df   AIC   BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)    
## cfa5 125 22718 22911   353                                        
## cfa4 127 22574 22758   213       -140   0.0       2          1    
## cfa1 135 22914 23065   570        357   0.3       8     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


cfa4mod <- "
Homeostasis =~ Hungry + Thirsty + Sex_arousal 
Skin =~ Tickle + Itch + Bruise 
Elimination =~ Defecate + Sneeze + Vomit 
Muscle =~ Breathing + Heart + Muscles
" |> 
  lavaan::cfa(data=data)

cfa5short <- "
Skin =~ Tickle + Itch
Expulsion =~ Sneeze
Interoception =~ Heart + Hungry + Breathing + Thirsty + Sex_arousal
Nociception =~ Muscles + Bruise
Elimination =~ Defecate + Vomit
" 

anova(cfa4mod, lavaan::cfa(cfa5short, data=data))
## 
## Chi-Squared Difference Test
## 
##                                     Df   AIC   BIC Chisq Chisq diff RMSEA
## lavaan::cfa(cfa5short, data = data) 45 15231 15369 158.5                 
## cfa4mod                             48 15140 15266  73.4      -85.1     0
##                                     Df diff Pr(>Chisq)
## lavaan::cfa(cfa5short, data = data)                   
## cfa4mod                                   3          1

cfa4h <- "
Homeostasis =~ Hungry + Thirsty + Sex_arousal 
Skin =~ Tickle + Itch + Bruise 
Elimination =~ Defecate + Sneeze + Vomit 
Muscle =~ Breathing + Heart + Muscles
Interoception =~ Homeostasis + Skin + Elimination + Muscle
" |> 
  lavaan::cfa(data=data)

anova(cfa4mod, cfa4h)
## 
## Chi-Squared Difference Test
## 
##         Df   AIC   BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)  
## cfa4mod 48 15140 15266  73.4                                      
## cfa4h   50 15142 15259  79.3       5.85 0.063       2      0.054 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_perf(cfa4mod)
##   Index Score Threshold Interpretation
## 1  AGFI 0.961      0.90   satisfactory
## 2   CFI 0.982      0.90   satisfactory
## 3  NNFI 0.976      0.90   satisfactory
## 4 RMSEA 0.033      0.05   satisfactory
## 5  SRMR 0.033      0.08   satisfactory

check_perf(cfa4h)
##   Index  Score Threshold Interpretation
## 1  AGFI 0.9596      0.90   satisfactory
## 2   CFI 0.9796      0.90   satisfactory
## 3  NNFI 0.9730      0.90   satisfactory
## 4 RMSEA 0.0347      0.05   satisfactory
## 5  SRMR 0.0364      0.08   satisfactory
check_perf(update(cfa4h, data=standardize(df)))
##   Index  Score Threshold Interpretation
## 1  AGFI 0.9162      0.90   satisfactory
## 2   CFI 0.9020      0.90   satisfactory
## 3  NNFI 0.8706      0.90           poor
## 4 RMSEA 0.0677      0.05           poor
## 5  SRMR 0.0556      0.08   satisfactory
check_perf(update(cfa4h, data=standardize(df2)))
##   Index  Score Threshold Interpretation
## 1  AGFI 0.9178      0.90   satisfactory
## 2   CFI 0.9315      0.90   satisfactory
## 3  NNFI 0.9096      0.90   satisfactory
## 4 RMSEA 0.0660      0.05           poor
## 5  SRMR 0.0504      0.08   satisfactory
graph_data <- tidySEM::prepare_graph(cfa4h, 
                                     variance_diameter=NA) 

nodes <- graph_data$nodes |> 
  mutate(fill = ifelse(shape == "rect", "grey", "#F06292"),
         textsize = ifelse(shape == "rect", rel(3.5), rel(4)),
         fontface = ifelse(shape == "oval", "bold", "plain"),
         hjust = ifelse(shape == "oval", 0.5, -0.3),
         textcolor = ifelse(shape == "oval", "white", "black"),
         label = str_replace(label, "Sex_arousal", "Arousal")) |> 
  select(-x, -y) 

edges <- graph_data$edges |>
  mutate(est = as.numeric(est))

p_graph <- tidygraph::tbl_graph(nodes = nodes, edges = edges) |> 
  ggraph(layout = "kk") +
  geom_edge_arc(aes(filter=lhs!="Interoception",
                    edge_width = est,
                    label=est_std),
                strength = 0.15,
                color="#2196F3",
                angle_calc="along",
                label_dodge=unit(-0.02, "npc"),
                label_size = rel(4)) +
  geom_edge_link(aes(filter=lhs=="Interoception",
                     edge_width = est, 
                     # color = color,
                     label=est_std), 
                 color="#3F51B5",
                 angle_calc="along", 
                 label_dodge=unit(-0.02, "npc"),
                 label_size = rel(5)) +
  geom_node_point(aes(shape = shape, color=name, size=shape)) +
  geom_node_text(aes(label = label, fontface=shape),
                 size=nodes$textsize,
                 fontface = nodes$fontface,
                 hjust = nodes$hjust,
                 color = nodes$textcolor) +
  # scale_y_continuous(expand = expansion(add=c(0.1, 0.1))) +
  scale_x_continuous(expand = expansion(add=c(0.3, 0.4))) +
  scale_edge_width_continuous(range=c(0.5, 2), guide="none") +
  scale_color_manual(values=c("Interoception" = "red",
                              "Muscles" = "#E91E63",
                              "Heart" = "#E91E63",
                              "Breathing" = "#E91E63",
                              "Muscle" = "#E91E63",
                              "Skin" = "#FF5722",
                              "Tickle" = "#FF5722",
                              "Itch" = "#FF5722",
                              "Bruise" = "#FF5722",
                              "Elimination" = "#FFC107",
                              "Defecate" = "#FFC107",
                              "Sneeze" = "#FFC107",
                              "Vomit" = "#FFC107",
                              "Homeostasis" = "#9C27B0",
                              "Hungry" = "#9C27B0",
                              "Thirsty" = "#9C27B0",
                              "Sex_arousal" = "#9C27B0"), guide="none") +
  scale_shape_manual(values=c("oval"="circle", "rect"="square"), guide="none") +
  scale_size_manual(values = c(35, 5), guide="none") +
  labs(title = "C. Structural Equation Model")  +
  theme_graph() +
  theme(plot.title = element_text(hjust=0.5))
p_graph

Omega

semTools::compRelSEM(cfa4h, higher="Interoception")
##   Homeostasis          Skin   Elimination        Muscle Interoception 
##         0.713         0.689         0.670         0.734         0.716
scores <- as.data.frame(predict(cfa4h))
scores$Participant <- df3$Participant
write.csv(scores, "IASR_scores.csv", row.names = FALSE)

Comparison with MAIA

df <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |> 
  merge(read.csv("IASR_scores.csv"), by="Participant") |> 
  select(-ends_with("_RT"), -ends_with("_SD"))

df$Average <- rowMeans(select(df, starts_with("Item_IAS")))

vars <- c("Homeostasis", "Skin", "Elimination", "Muscle", "Interoception")


make_correlation <- function(x, y) {
  set.seed(3)
  cor <- correlation::correlation(
    x, y,
    bayesian = TRUE,
    bayesian_prior = "medium.narrow"
  ) |>
    datawizard::data_remove(c("ROPE_Percentage"))

  cor$BF_Spearman <- format_bf(
    correlation::correlation(
      x, y,
      bayesian = TRUE,
      ranktransform = TRUE,
      bayesian_prior = "medium.narrow"
    )$BF,
    name = NULL, stars = FALSE
  )
  
  r <- filter(cor, BF >= 10)

  if(nrow(r) == 0) return(list(data = cor, print = "No significant (BF >= 10) correlations."))

  
  to_print <- r |>
    arrange(Parameter1, desc(BF)) |>
    format_table() |>
    mutate(Coefficient = paste(rho, `95% CI`)) |>
    select(Parameter1, Parameter2, Coefficient, BF, `BF (Spearman)` = BF_Spearman)
  
  list(data = cor, print = to_print)
}

maia <- make_correlation(df[vars], select(df, starts_with("MAIA_"))) 
export_table(maia$print)
## NULL

maia$data |> 
  arrange(Parameter2) |> 
  tidygraph::as_tbl_graph() |> 
  mutate(Type = ifelse(str_detect(name, "MAIA_"), "MAIA", "IAS")) |> 
  tidygraph::activate("edges") |> 
  filter(BF > 10) |> 
  ggraph(layout = 'stress', weight=rho) + 
  # geom_edge_arc(aes(edge_width = rho),
  #               strength=0) +
  geom_edge_link(aes(edge_alpha = after_stat(index), edge_width = rho), 
                 edge_colour = "black") +
  geom_node_point(aes(color = Type)) +
  geom_node_text(aes(label=name)) +
  scale_edge_width_continuous(range = c(0.1, 3)) +
  coord_fixed() +
  theme_graph()

Study 3

Personality

ipip <- make_correlation(df[vars], select(df, starts_with("IPIP6_"))) 
export_table(ipip$print)
## Parameter1    |              Parameter2 |          Coefficient |       BF | BF (Spearman)
## -----------------------------------------------------------------------------------------
## Homeostasis   |   IPIP6_HonestyHumility | -0.21 [-0.29, -0.13] | 6.82e+03 |        623.03
## Homeostasis   | IPIP6_Conscientiousness |  0.17 [ 0.09,  0.26] |   250.18 |        103.74
## Homeostasis   |       IPIP6_Neuroticism | -0.16 [-0.25, -0.07] |    99.00 |          5.14
## Homeostasis   |     IPIP6_Agreeableness |  0.14 [ 0.05,  0.22] |    17.91 |         39.05
## Interoception |     IPIP6_Agreeableness |  0.15 [ 0.07,  0.23] |    48.73 |         76.28
## Interoception | IPIP6_Conscientiousness |  0.14 [ 0.05,  0.22] |    15.14 |          9.30
## Interoception |   IPIP6_HonestyHumility | -0.13 [-0.22, -0.05] |    12.19 |          4.32
## Interoception |       IPIP6_Neuroticism | -0.13 [-0.21, -0.04] |    10.74 |          1.99
## Muscle        |     IPIP6_Agreeableness |  0.16 [ 0.07,  0.24] |    98.21 |        374.28
## Muscle        |          IPIP6_Openness |  0.14 [ 0.06,  0.22] |    20.00 |         35.95

pid <- make_correlation(df[vars], select(df, starts_with("PID5_"))) 
export_table(pid$print)
## Parameter1  |          Parameter2 |          Coefficient |    BF | BF (Spearman)
## --------------------------------------------------------------------------------
## Homeostasis |     PID5_Detachment | -0.16 [-0.24, -0.07] | 77.71 |         20.56
## Homeostasis |   PID5_Psychoticism | -0.16 [-0.24, -0.07] | 71.17 |         25.13
## Homeostasis | PID5_NegativeAffect | -0.13 [-0.22, -0.05] | 11.83 |          3.83
## Muscle      | PID5_NegativeAffect | -0.14 [-0.22, -0.05] | 13.36 |          3.42

Psychopathology

asq <- make_correlation(df[vars], select(df, starts_with("ASQ_")))
export_table(asq$print)
## Parameter1    |       Parameter2 |       Coefficient |     BF | BF (Spearman)
## -----------------------------------------------------------------------------
## Elimination   |    ASQ_Switching | 0.16 [0.07, 0.25] |  98.16 |        > 1000
## Homeostasis   |    ASQ_Switching | 0.17 [0.08, 0.26] | 193.29 |        312.15
## Homeostasis   | ASQ_SocialSkills | 0.15 [0.06, 0.23] |  27.59 |         28.62
## Interoception |    ASQ_Switching | 0.17 [0.08, 0.25] | 155.71 |        542.73
## Muscle        |    ASQ_Switching | 0.13 [0.05, 0.22] |  10.26 |         47.62
## Skin          |     ASQ_Patterns | 0.18 [0.09, 0.26] | 342.90 |        177.31
## Skin          |    ASQ_Switching | 0.15 [0.06, 0.23] |  28.41 |         70.39
## Skin          | ASQ_SocialSkills | 0.14 [0.05, 0.22] |  17.16 |         11.83
spq <- make_correlation(df[vars], select(df, starts_with("SPQ_"))) 
export_table(spq$print)
## Parameter1    |            Parameter2 |          Coefficient |     BF | BF (Spearman)
## -------------------------------------------------------------------------------------
## Homeostasis   |     SPQ_SocialAnxiety | -0.17 [-0.25, -0.09] | 214.40 |         80.54
## Homeostasis   |    SPQ_NoCloseFriends | -0.17 [-0.25, -0.09] | 171.05 |         69.36
## Homeostasis   |         SPQ_OddSpeech | -0.15 [-0.23, -0.06] |  34.50 |          6.37
## Homeostasis   | SPQ_ConstrictedAffect | -0.14 [-0.23, -0.05] |  19.14 |          8.61
## Interoception |    SPQ_NoCloseFriends | -0.16 [-0.24, -0.07] |  77.13 |         65.39
## Interoception |     SPQ_SocialAnxiety | -0.15 [-0.23, -0.06] |  29.14 |         36.91
## Muscle        |    SPQ_NoCloseFriends | -0.15 [-0.24, -0.06] |  32.95 |         24.67
## Muscle        |     SPQ_SocialAnxiety | -0.14 [-0.23, -0.06] |  14.36 |         17.12
bpd <- make_correlation(df[vars], select(df, starts_with("BPD"))) 
export_table(bpd$print)
## NULL
phq <- make_correlation(df[vars], select(df, starts_with("PHQ4_"))) 
export_table(phq$print)
## Parameter1  |   Parameter2 |          Coefficient |    BF | BF (Spearman)
## -------------------------------------------------------------------------
## Homeostasis | PHQ4_Anxiety | -0.14 [-0.23, -0.05] | 16.03 |          7.73

Other

pi <- make_correlation(df[vars], select(df, starts_with("PI_")))
export_table(pi$print)
## NULL
lie <- make_correlation(df[vars], select(df, starts_with("LIE_")))
export_table(lie$print)
## NULL
gcbs <- make_correlation(df[vars], select(df, starts_with("GCBS_"))) 
export_table(gcbs$print)
## Parameter1    |              Parameter2 |       Coefficient |    BF | BF (Spearman)
## -----------------------------------------------------------------------------------
## Elimination   | GCBS_GlobalConspiracies | 0.14 [0.05, 0.22] | 14.82 |         22.87
## Homeostasis   | GCBS_GlobalConspiracies | 0.15 [0.07, 0.24] | 55.63 |        494.37
## Homeostasis   |  GCBS_PersonalWellbeing | 0.14 [0.06, 0.23] | 23.58 |         51.00
## Interoception | GCBS_GlobalConspiracies | 0.14 [0.06, 0.23] | 27.79 |         51.03
## Muscle        | GCBS_GlobalConspiracies | 0.13 [0.05, 0.22] | 13.74 |         22.97
age <- make_correlation(df[vars], data.frame(Age = df$Age)) 
export_table(age$print)
## NULL
sex <- make_correlation(df[vars], data.frame(Sex = ifelse(df$Sex == "Male", 0, 1)) )
export_table(sex$print)
## Parameter1  | Parameter2 |          Coefficient |    BF | BF (Spearman)
## -----------------------------------------------------------------------
## Homeostasis |        Sex | -0.14 [-0.23, -0.05] | 15.62 |          4.59
p_cor <- rbind(
  mutate(maia$data, Domain = "MAIA"),
  mutate(ipip$data, Domain = "Normal Personality"),
  mutate(pid$data, Domain = "Maladaptive Personality"),
  mutate(asq$data, Domain = "Autistic Traits"),
  mutate(spq$data, Domain = "Schizotypic Traits"),
  mutate(phq$data, Domain = "Mood"),
  mutate(bpd$data, Domain = "Mood"),
  mutate(gcbs$data, Domain = "Conspiracy Beliefs"),
  mutate(pi$data, Domain = "Primal World Beliefs"),
  mutate(lie$data, Domain = "Lying Profile")
  # mutate(age$data, Domain = "Demographic"),
  # mutate(sex$data, Domain = "Demographic")
) |> 
  mutate(alpha = ifelse(BF >= 10, "BF >= 10", "BF < 10"),
         Parameter1 = fct_relevel(Parameter1, "Interoception", "Muscle", "Homeostasis", "Elimination", "Skin"),
         Domain = fct_relevel(Domain, "MAIA", "Normal Personality", "Maladaptive Personality", "Autistic Traits", "Schizotypic Traits", "Mood", "Lying Profile", "Conspiracy Beliefs", "World Beliefs"),
         Parameter2 = str_remove_all(Parameter2, "MAIA_|IPIP6_|PID5_|ASQ_|SPQ_|GCBS_|PHQ4_|LIE_|PI_"),
         Parameter2 = str_replace(Parameter2, "Sex", "Sex (Female vs. Male)"),
         Parameter2 = str_replace(Parameter2, "HonestyHumility", "Honesty-Humility"),
         Parameter2 = str_replace(Parameter2, "AttentionRegulation", "Attention Regulation"),
         Parameter2 = str_replace(Parameter2, "EmotionalAwareness", "Emotional Awareness"),
         Parameter2 = str_replace(Parameter2, "SelfRegulation", "Self-Regulation"),
         Parameter2 = str_replace(Parameter2, "BodyListening", "Body Listening"),
         Parameter2 = str_replace(Parameter2, "NotWorrying", "Not Worrying"),
         Parameter2 = str_replace(Parameter2, "NotDistracting", "Not Distracting"),
         Parameter2 = str_replace(Parameter2, "NegativeAffect", "Negative Affect"),
         Parameter2 = str_replace(Parameter2, "Patterns", "Patterns and Numbers"),
         Parameter2 = str_replace(Parameter2, "Routine", "Routines"),
         Parameter2 = str_replace(Parameter2, "LackSocialSkills", "Lack of Social Skills"),
         Parameter2 = str_replace(Parameter2, "LowAttentionalSwitching", "Low Attentional Switching"),
         Parameter2 = str_replace(Parameter2, "MagicalThinking", "Magical Thinking"),
         Parameter2 = str_replace(Parameter2, "UnusualPerceptions", "Unusual Perceptions"),
         Parameter2 = str_replace(Parameter2, "OddSpeech", "Odd Speech"),
         Parameter2 = str_replace(Parameter2, "ConstrictedAffect", "Constricted Affect"),
         Parameter2 = str_replace(Parameter2, "SocialAnxiety", "Social Anxiety"),
         Parameter2 = str_replace(Parameter2, "NoCloseFriends", "No Close Friends"),
         Parameter2 = str_replace(Parameter2, "GlobalConspiracies", "Global Conspiracies"),
         Parameter2 = str_replace(Parameter2, "PersonalWellbeing", "Personal Wellbeing"),
         Parameter2 = str_replace(Parameter2, "InformationControl", "Information Control"),
         Parameter2 = str_replace(Parameter2, "GovernmentMalfeasance", "Government Malfeasance"),
         Parameter2 = fct_reorder(Parameter2, rho),
         color = as.character(sign(rho))) |>
  filter(Parameter2 != "Total") |> 
  ggplot(aes(x = rho, y = Parameter2)) +
  geom_bar(aes(alpha = alpha, fill = rho), stat = "identity") +
  geom_linerange(aes(alpha = alpha, xmin = CI_low, xmax = CI_high, color = color)) +
  geom_vline(xintercept = 0) +
  facet_grid(Domain~Parameter1, scales="free_y", switch ="y") +
  scale_alpha_manual(values = c(0.3, 1)) +
  scale_fill_gradient2(low = "#2196F3", high = "#FF5722", limits = c(-0.4, 0.4)) +
  scale_colour_manual(values = c("#2196F3", "#FF5722")) +
  scale_x_continuous(expand = c(0, 0), breaks = c(-0.3, 0, 0.3)) +
  coord_cartesian(xlim = c(-0.45, 0.45)) +
  labs(title = "Correlates of IAS-R Dimensions") +
  guides(alpha = "none", fill = "none", color="none") +
  theme_minimal() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        plot.title = element_text(face = "bold", hjust = 0.5),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        strip.placement.y = "outside",
        strip.background.x = element_rect(fill = "#E0E0E0", color = "white"),
        strip.text.x = element_text(face = "bold", size = rel(1.1)),
        strip.text.y = element_text(face = "bold", size = rel(0.95)))
p_cor

Figures

p_dens <- rbind(
  mutate(dens1, Study = "Study 1a"),
  filter(mutate(dens2, Study = "Study 1b"), Parameter %in% dens1$Parameter),
  filter(mutate(dens3, Study = "Study 2"), Parameter %in% dens1$Parameter)
) |> 
  ggplot(aes(x=x, y=y, color=Parameter)) +
  geom_line() +
  facet_wrap(~Study, scales = "free", ncol = 1, strip.position = "right") +
  scale_color_metro_d() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(y = "Density", x = "Answer", title = "A. Distribution of Answers") +
  guides(color = guide_legend(title="Item", ncol=1, byrow=TRUE, keyheight=0)) +
  theme_modern(axis.title.space = 5) +
  theme(axis.text.y = element_blank(),
        axis.line.y = element_blank(),
        axis.title.y = element_blank(),
        plot.title = element_text(face = "bold", hjust = 0),
        legend.text = element_text(size = rel(0.8)),
        legend.spacing.y = unit(0.1, 'cm'),
        strip.background = element_rect(fill = "#E0E0E0", color = "white")) 

p_efa <- rbind(
  attributes(efa4)$dataset[names(data)] |> 
    parameters::factor_analysis(n=4, rotation = "oblimin", sort = TRUE, fm="ml") |> 
    data_plot() |> 
    mutate(Study = "Study 1a"),
  attributes(efa4b)$dataset[names(data)] |> 
    parameters::factor_analysis(n=4, rotation = "oblimin", sort = TRUE, fm="ml") |> 
    data_plot() |> 
    mutate(Study = "Study 1b"),
  efa4c |> 
    data_plot() |> 
    mutate(Study = "Study 2") 
) |> 
  mutate(Variable = fct_relevel(Variable, sort(levels(Variable), decreasing = TRUE)),
         fill = as.factor(sign(y)),
         x = abs(y),
         Component = str_replace_all(Component, "ML", "Factor ")) |> 
  ggplot(aes(x = x, y = Variable, fill = fill)) +
  geom_bar(stat = "identity") + 
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 0.25, 0.5, 0.75), labels = c("0", ".25", ".50", ".75")) +
  scale_fill_manual(values = c("#E91E63", "#4CAF50")) +
  labs(x = "Loading", title = "B. Exploratory Factor Analysis") +
  guides(fill = "none") +
  theme_modern(axis.title.space = 5) +
  theme(panel.grid.major.y = element_line(),
        axis.title.y = element_blank(),
        axis.line.x = element_blank(),
        axis.text.y = element_text(size = rel(0.5)),
        axis.text.x = element_text(size = rel(0.8)),
        plot.title = element_text(face = "bold", hjust = 0),
        strip.placement = "outside",
        strip.background.y = element_rect(fill = "#E0E0E0", color = "white"),
        strip.background.x = element_rect(fill = "#F8BBD0", color = "white")) +
  facet_grid(Study ~ Component, switch = "y")

fig <- ((wrap_elements(p_dens) + wrap_elements(p_efa)) / p_graph)
ggsave("figures/Figure1.png", fig, width=fig.width, height=fig.width*1.3, dpi=200)
ggsave("figures/Figure2.png", p_cor, width=21*0.35, height=29.7*0.44, dpi=200, bg="white")