Mega-analysis of the Interoceptive Accuracy Scale (IAS)

Study 1

Code
library(tidyverse)
library(easystats)
library(patchwork)
library(lavaan)
library(ggraph)
library(tidySEM)
library(EGAnet)
library(haven)

ias_vars <- c(
  "Heart", "Hungry", "Breathing", "Thirsty", "Urinate", "Defecate", "Taste", "Vomit", "Sneeze", "Cough", "Temp",
  "Sex_arousal", "Wind", "Burp", "Muscles", "Bruise", "Pain", "Blood_Sugar", "Affective_touch", "Tickle", "Itch"
)

Data Preparation

  • Murphy (2020)
    • https://osf.io/3m5nh
Code
df1a <- 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"))))

df1b <- 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"))))
  • Gaggero (2021)
    • https://osf.io/5x9sg
    • BPQ_A = autonomic symptoms
Code
load("../data/Gaggero2021/data.rda")
df2 <- data |> 
  mutate(Gender = as.character(gender)) |> 
  mutate(across(starts_with("IAS "), as.numeric)) |>
  rename(
    Age=age, Heart = `IAS 1`, Hungry = `IAS 2`, Breathing = `IAS 3`, Thirsty = `IAS 4`,
    Urinate = `IAS 5`, Defecate = `IAS 6`, Taste = `IAS 7`, Vomit = `IAS 8`, Sneeze = `IAS 9`,
    Cough = `IAS 10`, Temp = `IAS 11`, Sex_arousal = `IAS 12`, Wind = `IAS 13`, Burp = `IAS 14`,
    Muscles = `IAS 15`, Bruise = `IAS 16`, Pain = `IAS 17`, Blood_Sugar = `IAS 18`,
    Affective_touch = `IAS 19`, Tickle = `IAS 20`, Itch = `IAS 21`,
    MAIA_Trusting = MAIA_BodyTrusting)
rm(data)


df2$BPQ_BodyAwareness <- rowMeans(select(df2, `BPQ 1`:`BPQ 26`), na.rm = TRUE)
df2$BPQ_AutonomicReactivity <- rowMeans(select(df2, `BPQ 27`:`BPQ 46`), na.rm = TRUE)
  • Campos (2022) - Study 1
    • https://osf.io/j6ef3
Code
df3 <- haven::read_sav("../data/Campos2022/Dataset_Test.sav") |>
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(Sex == 1, "Male", ifelse(Sex == 0, "Female", "Other")))) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21,
    TAS_EOT = TAS_ExternallyOrientedThinking, TAS_DIF = TAS_IdentifyingFeelings, TAS_DDF = TAS_DescribingFeelings
  )


df3$BPQ_BodyAwareness <- rowMeans(select(df3, BPQ1:BPQ26), na.rm = TRUE)
df3$BPQ_AutonomicReactivity <- rowMeans(select(df3, BPQ27:BPQ46), na.rm = TRUE)
  • Todd (2022)
    • https://osf.io/ms354/
    • MAIA does not include items for not-worrying and not-distracting
    • BPQ 26 - body awareness
Code
df4 <- haven::read_sav("../data/Todd2022/CompleteData_new.sav") |>
  rename_with(\(x) str_remove(x, "IAS_"), .cols = starts_with("IAS_")) |>
  rename(
    Heart = IAS1, Hungry = IAS2, Breathing = IAS3, Thirsty = IAS4,
    Urinate = IAS5, Defecate = IAS6, Taste = IAS7, Vomit = IAS8, Sneeze = IAS9,
    Cough = IAS10, Temp = IAS11, Sex_arousal = IAS12, Wind = IAS13, Burp = IAS14,
    Muscles = IAS15, Bruise = IAS16, Pain = IAS17, Blood_Sugar = IAS18,
    Affective_touch = IAS19, Tickle = IAS20, Itch = IAS21
  )

df4$MAIA_AttentionReg <- rowMeans(select(df4, starts_with("MAIA_AR")), na.rm = TRUE)
df4$MAIA_BodyListening <- rowMeans(select(df4, starts_with("MAIA_BL")), na.rm = TRUE)
df4$MAIA_Trusting <- rowMeans(select(df4, starts_with("MAIA_TR")), na.rm = TRUE)
df4$MAIA_EmoAwareness <- rowMeans(select(df4, starts_with("MAIA_EA")), na.rm = TRUE)
df4$MAIA_SelfReg<- rowMeans(select(df4, starts_with("MAIA_SR")), na.rm = TRUE)
df4$MAIA_Noticing <- rowMeans(select(df4, starts_with("MAIA_NT")), na.rm = TRUE)

df4$BPQ_BodyAwareness <- rowMeans(select(df4, starts_with("BPQ_")), na.rm = TRUE)
  • Arslanova (2022)
    • https://osf.io/mp3cy/
    • note: Final.xlsx includes information about which participants passed the attention checks
Code
df5_attention <- readxl::read_excel("../data/Arslanova2022/Participants2.xlsx", 
    sheet = "FINAL") |>
  filter(ActorHR == 1) |>
  select(Subj.ID, Gender, Age) |>
  mutate(Gender = as.numeric(Gender), Age = as.numeric(Age)) |>
  filter(!is.na(Age)) |> #error on the documentation of this participant - attention check failed but not noted
  mutate(Gender = case_when(
    Gender== 0 ~ "Male",
    Gender== 1 ~ "Female" # based on paper reporting 65 women
  ))
New names:
• `Why No?` -> `Why No?...4`
• `Why No?` -> `Why No?...6`
• `` -> `...12`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`
Code
df5 <- readxl::read_excel("../data/Arslanova2022/Murphy_Iacc_scale.xlsx")  |>
  filter(str_detect(Question.Key, "IAC_")) |> 
  filter(str_detect(Question.Key, "-quantised")) |> 
  pivot_wider(names_from = Question.Key, values_from = Response) |>
  rename_with(\(x) str_remove(x, "IAC_"), .cols = starts_with("IAC_")) |> 
  rename_with(\(x) str_remove(x, "-quantised"), .cols = everything()) |> 
  rename(Sex_arousal = SexuallyAroused, Itch = Itchy, Sex_arousal = SexuallyAroused,
         Temp = HotCold, Tickle = Ticklish, Breathing= BreatheFast, 
         Affective_touch= PleasantAffectionate, Blood_Sugar= BloodSugar, 
         Muscles=TiredMuscles, Heart= FastHeart, Taste=Tastes) |> 
  rename(Subj.ID = "Participant.Public.ID") |>
  select(-starts_with("Participant")) 

df5 <- df5 |>
  filter(Subj.ID %in% df5_attention$Subj.ID) |>
  select(-Subj.ID) |>
  mutate(across(everything(), as.numeric))
  • Brand (2022)
    • https://osf.io/xwz6g/
Code
df6 <- haven::read_sav("../data/Brand2022/LatentVariableApproach.sav") |> 
  select(-SERIAL) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    PHQ4_Depression= PHQ2_sum,
    PHQ4_Anxiety = GAD2_sum
  ) 
  • Brand (2023)
    • https://osf.io/3f2h6/
Code
df7a <- haven::read_sav("../data/Brand2023/Data_Giessen.sav") |> 
  select(-COUNTRY_OTHER) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(SEX == 1, "Male", ifelse(SEX == 2, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    MAIA_Noticing = MAIA_NOTICING, MAIA_NotDistracting = MAIA_NON_DISTRACTING,
    MAIA_NotWorrying = MAIA_NOT_WORRYING, MAIA_AttentionReg = mAIA_ATTENTION_REGULATION,
    MAIA_SelfReg = MAIA_SELF_REGULATION, MAIA_BodyListening = MAIA_BODY_LISTENING,
    MAIA_Trusting = MAIA_TRUSTING, MAIA_EmoAwareness = MAIA_EMOTIONAL_AWARENESS,
    STAIT_Sum = STAI_T_Sumscore
  ) |> 
  mutate(across(starts_with("MAIA_"), \(x) x - 1),
         across(starts_with("NEOFFI_"),  \(x) datawizard::rescale(x, to = c(0, 1)))) # Original scale 1-6

df7a$STAIT_Mean <- df7a$STAIT_Sum/20
df7a$Neuroticism <- rowMeans(select(df7a,NEOFFI_01:NEOFFI_12))


df7b <- haven::read_sav("../data/Brand2023/Data_Mainz.sav") |> 
  select(-COUNTRY_OTHER, -EDUCATION_OTHER, -Sample) |> 
  mutate_all(as.numeric) |> 
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IAS_01, Hungry = IAS_02, Breathing = IAS_03, Thirsty = IAS_04,
    Urinate = IAS_05, Defecate = IAS_06, Taste = IAS_07, Vomit = IAS_08, Sneeze = IAS_09,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    MAIA_Noticing = maiaNotice, MAIA_NotDistracting = maiaNonDistracting,
    MAIA_NotWorrying = maiaNotWorrying, MAIA_AttentionReg = maiaAttentionRegulation,
    MAIA_SelfReg = maiaSelfRegulation, MAIA_BodyListening = maiaBodyListening,
    MAIA_Trusting = maiaTrusting, MAIA_EmoAwareness = maiaEmotionalAwareness) |> 
  mutate(across(starts_with("MAIA_"), \(x) datawizard::rescale(x, to = c(0, 5), range = c(1, 5))),# Original scale 1-5
          across(starts_with("PHQ15_"), \(x) datawizard::rescale(x, to = c(0,2), range = c(1,3))),
         
           across(PHQ9_01:PHQ9_09, \(x) datawizard::rescale(x, to = c(0,3), range =c(1,4))))

df7b$BPQ_BodyAwareness <- rowMeans(select(df7b, BPQ_01:BPQ_12), na.rm = TRUE)
df7b$PHQ15_Sum <-rowMeans(select(df7b,PHQ15_01:PHQ15_15recode), na.rm = TRUE)
df7b$PHQ9_Sum <- rowMeans(select(df7b, PHQ9_01:PHQ9_09), na.rm = TRUE)
df7b$STAIT_Mean <- df7b$STAIT_Sum/20
  
  
df7c <- haven::read_sav("../data/Brand2023/Data_PotVie.sav") |> 
  select(-ID) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(GENDER == 2, "Male", ifelse(GENDER == 1, "Female", "Other")))) |> 
  rename(
    Age=AGE, Heart = IA02_01, Hungry = IA02_02, Breathing = IA02_03, Thirsty = IA02_04,
    Urinate = IA02_05, Defecate = IA02_06, Taste = IA02_07, Vomit = IA02_08, Sneeze = IA02_09,
    Cough = IA02_10, Temp = IA02_11, Sex_arousal = IA02_12, Wind = IA02_13, Burp = IA02_14,
    Muscles = IA02_15, Bruise = IA02_16, Pain = IA02_17, Blood_Sugar = IA02_18,
    Affective_touch = IA02_19, Tickle = IA02_20, Itch = IA02_21,
    MAIA_Noticing = MAIA_noticing, MAIA_NotDistracting = MAIA_notdis, 
    MAIA_NotWorrying = MAIA_notwor, MAIA_AttentionReg = MAIA_attreg,
    MAIA_SelfReg = MAIA_selfreg, MAIA_BodyListening = MAIA_bodlis,
    MAIA_Trusting = MAIA_trust, MAIA_EmoAwareness = MAIA_emoaw,
    ICQ = ICQ_sum,
    STAIT_Sum = STAIT_sum,
    TAS_DIF = TAS_ident, TAS_EOT = TAS_ext, TAS_DDF = TAS_des) |> 
  mutate(across(starts_with("MAIA_"), \(x) datawizard::rescale(x, to = c(0, 5), range = c(0, 6)))) # Original scale 0-6

df7c$BPQ_BodyAwareness <- df7c$BPQ_sum / 26
df7c$STAIT_Mean <- df7c$STAIT_Sum/20
  • Lin (2023)
    • https://osf.io/3eztd
    • Note: No tickle because same Chinese character
    • missing phq9_9 item
Code
df8a <- haven::read_sav("../data/Lin2023/Study 1 & 3.sav") |>
  select(-sex) |> 
  mutate_all(as.numeric) |>
  mutate(Gender = as.character(ifelse(sex_dummy == 1, "Male", ifelse(sex_dummy == 0, "Female", "Other")))) |>
  rename(
    Age = age,
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH
  ) 

df8a$BPQ_BodyAwareness <- rowMeans(select(df8a,BPQ1:BPQ12), na.rm = FALSE)


df8b <- haven::read_sav("../data/Lin2023/Study 2.sav") |>
  mutate(Gender = as.character(ifelse(Sex == "男", "Male", ifelse(Sex == "女", "Female", "Other")))) |>
  rename(
    Heart = Heart, Hungry = HUNGR, Breathing = BREAT, Thirsty = Thirs,
    Urinate = URINA, Defecate = Defec, Taste = TASTE, Vomit = VOMIT, Sneeze = Sneez,
    Cough = COUGH, Temp = TEMPE, Sex_arousal = SEXAR, Wind = WIND, Burp = Burp,
    Muscles = MUSCL, Bruise = Bruis, Pain = PAIN, Blood_Sugar = BloSu,
    Affective_touch = Touch, Itch = ITCH,
    PHQ15_Sum = PHQ15_sum, PHQ9_sum = PHQ9_sum)

df8b$PHQ9_Sum <- df8b$PHQ9_sum / 9
df8b$PHQ15_Sum <- df8b$PHQ15_Sum/ 15
  • VonMohr (2023) - Study 3
    • https://osf.io/7p9u5/
Code
df9 <- read.csv("../data/VonMohr2023/DataSet_study3.csv")
df9 <- df9[complete.cases(select(df9, starts_with("IAS_"))), ]
df9 <- df9 |>
  mutate(Gender = as.character(ifelse(GENDER == "Man", "Male", ifelse(GENDER == "Woman", "Female", "Other")))) |>
  rename(
    Age=AGE, Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |>
  mutate(across(starts_with("BFI_"), \(x) datawizard::rescale(x, to = c(0, 1)))) |>
  mutate(
    Neuroticism = (BFI_5 + BFI_10 + BFI_15)/3,
    Agreeableness = (BFI_3 + BFI_6 + BFI_13)/3,
    Extraversion = (BFI_2 + BFI_8 + BFI_12)/3,
    Conscientiousness = (BFI_1 + BFI_7 + BFI_11)/3,
    Openness = (BFI_4 + BFI_9 + BFI_14)/3)
  • Makowski
    • Note: sample has some missing items (No Taste, Cough, Blood_Sugar)
    • Note: Analog scales
Code
df10 <- read.csv("https://raw.githubusercontent.com/RealityBending/IllusionGameReliability/main/data/preprocessed_questionnaires.csv") |>
  rename(
    Gender = Sex, Heart = Item_IAS_Interoception_1, Hungry = Item_IAS_Interoception_2,
    Breathing = Item_IAS_Interoception_3, Thirsty = Item_IAS_Interoception_4,
    Temp = Item_IAS_Interoception_5, Sex_arousal = Item_IAS_Interoception_6,
    Urinate = Item_IAS_Elimination_1, Defecate = Item_IAS_Elimination_2,
    Vomit = Item_IAS_Elimination_3, Wind = Item_IAS_Expulsion_1,
    Burp = Item_IAS_Expulsion_2, Sneeze = Item_IAS_Expulsion_3,
    Muscles = Item_IAS_Nociception_1, Bruise = Item_IAS_Nociception_2,
    Pain = Item_IAS_Nociception_3, Affective_touch = Item_IAS_Skin_1,
    Tickle = Item_IAS_Skin_2, Itch = Item_IAS_Skin_3,
    MAIA_EmoAwareness = MAIA_EmotionalAwareness, MAIA_AttentionReg = MAIA_AttentionRegulation,
    MAIA_SelfReg = MAIA_SelfRegulation,
  ) |>
  filter(!is.na(Urinate)) |> 
  datawizard::rescale(select=ias_vars, to = c(1, 5), range = c(0, 1)) |> 
  mutate(across(starts_with("MAIA_"), \(x) datawizard::rescale(x, to = c(0, 5))),
         across(starts_with("PI_"), \(x) datawizard::rescale(x, to = c(0, 1))),
         across(starts_with("Item_PHQ4_Depression_"), \(x) datawizard::rescale(x, to = c(1,4))),
        across(starts_with("Item_PHQ4_Anxiety_"), \(x) datawizard::rescale(x, to = c(1,4))),
        across(starts_with("IPIP6_"), \(x) datawizard::rescale(x, to = c(0,1))),
        PHQ4_Depression = Item_PHQ4_Depression_3 + Item_PHQ4_Depression_4,
          PHQ4_Anxiety = Item_PHQ4_Anxiety_1 + Item_PHQ4_Anxiety_2) |>
  rename(
     Agreeableness = IPIP6_Agreeableness,
    Conscientiousness = IPIP6_Conscientiousness,
    Extraversion = IPIP6_Extraversion,
    HonestyHumility = IPIP6_HonestyHumility,
    Neuroticism = IPIP6_Neuroticism,
    Openness = IPIP6_Openness
  )
Code
df11  <- read.csv("https://raw.githubusercontent.com/DominiqueMakowski/PHQ4R/main/study2/data/data.csv") |>
  rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    typeSample = Sample,
    BDI_Sum = BDI2_Total) |>
  mutate(
    MAIA2_NotDistracting_1_R = 1 - MAIA2_NotDistracting_1_R,
    MAIA2_NotDistracting_2_R = 1 - MAIA2_NotDistracting_2_R,
    MAIA2_NotDistracting_3_R = 1 - MAIA2_NotDistracting_3_R,
    MAIA2_NotDistracting_4_R = 1 - MAIA2_NotDistracting_4_R,
    MAIA2_NotDistracting_5_R = 1 - MAIA2_NotDistracting_5_R,
    MAIA2_NotDistracting_6_R = 1 - MAIA2_NotDistracting_6_R,
    MAIA2_NotWorrying_1_R = 1 - MAIA2_NotWorrying_1_R,
    MAIA2_NotWorrying_2_R = 1 - MAIA2_NotWorrying_2_R,
    MAIA2_NotWorrying_5_R = 1 - MAIA2_NotWorrying_5_R) |> 
  mutate(across(starts_with("MAIA2_"), \(x) datawizard::rescale(x, to = c(0, 5))),
         across(starts_with("PHQ4_Depression_"), \(x) datawizard::rescale(x, to = c(1,4))),
         across(starts_with("PHQ4_Anxiety_"), \(x) datawizard::rescale(x, to = c(1,4))),
         across(starts_with("STAI5_"), \(x) datawizard::rescale(x, to = c(1, 4)))) |>
   mutate(PHQ4_Depression = PHQ4_Depression_3 + PHQ4_Depression_4,
          PHQ4_Anxiety = PHQ4_Anxiety_1 + PHQ4_Anxiety_2,
          STAIT_Mean = (STAI5_1 + STAI5_2 + STAI5_3 + STAI5_4 + STAI5_4 + STAI5_5)/5) |>
  datawizard::rescale(select=ias_vars, to = c(1, 5), range = c(0, 1))

# Compute subscales

df11$MAIA_AttentionReg <- rowMeans(select(df11 , starts_with("MAIA2_Attention")))
df11$MAIA_EmoAwareness <- rowMeans(select(df11 , starts_with("MAIA2_Emotional")))
df11$MAIA_Noticing <- rowMeans(select(df11 , starts_with("MAIA2_Noticing")))
df11$MAIA_BodyListening <- rowMeans(select(df11 , starts_with("MAIA2_BodyListening")))
df11$MAIA_NotDistracting<- rowMeans(select(df11 , starts_with("MAIA2_NotDistracting")))
df11$MAIA_NotWorrying<- rowMeans(select(df11 , starts_with("MAIA2_NotWorrying")))
df11$MAIA_SelfReg<- rowMeans(select(df11 , starts_with("MAIA2_SelfRegulation")))
df11$MAIA_Trusting<- rowMeans(select(df11 , starts_with("MAIA2_Trusting")))
Code
df12  <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionPrimals/refs/heads/main/data/data_participants.csv") |>
  select("Participant" = "participant_id", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21,
    MAIA_EmoAwareness = MAIA_EmotionalAwareness, 
    MAIA_AttentionReg = MAIA_AttentionRegulation,
    MAIA_SelfReg = MAIA_SelfRegulation
  ) |>
  filter(!is.na(Heart)) |>  # participant is outlier and had total IAS scores below 0.4 
  datawizard::rescale(select=ias_vars, to = c(1, 5), range = c(0, 1)) |> 
  mutate(across(starts_with("MAIA_"), \(x) datawizard::rescale(x, to = c(0, 5)))) 
Code
df13  <- read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/rawdata_participants.csv") |>
  select("Participant", "Gender", "Age", "Ethnicity", starts_with(c("IAS","MAIA","PI", "PHQ4"))) |>
   rename(
    Heart = IAS_1, Hungry = IAS_2, Breathing = IAS_3, Thirsty = IAS_4,
    Urinate = IAS_5, Defecate = IAS_6, Taste = IAS_7, Vomit = IAS_8, Sneeze = IAS_9,
    Cough = IAS_10, Temp = IAS_11, Sex_arousal = IAS_12, Wind = IAS_13, Burp = IAS_14,
    Muscles = IAS_15, Bruise = IAS_16, Pain = IAS_17, Blood_Sugar = IAS_18,
    Affective_touch = IAS_19, Tickle = IAS_20, Itch = IAS_21
  ) |> 
  filter(Participant %in% read.csv("https://raw.githubusercontent.com/RealityBending/InteroceptionScale/refs/heads/main/study2/data/data_participants.csv")$Participant) |>
  mutate(
    MAIA_NotDistracting_1_R = 1 - MAIA_NotDistracting_1_R,
    MAIA_NotDistracting_2_R = 1 - MAIA_NotDistracting_2_R,
    MAIA_NotDistracting_3_R = 1 - MAIA_NotDistracting_3_R,
    MAIA_NotDistracting_4_R = 1 - MAIA_NotDistracting_4_R,
    MAIA_NotDistracting_5_R = 1 - MAIA_NotDistracting_5_R,
    MAIA_NotDistracting_6_R = 1 - MAIA_NotDistracting_6_R,
    MAIA_NotWorrying_1_R = 1 - MAIA_NotWorrying_1_R,
    MAIA_NotWorrying_2_R = 1 - MAIA_NotWorrying_2_R,
    MAIA_NotWorrying_5_R = 1 - MAIA_NotWorrying_5_R) |>
    mutate(across(starts_with("MAIA_"), \(x) datawizard::rescale(x, to = c(0, 5))),
           across(starts_with("PI18_"), \(x) datawizard::rescale(x, to = c(0,1)))) |>
  mutate(
    PI_Alive = (PI18_GA_2 + PI18_GA_14 + (1-PI18_A_12_R) + PI18_A_16 + PI18_A_4)/5, 
    PI_Good = ((1-PI18_GE_5_R) + (1-PI18_GS_13_R) + PI18_GA_14 + PI18_GE_8 + PI18_GA_2 + PI18_GE_1 + (1-PI18_GE_9_R) + (1-PI18_GS_10_R) + (1-PI18_GE_6_R) + PI18_GS_15 + PI18_GS_17 + (1-PI18_GS_11_R) + PI18_GS_3 + PI18_GE_18 +  PI18_GE_7)/ 15,
    PI_Safe = ((1-PI18_GS_13_R) + (1-PI18_GS_10_R) + PI18_GS_15 + PI18_GS_17 + (1-PI18_GS_11_R) + PI18_GS_3)/ 6,
    PI_Enticing = ((1-PI18_GE_5_R) + PI18_GE_8 + PI18_GE_1 +(1-PI18_GE_9_R) + (1-PI18_GE_6_R) + PI18_GE_18 + PI18_GE_7)/7)

  

# Compute subscales
df13$MAIA_AttentionReg <- rowMeans(select(df13 , starts_with("MAIA_Attention")))
df13$MAIA_EmoAwareness <- rowMeans(select(df13 , starts_with("MAIA_Emotional")))
df13$MAIA_Noticing <- rowMeans(select(df13 , starts_with("MAIA_Noticing")))
df13$MAIA_BodyListening <- rowMeans(select(df13 , starts_with("MAIA_BodyListening")))
df13$MAIA_NotDistracting<- rowMeans(select(df13 , starts_with("MAIA_NotDistracting")))
df13$MAIA_NotWorrying<- rowMeans(select(df13 , starts_with("MAIA_NotWorrying")))
df13$MAIA_SelfReg<- rowMeans(select(df13 , starts_with("MAIA_SelfRegulation")))
df13$MAIA_Trusting<- rowMeans(select(df13 , starts_with("MAIA_Trusting")))



df13$PHQ4_Depression <- (df13$PHQ4_Depression_4 + df13$PHQ4_Depression_3)/2
df13$PHQ4_Anxiety <- (df13$PHQ4_Anxiety_2 + df13$PHQ4_Anxiety_1)/2
Code
# pi_vars <- c("PI_Enticing", "PI_Alive", "PI_Safe", "PI_Good", "PI_Changing",
#              "PI_Hierarchical", "PI_Understandable", "PI_AboutMe", "PI_Abundant",
#              "PI_Acceptable", "PI_Beautiful", "PI_Cooperative", "PI_Funny",
#              "PI_Harmless", "PI_Improvable", "PI_Intentional", "PI_Interconnected",
#              "PI_Interesting", "PI_Just", "PI_Meaningful", "PI_NeedsMe",
#              "PI_Pleasurable", "PI_Progressing", "PI_Regenerative", "PI_Stable",
#              "PI_WorthExploring")
# 
# 
# correlation::correlation(
#   df12 [, "IAS_Total", drop = FALSE],
#   data2 = select(df12 , all_of(pi_vars)), p_adjust = "none"
# ) |>
#   # Formatting correlation results
#   mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE))) |>
# 
#   # Releveling factors for visualization
#   mutate(Parameter2 = fct_relevel(Parameter2, rev(pi_vars)),
#          Parameter1 = fct_relevel(Parameter1, "IAS_Total")) |>
# 
#   # Plotting the correlation matrix
#   ggplot(aes(x = Parameter1, y = Parameter2)) +
#   geom_tile(aes(fill = r), color = "white") +
#   geom_text(aes(label = val), size = 3) +
#   labs(title = "Correlation Matrix") +
#   scale_fill_gradient2(low = "blue", high = "red", mid = "white", limit = c(-1, 1), midpoint = 0, guide = guide_colourbar(ticks = FALSE)) +
#   theme_minimal() +
#   theme(
#     legend.title = element_blank(),
#     axis.title.x = element_blank(),
#     axis.title.y = element_blank(),
#     axis.text.x = element_text(hjust = 1)
#   )
  • Poerio et al.,
    • Study 1: https://osf.io/49wbv
Code
df14  <- haven::read_sav("../data/Giulia/Interoceptive Attention ESM.sav") |> 
  select(AGE, GENDER, contains("IAS_ACC"), -IAS_ACC_20, starts_with(c("MAIA", "TAS20", "DEP", "ANX")),  -ANX_21, -TAS20_21) |> #anx_21 , tas20_21 are attention checks
  mutate_all(as.numeric) |>
  rename(
    Age = AGE, Gender = GENDER,
    Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4,
    Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7, Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9,
    Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
    Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18,
    Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_21, Itch = IAS_ACC_22,
    MAIA_EmoAwareness =MAIA_Emotional, MAIA_SelfReg = MAIA_SelfRegulation, MAIA_AttentionReg = MAIA_Attention,
    STAIT_Sum = ANX
  ) |>
   mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender %in% c(3, 5) ~ "Other",
    Gender == 4 ~ "NA"),
    TAS_DIF = TAS20_1 + TAS20_3 + TAS20_6 + TAS20_7 + TAS20_9 + TAS20_13 + TAS20_14,
    TAS_DDF = TAS20_2 + (6-TAS20_4) + TAS20_11 + TAS20_12 + TAS20_17,
    TAS_EOT = (6-TAS20_5) + TAS20_8 + (6-TAS20_10) + TAS20_15 + TAS20_16 + (6-TAS20_18) + (6-TAS20_19) + TAS20_20) 

df14$STAIT_Mean <- df14$STAIT_Sum/20
  • Poerio: unpublished
Code
df15   <- haven::read_sav("../data/Giulia/Sleep and Exteroceptive Interoceptive Sensitivity.sav") |> # 165 participants
  filter(Total_Attention_Fail == 0) |> # 32 participants removed based on failing one or more checks
  select(AGE, GENDER, contains("IAS_ACC"), ) |> 
  mutate_all(as.numeric) |> 
  rename( Age = AGE, Gender = GENDER, 
          Heart = IAS_ACC_1, Hungry = IAS_ACC_2, Breathing = IAS_ACC_3, Thirsty = IAS_ACC_4, Urinate = IAS_ACC_5, Defecate = IAS_ACC_6, Taste = IAS_ACC_7,
          Vomit = IAS_ACC_8, Sneeze = IAS_ACC_9, Cough = IAS_ACC_10, Temp = IAS_ACC_11, Sex_arousal = IAS_ACC_12, Wind = IAS_ACC_13, Burp = IAS_ACC_14,
          Muscles = IAS_ACC_15, Bruise = IAS_ACC_16, Pain = IAS_ACC_17, Blood_Sugar = IAS_ACC_18, Affective_touch = IAS_ACC_19, Tickle =IAS_ACC_20, 
          Itch = IAS_ACC_21 )  |> 
  mutate(Gender = case_when(
    Gender == 1 ~ "Male", 
    Gender == 2 ~"Female", 
    Gender %in% c(3, 5) ~ "Other", 
    Gender == 4 ~ "NA", )) |>
  na.omit()
  • Raquel Arjona and Robyn Scharte data upublished
Code
df16 <- read.csv("../data/Scharte2025/data_IAS_rs.csv") |>
  rename(Age= age, Gender = gender_identity, Education = education,
         Heart = ias1, Hungry = ias2, Breathing = ias3, Thirsty = ias4, Urinate = ias5, 
         Defecate=ias6, Taste = ias7,Vomit = ias8, Sneeze = ias9, Cough = ias10, 
         Temp = ias11, Sex_arousal = ias12, Wind = ias13, Burp = ias14,Muscles = ias15, 
         Bruise = ias16, Pain = ias17, Blood_Sugar = ias18, Affective_touch = ias19, 
         Tickle =ias20, Itch = ias21,
         ) |>
  mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender %in% c(3, 4, 5) ~ "Other",
    Gender == 6 ~ "NA"),
    Education = case_when(
      Education == 6 ~ "Bachelor",
      Education == 7 ~ "Master",
      Education == 8 ~ "Doctorate",
      Education == 2 | Education == 3 ~ "High School", # GCSE + A-levels (16-18 years old)
      Education == 4 | Education == 5 ~ "Other",       # NVQ and foundation year
      Education == 1 ~ "No education"),
    MAIA_Noticing = (maia1 + maia2 + maia3 + maia4)/4,
    MAIA_NotDistracting = (maia5_rev + maia6_rev + maia7_rev +maia8_rev + maia9_rev + maia10_rev)/6, 
    MAIA_NotWorrying = (maia11_rev +maia12_rev  + maia13)/3,
    MAIA_AttentionReg = (maia16 + maia17 + maia18 + maia19 + maia20 + maia21 + maia22)/7,
    MAIA_EmoAwareness = (maia23 + maia24 + maia25 + maia26 +maia27)/5,
    MAIA_SelfReg = (maia28 + maia29 + maia30 + maia31)/4,
    MAIA_BodyListening = (maia32 + maia33 +maia34)/3,
    MAIA_Trusting = (maia35 + maia36 + maia37)/3 ,
    GAD7_Sum = rowSums(across(starts_with("gad1_")), na.rm = TRUE),
    PHQ9_Sum = rowMeans(across(starts_with("phq1_")), na.rm = TRUE)) |>
  filter (!Age == 2004.00, !Age == 1999.00, # misreported age
          att == 2) # only those that select 'no' the attention check
  • Petzke et al 2024
  • STAI - is the stai sf with 6 items for state anxiety (we use stait-t)
  • https://osf.io/ntpu4
Code
df17 <- haven::read_sav("../data/Petzke2024/ETUDE1_main study.sav") |>
  mutate_all(as.numeric) |>
  select(age_to_be_completed_that_year, Gender, starts_with(c("IA0", "PH01")), PHQ_15_total, DIF_total, DDF_total, EOT_total, PHQ_4_total, BPQ_total, starts_with("BP01_")) |>
  rename(Age= age_to_be_completed_that_year,
         Heart = IA01_01, Hungry = IA01_02, Breathing = IA01_03, Thirsty = IA01_04, Urinate = IA01_05, 
         Defecate=IA01_06, Taste = IA01_07,Vomit = IA01_08, Sneeze = IA01_09, Cough = IA01_10, 
         Temp = IA01_11, Sex_arousal = IA01_12, Wind = IA01_13, Burp = IA01_14,Muscles = IA01_15, 
         Bruise = IA01_16, Pain = IA01_17, Blood_Sugar = IA01_18, Affective_touch = IA01_19, 
         Tickle =IA01_20, Itch = IA01_21,
         TAS_DIF = DIF_total, TAS_DDF =  DDF_total, TAS_EOT = EOT_total,
         PHQ15_Sum = PHQ_15_total,
         PHQ4_Sum = PHQ_4_total,
         BPQ_BodyAwareness = BPQ_total,
          ) |>
  mutate(Gender = case_when(
    Gender == 1 ~ "Male",
    Gender == 2 ~ "Female",
    Gender == 3 ~ "Other"))|>
   mutate(
    Age_code = zap_labels(Age),      # remove SPSS labels, get raw numeric codes
    Age = case_when(
      Age_code == -9 ~ NA_real_,
      TRUE ~ 2024 - (2000 + (21 - Age_code))  # compute actual age
    ))
Warning: There were 6 warnings in `mutate()`.
The first warning was:
ℹ In argument: `id = .Primitive("as.double")(id)`.
Caused by warning:
! NAs introduced by coercion
ℹ Run `dplyr::last_dplyr_warnings()` to see the 5 remaining warnings.
Code
# df17$BPQ_BodyAwarenessManual <- rowMeans(df17[, grep("^BP01_", names(df17))], na.rm = TRUE)

#Ranges
df17$PHQ15_Sum <- df17$PHQ15_Sum/15
df17$BPQ_BodyAwareness <- df17$BPQ_BodyAwareness/12

df17$PHQ4_Depression <- df17$PH01_01 + df17$PH01_02
df17$PHQ4_Anxiety <- df17$PH01_03 + df17$PH01_04

Participants

Code
# Save all
data <- list(df1a=df1a, df1b=df1b, df2=df2, df3=df3, df4=df4, df5=df5, df6=df6, df7a=df7a, df7b=df7b, df7c=df7c, df8a=df8a, df8b=df8b, df9=df9, df10=df10, df11 =df11 , df12 =df12 , df13  = df13 , df14 =df14 , df15  =df15, df16 = df16, df17 = df17)
save(data, file = "../data/data.RData")
  • Sample 1a: Data from Murphy’s (2020) study 1, 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).
  • Sample 1b: Data from Murphy’s (2020) study 6, 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).
  • Sample 2: Data from Gaggero’s(2020) study, downloaded from OSF, included 814 participants (Mean age = 24.9, SD = 5.3, range: [18, 58], 0.2% missing; Gender: 60.3% women, 39.4% men, 0.25% non-binary).
  • Sample 3: Data from Campos’s(2022) study, downloaded from OSF, included 515 participants (Mean age = 30.7, SD = 10.5, range: [18, 72]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 59.6% women, 40.4% men, 0.00% non-binary).
  • Sample 4: Data from Todd’s(2022) study, downloaded from OSF, included 802 participants ().
  • Sample 5: Data from Arslanova (2022) study, downloaded from OSF, included 143 participants (Mean age = 28.5, SD = 7.6, range: [18, 73]; Gender: 45.5% women, 54.5% men, 0.00% non-binary).
  • Sample 6: Data from Brand’s(2022) study, downloaded from OSF, included 619 participants (Mean age = 43.9, SD = 14.5, range: [18, 78]; Gender: 78.7% women, 20.2% men, 1.13% non-binary).
  • Sample 7a: Data from Brand’s(2023) study, downloaded from OSF, included 522 participants (Mean age = 23.4, SD = 6.7, range: [18, 79]; Gender: 79.5% women, 19.7% men, 0.77% non-binary).
  • Sample 7b: Data from Brand’s(2023) study, downloaded from OSF, included 1993 participants (Mean age = 32.0, SD = 12.6, range: [16, 81]; Gender: 77.7% women, 21.7% men, 0.60% non-binary).
  • Sample 7c: Data from Brand’s(2023) study, downloaded from OSF, included 808 participants (Mean age = 27.3, SD = 9.3, range: [18, 72], 0.5% missing; Gender: 68.9% women, 30.2% men, 0.87% non-binary).
  • Sample 8a: Data from Lin’s(2023) study, downloaded from OSF, included 1166 participants (Mean age = 32.5, SD = 8.4, range: [16, 60]; Gender: 57.0% women, 43.0% men, 0.00% non-binary).
  • Sample 8b: Data from Lin’s(2023) study, downloaded from OSF, included 500 participants (Mean age = 37.4, SD = 7.4, range: [20, 60]; Sex: 0.0% females, 0.0% males, 100.0% other; Gender: 56.2% women, 43.8% men, 0.00% non-binary).
  • Sample 9: Data from VonMohr’s(2023) study 3, downloaded from OSF, included 21843 participants (Mean age = 56.5, SD = 14.4, range: [18, 93], 0.2% missing; Gender: 73.2% women, 25.1% men, 1.55% non-binary, 0.15% missing).
  • Sample 10: Data from [Makowski’s(2023)] study , downloaded from OSF, included 485 participants (Mean age = 30.1, SD = 10.1, range: [18, 73]; Gender: 50.3% women, 49.7% men, 0.00% non-binary; Education: Bachelor, 45.15%; Doctorate, 1.86%; High school, 34.43%; Master, 15.88%; Other, 2.47%; Prefer not to say, 0.21%).
  • Sample 11: Data from [Makowski’s(2023)] study , downloaded from OSF, included 836 participants (Mean age = 25.1, SD = 11.3, range: [18, 76], 0.1% missing; Gender: 73.8% women, 22.6% men, 2.87% non-binary, 0.72% missing; Education: Bachelor, 22.85%; Doctorate, 2.15%; High School, 63.52%; Master, 6.22%; missing, 0.24%; Other, 5.02%).
  • Sample 12: Data from [Makowski’s(2023)] study , downloaded from OSF, included 146 participants (Mean age = 21.1, SD = 4.3, range: [18, 50], 0.7% missing; Gender: 80.8% women, 15.1% men, 2.74% non-binary, 1.37% missing).
  • Sample 13: Data from [Makowski’s(2023)] study , downloaded from OSF, included 737 participants (Mean age = 36.8, SD = 14.9, range: [17, 87]; Gender: 57.3% women, 41.1% men, 1.63% non-binary).
  • Sample 14: Data from Poerio’s(2024) study , included 107 participants (Mean age = 26.8, SD = 9.2, range: [18, 57]; Gender: 74.8% women, 23.4% men, 1.87% non-binary)
  • Sample 15 : Data from [Poerio’s] study , included 131 participants (Mean age = 30.9, SD = 12.0, range: [18, 60]; Gender: 76.3% women, 22.9% men, 0.76% non-binary)
  • Sample 16 : Data from [Arjona’s] study , included 279 participants (Mean age = 26.4, SD = 13.2, range: [18, 79]; Gender: 67.7% women, 24.7% men, 6.09% non-binary, 1.43% missing; Education: Bachelor, 18.28%; Doctorate, 1.08%; High School, 62.72%; Master, 10.39%; missing, 0.36%; No education, 0.72%; Other, 6.45%)
  • Sample 17 : Data from Petzke’s (2024) study, downloaded from OSF, included 254 participants (Mean age = 31.5, SD = 10.7, range: [22, 69]; Gender: 68.5% women, 30.7% men, 0.79% non-binary).

Total N = 33526 .

Code
desc <- function(df){
  summary <- df |> summarise(
      n = n(),
      mean = mean(Age, na.rm = TRUE),
      sd = sd(Age, na.rm = TRUE),
      min = min(Age, na.rm = TRUE),
      max = max(Age, na.rm = TRUE),
      female = sum(Gender == "Female", na.rm = TRUE)
    )
    return(summary)
}

desc_df1a <- desc(df1a)
desc_df1b <- desc(df1b)
desc_df2 <- desc(df2)
desc_df3 <- desc(df3)
#desc_df4 <- desc(df4) # no demographic data available
desc_df5 <- desc(df5_attention)
desc_df6 <- desc(df6)
desc_df7a <- desc(df7a)
desc_df7b <- desc(df7b)
desc_df7c <- desc(df7c)
desc_df8a <- desc(df8a)
desc_df8b <- desc(df8b)
desc_df9 <- desc(df9)
desc_df10 <- desc(df10)
desc_df11 <- desc(df11)
desc_df12 <- desc(df12)
desc_df13 <- desc(df13)
desc_df14 <- desc(df14)
desc_df15 <- desc(df15)
desc_df16 <- desc(df16)
desc_df17 <- desc(df17)


desc_total <- rbind(desc_df1a, desc_df1b, desc_df2, desc_df3, desc_df5, desc_df6, desc_df7a, desc_df7b, desc_df7c, desc_df8a, desc_df8b, desc_df9, desc_df10, desc_df11 , desc_df12 , desc_df13 , desc_df14 , desc_df15, desc_df16, desc_df17)

## calculate weighted mean

total_mean <- sum(desc_total$mean * desc_total$n) / sum(desc_total$n)
total_sd <- sum(desc_total$sd * desc_total$n) / sum(desc_total$n)
perce_female <- sum(desc_total$female)/sum(desc_total$n) * 100
Code
library(gt)

# APA style ####
gt_apastyle <- function(gt_table, font.size=12) {
  gt_table  |> 
    gt::opt_table_lines(extent = "none") %>%
    gt::tab_options(
      heading.border.bottom.width = 2,
      heading.border.bottom.color = "black",
      heading.border.bottom.style = "solid",
      table.border.top.color = "black",
      table.border.top.style = "solid",
      table.border.top.width = 2,  
      table_body.hlines.color = "white",
      table_body.border.top.color = "black",
      table_body.border.top.style = "solid",
      table_body.border.top.width = 2,
      heading.title.font.size = font.size,
      table.font.size = font.size,
      heading.subtitle.font.size = font.size,
      table_body.border.bottom.color = "black",
      table_body.border.bottom.width = 2,
      table_body.border.bottom.style = "solid",
      column_labels.border.bottom.color = "black",
      column_labels.border.bottom.style = "solid",
      column_labels.border.bottom.width = 1,
      latex.use_longtable = FALSE
    ) |> 
      gt::opt_table_font(font = "times")
}

make_age <- function(age) {
  age <- as.numeric(age) 
  
  mean(age, na.rm = TRUE) |> 
    insight::format_value(digits=1) |> 
    paste0(" ± ", 
           insight::format_value(sd(age, na.rm = TRUE), digits=1)) 
} 

# Create the dataframe
table <- data.frame(
  Sample = c('', 'Sample 1a', 'Sample 1b', 'Sample 2', 'Sample 3', 'Sample 4', 'Sample 5', 'Sample 6', '', 'Sample 7a', 'Sample 7b', 'Sample 7c', '', 'Sample 8a', 'Sample 8b', 'Sample 9', 'Sample 10', 'Sample 11', 'Sample 12', 'Sample 13', 'Sample 14', 'Sample 15', 'Sample 16', 'Sample 17', 'Sample 18'),   
  Reference = c('Murphy et al., (2020)', '', '', 'Gaggero et al., (2021)', 'Campos et al., (2022)', 'Todd et al., (2022)', 'Arslanova et al., (2022)', 'Brand et al., (2022)', 'Brand et al., (2023)', '', '', '', 'Lin et al., (2023)', '', '', 'VonMohr et al., (2023)', 'Makowski et al., (2023a)', 'Makowski et al., (2023b)', 'Makowski et al., (2023c)', 'Makowski et al., (2024)', 'Poerio et al., (2024)', 'Poerio et al., unpublished', 'Arjona et al., unpublished', 'Petzke et al., (2024)',  'Total'),
  Language = c('', 'English', 'English' , 'English and Italian', 'Portuguese', 'English', 'English', 'German', '', 'German', 'German', 'German', '', 'Chinese', 'Chinese', 'English', 'English', 'English', 'English', 'English', 'English', 'English','English', 'German', ''),
  N = c('', nrow(df1a), nrow(df1b), nrow(df2), nrow(df3), nrow(df4), nrow(df5_attention), nrow(df6),'', nrow(df7a), nrow(df7b), 802,'', nrow(df8a), nrow(df8b), nrow(df9), nrow(df10), nrow(df11), nrow(df12), nrow(df13), nrow(df14), nrow(df15), nrow(df16), nrow(df17), 33526),
  'Difference' = c('','', '', '', '', '', '', '', '', '', '', '', '', 'Collapsed "Itch" and "Tickling"', 'Collapsed "Itch" and "Tickling"', '', 'Analog scales. No Temperature, Blood sugar and Cough items', 'Analog scales', 'Analog scales', '', '', '', '', '', ''),
  Age = c('', make_age(df1a$Age), make_age(df1b$Age), make_age(df2$Age), make_age(df3$Age), "48.6.6 ± 14.1*", make_age(df5_attention$Age), make_age(df6$Age),'', make_age(df7a$Age), make_age(df7b$Age), make_age(df7c$Age),'', make_age(df8a$Age), make_age(df8b$Age), make_age(df9$Age), make_age(df10$Age), make_age(df11$Age), make_age(df12$Age), make_age(df13$Age), make_age(df14$Age), make_age(df15$Age), make_age(df16$Age), make_age(df17$Age), '47.96 ± 13.1'),
  Range = c('', '18-69', '18-91', '18-58', '18-72', '18-92*', '18-73', '18-78', '', '18-79', '16-81', '18-72','', '16-60', '20-60', '18-93', '18-73', '17-76','18-50', '17, 87', '18-57', '18-60', '18-79', '22-69', '17-93'),
  Female_Percentage = c('', '69.4%', '70.1%', '60.3%', '59.6%', '50%*', '46.8%', '78.7%', '', '79.5%', '77.7%', '68.9%', '', '57.0%', '56.2%', '73.2%', '50.3%', '53.0%', '76%', '57.3%', '74.8%', '75.9%', '67.7%', '68.5%', '71.3%'),
  Availability= c('osf.io/3m5nh', '', '', 'osf.io/5x9sg', 'osf.io/j6ef3', 'osf.io/ms354', 'osf.io/mp3cy', 'osf.io/xwz6g', 'osf.io/3f2h6', '', '', '','osf.io/3eztd','', '', 'osf.io/7p9u5', 'github.com/RealityBending/IllusionGameReliability', 'github.com/DominiqueMakowski/PHQ4R', 'github.com/RealityBending/InteroceptionPrimals', 
'github.com/RealityBending/InteroceptionScale',
'osf.io/49wbv', '','', 'osf.io/seru4', '')
) 

table_apa <- table |> 
  gt() |>
  cols_align(align = c("right"), columns = "Age") |> 
  cols_label(Age = "Age (Mean  ± SD)", Female_Percentage = "Female %") |> 
  tab_footnote("Note.*Information taken from the sample description of relevant paper rather than recomputed.")
  
gt_apastyle(table_apa, font.size=12)
# gtsave(gt_apastyle(table_apa, font.size=9), "figures/table1.tex")
saveRDS(table, "figures/table1.rds")

# gtsave(table_apa, "figures/table1.png")

Item Selection

Distribution

The items with the differing distribution pattern (with a lower mode around 2/5) are “Affective Touch”, “Blood Sugar” and “Bruise”.

Code
dens1a <- estimate_density(select(df1a, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1a")
dens1b <- estimate_density(select(df1b, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 1b")
dens2 <- estimate_density(select(df2, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 2")
dens3 <- estimate_density(select(df3, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 3")
dens4 <- estimate_density(select(df4, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 4")
dens5 <- estimate_density(select(df5, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 5")
dens6 <- estimate_density(select(df6, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 6")
dens7a <- estimate_density(select(df7a, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7a")
dens7b <- estimate_density(select(df7b, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7b")
dens7c <- estimate_density(select(df7c, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 7c")
dens8a <- estimate_density(select(df8a, all_of(setdiff(ias_vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8a")
dens8b <- estimate_density(select(df8b, all_of(setdiff(ias_vars, "Tickle"))), method = "kernSmooth") |>
  mutate(Sample = "Sample 8b")
dens9 <- estimate_density(select(df9, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 9")
dens10 <- estimate_density(select(df10, all_of(setdiff(ias_vars, c("Taste", "Cough", "Blood_Sugar")))), method = "kernSmooth") |>
  mutate(Sample = "Sample 10")
dens11<- estimate_density(select(df11, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 11")
dens12  <- estimate_density(select(df12, all_of(ias_vars)), method = "kernSmooth") |>
    mutate(Sample = "Sample 12")
dens13  <- estimate_density(select(df13, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 13")
dens14  <- estimate_density(select(df14, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 14")
dens15  <- estimate_density(select(df15, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 15")
dens16  <- estimate_density(select(df16, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 16")
dens17  <- estimate_density(select(df17, all_of(ias_vars)), method = "kernSmooth") |>
  mutate(Sample = "Sample 17")


p1_data <- rbind(dens1a, dens1b, dens2, dens3, dens4, dens5, dens6, dens7a, dens7b, dens7c, dens8a, dens8b, dens9, dens10, dens11 , dens12, dens13, dens14, dens15, dens16, dens17) 

linetype <- setNames(rep("solid", length(ias_vars)), ias_vars)
linetype["Affective touch"] <- "dotted"
linetype["Blood sugar"] <- "dashed"
linetype["Sex arousal"] <- "solid"
linetype["Bruise"] <- "dashed"

# proportion of non-1 and non-5 values 
p1_data |>
  filter(x > 1 & x < 5) |>  # Keep only x values strictly between 1 and 5
  summarise(proportion = sum(y) / sum(p1_data$y))  # Normalize by total density
  proportion
1  0.9996777
Code
p1 <- p1_data |>
  mutate(Sample = fct_relevel(Sample, "Sample 1a", "Sample 1b", "Sample 2", "Sample 3", "Sample 4", "Sample 5", "Sample 6", "Sample 7a", "Sample 7b", "Sample 7c", "Sample 8a", "Sample 8b", "Sample 9", "Sample 10", "Sample 11", "Sample 12", "Sample 13", "Sample 14", "Sample 15", 'Sample 16', 'Sample 17'),
         Parameter = ifelse(Parameter == "Affective_touch", "Affective touch", Parameter),
         Parameter = ifelse(Parameter == "Blood_Sugar", "Blood sugar", Parameter),
         Parameter = ifelse(Parameter == "Sex_arousal", "Sex arousal", Parameter)) |>
  # mutate(Dashed = ifelse(Parameter %in% c("Affective_touch", "Blood_Sugar", "Bruise"), TRUE, FALSE)) |> 
  ggplot(aes(x = x, y = y, color = Parameter)) +
  geom_line(aes(linetype = Parameter), linewidth = 0.5) +
  scale_color_material_d() +
  scale_linetype_manual(values = linetype)  +
  labs(x = "Response", y = "Distribution", color = "Item", linetype = "Item", title = "Item Distribution") +
  guides(color = guide_legend(ncol = 1)) +
  facet_wrap(~Sample, scales = "free_y", nrow = 3) +
  theme_minimal() +
  theme(axis.text.y = element_blank(),
        axis.title.y = element_blank(),
        plot.title = element_text(face="bold"))

p1

Code
#  ONLY IAS DATA
data1a <- select(df1a, all_of(dens1a$Parameter))
data1b <- select(df1b, all_of(dens1b$Parameter))
data2 <- select(df2, all_of(dens2$Parameter))
data3 <- select(df3, all_of(dens3$Parameter))
data4 <- select(df4, all_of(dens4$Parameter))
data5 <- select(df5, all_of(dens5$Parameter))
data6 <- select(df6, all_of(dens6$Parameter))
data7a <- select(df7a, all_of(dens7a$Parameter))
data7b <- select(df7b, all_of(dens7b$Parameter))
data7c <- select(df7c, all_of(dens7c$Parameter))
data8a <- select(df8a, all_of(dens8a$Parameter))
data8b <- select(df8b, all_of(dens8b$Parameter))
data9 <- select(df9, all_of(dens9$Parameter))
data10 <- select(df10, all_of(dens10$Parameter))
data11  <- select(df11, all_of(dens11$Parameter))
data12  <- select(df12, all_of(dens12$Parameter))
data13  <- select(df13, all_of(dens13$Parameter))
data14  <- select(df14, all_of(dens14$Parameter))
data15   <- select(df15, all_of(dens15$Parameter))
data16   <- select(df16, all_of(dens16$Parameter))
data17   <- select(df17, all_of(dens17$Parameter))


data_all <- rbind(
  data1a, data1b, data2, data3, data4, data5, data6, data7a, data7b, data7c,
  mutate(data8a, Tickle = NA), 
  mutate(data8b, Tickle = NA), data9,
  mutate(data10, Taste = NA, Cough = NA, Blood_Sugar = NA),  
  data11,  data12 , data13, data14, data15, data16, data17  
) 

Correlations

An overall positive intercorrelation pattern, with no clear structure emerging.

Code
make_correlation <- function(df) {
  correlation::correlation(df, redundant = TRUE) |>
    correlation::cor_sort() |>
    # correlation::cor_lower() |>
    mutate(val = paste0(insight::format_value(r), format_p(p, stars_only = TRUE)),
           r = ifelse(Parameter1 == Parameter2, NA, r),
           lbl = ifelse(Parameter1 == Parameter2, "", format_value(r, zap_small = TRUE, lead_zero = FALSE)),
           lbl = ifelse(p > .05, "", lbl),
           Param1 = fct_relevel(str_replace(Parameter1, "_", " "), str_replace(levels(Parameter1), "_", " ")),
           Param2 = fct_relevel(str_replace(Parameter2, "_", " "), str_replace(levels(Parameter2), "_", " "))) |>
    # mutate(Parameter2 = fct_rev(Parameter2)) |>
    ggplot(aes(x = Param1, y = Param2)) +
    geom_tile(aes(fill = r), color = "white") +
    geom_text(aes(label = lbl), size = 3) +
    labs(title = "Correlation Matrix", subtitle = paste0("N = ", nrow(df))) +
    scale_fill_gradientn(
      colors = c("white", "#FFF9C4", "#FFF59D", "#FFEB3B", "#FFCA28", "#FF9800", "#FF5722",
        "#F44336", "#E91E63", "#9C27B0", "#673AB7", "#512DA8", "#311B92"),
      na.value = "#263238",
      limits = c(0.05, 0.75),
      guide = guide_colorbar(ticks.colour = NA)
    ) +
    # scale_fill_metro_c(limit = c(0, 0.75), guide = guide_colourbar(ticks = FALSE)) +
    theme_minimal() +
    theme(
      legend.title = element_blank(),
      legend.position = "none",
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.title = element_text(face = "bold")
    )
}

p2 <- make_correlation(data_all)
p2

Code
fig1 <- p1 / p2
ggsave("figures/Figure1.png", fig1, width=10, height=12, dpi=300, bg="white")
Code
make_correlation(data1a)

Code
make_correlation(data1b)

Code
make_correlation(data2)

Code
make_correlation(data3)

Code
make_correlation(data6)

Code
make_correlation(data5)

Code
make_correlation(data6)

Code
make_correlation(data7a)

Code
make_correlation(data7b)

Code
make_correlation(data7c)

Code
make_correlation(data8a)

Code
make_correlation(data8b)

Code
make_correlation(data9)

Code
make_correlation(data10)

Code
make_correlation(data11)

Code
make_correlation(data12)

Code
make_correlation(data13)

Code
make_correlation(data14)

Code
make_correlation(data15)

Code
make_correlation(data16)

Code
make_correlation(data17)

Unique Variable Analysis (UVA)

Unique Variable Analysis (Christensen, Garrido, & Golino, 2023) uses the weighted topological overlap measure (Nowick et al., 2009) on an estimated network. Values greater than 0.25 are determined to have considerable local dependence (i.e., redundancy) that should be handled (variables with the highest maximum weighted topological overlap to all other variables (other than the one it is redundant with) should be removed).

Code
uva0 <- EGAnet::UVA(data = data_all, cut.off = 0.3)
uva0
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.361

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.288
 Urinate Defecate 0.265

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Sneeze     Cough 0.228
  Heart Breathing 0.227
 Hungry   Thirsty 0.217
Code
uva0$keep_remove
$keep
[1] "Itch"

$remove
[1] "Tickle"
Code
uva1a <- EGAnet::UVA(data = data1a, cut.off = 0.3)
uva1a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.286
   Wind   Burp 0.275

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.244
 Urinate Defecate 0.241
Code
uva1a$keep_remove
NULL
Code
uva1b <- EGAnet::UVA(data = data1b, cut.off = 0.3)
uva1b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.350
 Sneeze  Cough 0.317
   Wind   Burp 0.309

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva1b$keep_remove
$keep
[1] "Cough" "Burp"  "Itch" 

$remove
[1] "Sneeze" "Wind"   "Tickle"
Code
uva2 <- EGAnet::UVA(data = data2, cut.off = 0.3)
uva2
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Tickle      Itch 0.283
  Heart Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Sneeze    Cough 0.230
 Urinate Defecate 0.223
  Hungry  Thirsty 0.200
Code
uva2$keep_remove
NULL
Code
uva3 <- EGAnet::UVA(data = data3, cut.off = 0.3)
uva3
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.445
 Sneeze  Cough 0.318

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
   Wind   Burp 0.256

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Bruise Blood_Sugar 0.219
 Urinate    Defecate 0.217
   Heart   Breathing 0.208
Code
uva3$keep_remove
$keep
[1] "Sneeze" "Itch"  

$remove
[1] "Cough"  "Tickle"
Code
uva4 <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva4
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237
Code
uva4$keep_remove
$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  
Code
uva5 <- EGAnet::UVA(data = data5, cut.off = 0.3)
Warning: Some variables did not belong to a dimension: Blood_Sugar, Affective_touch 

Use caution: These variables have been removed from the TEFI calculation
Code
uva5
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

    node_i node_j   wto
    Tickle   Itch 0.247
    Bruise   Itch 0.246
 Breathing  Vomit 0.207
Code
uva5$keep_remove
NULL
Code
uva6 <- EGAnet::UVA(data = data6, cut.off = 0.3)
uva6
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.416
  Sneeze    Cough 0.332
    Wind     Burp 0.314
  Tickle     Itch 0.303

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Hungry Thirsty 0.274

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.237
Code
uva6$keep_remove
$keep
[1] "Urinate" "Sneeze"  "Wind"    "Itch"   

$remove
[1] "Defecate" "Cough"    "Burp"     "Tickle"  
Code
uva7a <- EGAnet::UVA(data = data7a, cut.off = 0.3)
uva7a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.434
 Tickle   Itch 0.376

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.295
 Urinate Defecate 0.294

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva7a$keep_remove
$keep
[1] "Sneeze" "Tickle"

$remove
[1] "Cough" "Itch" 
Code
uva7b <- EGAnet::UVA(data = data7b, cut.off = 0.3)
uva7b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.317
  Tickle     Itch 0.311

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.281
 Sneeze     Cough 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.237
Code
uva7b$keep_remove
$keep
[1] "Defecate" "Itch"    

$remove
[1] "Urinate" "Tickle" 
Code
uva7c <- EGAnet::UVA(data = data7c, cut.off = 0.3)
uva7c
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.457
 Sneeze  Cough 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
    Wind     Burp 0.280
 Urinate Defecate 0.254

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
 Hungry   Thirsty 0.242
  Vomit    Sneeze 0.224
  Heart Breathing 0.218
Code
uva7c$keep_remove
$keep
[1] "Cough" "Itch" 

$remove
[1] "Sneeze" "Tickle"
Code
uva8a <- EGAnet::UVA(data = data8a, cut.off = 0.3)
uva8a
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.237
   Heart Breathing 0.218
  Hungry   Thirsty 0.213
Code
uva8a$keep_remove
NULL
Code
uva8b <- EGAnet::UVA(data = data8b, cut.off = 0.3)
uva8b
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.277
   Heart Breathing 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
   Wind   Burp 0.206
Code
uva8b$keep_remove
NULL
Code
uva9 <- EGAnet::UVA(data = data9, cut.off = 0.3)
uva9
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.379
    Wind     Burp 0.373
 Urinate Defecate 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.285
 Sneeze     Cough 0.278

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j  wto
 Hungry Thirsty 0.24
Code
uva9$keep_remove
$keep
[1] "Defecate" "Wind"     "Itch"    

$remove
[1] "Urinate" "Burp"    "Tickle" 
Code
uva10 <- EGAnet::UVA(data = data10, cut.off = 0.3)
uva10
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i node_j   wto
 Tickle   Itch 0.297

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i    node_j   wto
  Heart Breathing 0.209
Code
uva10$keep_remove
NULL
Code
uva11  <- EGAnet::UVA(data = data11 , cut.off = 0.3)
uva11 
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva11$keep_remove
NULL
Code
uva12  <- EGAnet::UVA(data = data12, cut.off = 0.3)
Warning: Some variables did not belong to a dimension: Cough 

Use caution: These variables have been removed from the TEFI calculation
Code
uva12 
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)
Code
uva12$keep_remove
NULL
Code
uva13 <- EGAnet::UVA(data = data13, cut.off = 0.3)
uva13 
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i    node_j   wto
  Heart Breathing 0.395

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i   node_j   wto
  Tickle     Itch 0.229
 Urinate Defecate 0.216
    Wind     Burp 0.209
  Hungry  Thirsty 0.204
Code
uva13$keep_remove
$keep
[1] "Breathing"

$remove
[1] "Heart"
Code
uva14  <- EGAnet::UVA(data = data14, cut.off = 0.3)
uva14 
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

  node_i    node_j   wto
 Urinate  Defecate 0.330
   Heart Breathing 0.314
    Wind      Burp 0.302

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i  node_j   wto
 Muscles    Pain 0.279
  Tickle    Itch 0.277
  Hungry Thirsty 0.255

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i node_j   wto
  Vomit Sneeze 0.201
Code
uva14$keep_remove
$keep
[1] "Heart"    "Defecate" "Burp"    

$remove
[1] "Breathing" "Urinate"   "Wind"     
Code
uva15   <- EGAnet::UVA(data = data15, cut.off = 0.3)
uva15  
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i    node_j   wto
 Hungry Breathing 0.252

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j   wto
   Wind    Burp 0.214
   Temp Muscles 0.202
Code
uva15$keep_remove
NULL
Code
uva16  <- EGAnet::UVA(data = data16, cut.off = 0.3)
uva16  
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j  wto
   Wind   Burp 0.36

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

 node_i  node_j   wto
 Tickle    Itch 0.280
 Hungry Thirsty 0.261

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

  node_i      node_j   wto
  Sneeze       Cough 0.246
 Urinate    Defecate 0.234
  Bruise Blood_Sugar 0.214
Code
uva16$keep_remove
$keep
[1] "Wind"

$remove
[1] "Burp"
Code
uva17  <- EGAnet::UVA(data = data17, cut.off = 0.3)
uva17  
Variable pairs with wTO > 0.30 (large-to-very large redundancy)

 node_i node_j   wto
 Sneeze  Cough 0.341
 Tickle   Itch 0.341

----

Variable pairs with wTO > 0.25 (moderate-to-large redundancy)

  node_i   node_j   wto
 Urinate Defecate 0.267

----

Variable pairs with wTO > 0.20 (small-to-moderate redundancy)

 node_i  node_j   wto
   Wind    Burp 0.237
 Hungry Thirsty 0.209
Code
uva17$keep_remove
$keep
[1] "Cough"  "Tickle"

$remove
[1] "Sneeze" "Itch"  
Code
rez_uva <- c(uva1a$keep_remove$remove,
  uva1b$keep_remove$remove,
  uva2$keep_remove$remove,
  uva3$keep_remove$remove,
  uva4$keep_remove$remove,
  uva5$keep_remove$remove,
  uva6$keep_remove$remove,
  uva7a$keep_remove$remove,
  uva7b$keep_remove$remove,
  uva7c$keep_remove$remove,
  uva8a$keep_remove$remove,
  uva8b$keep_remove$remove,
  uva9$keep_remove$remove,
  uva10$keep_remove$remove,
  uva11$keep_remove$remove,
  uva12$keep_remove$remove,
  uva13$keep_remove$remove,
  uva14$keep_remove$remove,
  uva15$keep_remove$remove,
  uva16$keep_remove$remove,
  uva17$keep_remove$remove)

sort(table(rez_uva), decreasing = TRUE) / length(rez_uva) 
rez_uva
    Tickle       Burp      Cough     Sneeze    Urinate   Defecate       Itch 
0.24137931 0.13793103 0.13793103 0.10344828 0.10344828 0.06896552 0.06896552 
      Wind  Breathing      Heart 
0.06896552 0.03448276 0.03448276 

Structure Refinement

Discarded the following items: - Tickle: not present in the same dataset and consistently flagged as redundant in UVA analysis.

Code
data_all <- select(data_all, -Tickle)
data1a <- select(data1a, -Tickle)
data1b <- select(data1b, -Tickle)
data2 <- select(data2, -Tickle)
data3 <- select(data3, -Tickle)
data6 <- select(data6, -Tickle)
data5 <- select(data5, -Tickle)
# data6 <- select(data6, -Tickle)
data7a <- select(data7a, -Tickle)
data7b <- select(data7b, -Tickle)
data7c <- select(data7c, -Tickle)
# data8a <- select(data8a, -Tickle)
# data8b <- select(data8b, -Tickle)
data9 <- select(data9, -Tickle)
data10 <- select(data10, -Tickle)
data11 <- select(data11 , -Tickle)
data12 <- select(data12 , -Tickle)
data13 <- select(data13 , -Tickle)
data14 <- select(data14 , -Tickle)
data15 <- select(data15 , -Tickle)
data16 <- select(data16 , -Tickle)
data17 <- select(data17 , -Tickle)

Initial

Code
colors <- c("Heart" = "#F44336", "Breathing" = "#FF5722",
            "Hungry" = "#FF9800", "Thirsty" = "#FFC107",
            "Burp" = "#4CAF50", "Wind" = "#009688",
            "Urinate" = "#63824C", "Defecate" = "#795548",
            "Cough" = "#8BC34A", "Sneeze" = "#CDDC39", 
            "Bruise" = "#673AB7", "Blood Sugar" = "#3F51B5",
            "Muscles" = "#2196F3", "Pain" = "#00BCD4",
            "Sex arousal" = "#FF4081", "Temp" = "#9C27B0",
            "Vomit" = "#76FF03", "Taste" = "#00E676",
            "Affective touch" = "#FF1744", "Itch" = "#D500F9")
Code
make_hclust <- function(data) {
  rez_pvclust <- pvclust::pvclust(data,
                 method.hclust = "complete",
                 method.dist = "correlation",
                 nboot = 1000, quiet=TRUE, parallel=TRUE)
  
  
  dendrogram <- as.dist(1 - cor(data, use = "pairwise.complete.obs")) |> 
  hclust(method = "complete") |> 
  tidygraph::as_tbl_graph() 


  # Process Nodes
  nodes <- as.list(dendrogram)$nodes |> 
    mutate(
      Size = ifelse(label != "", 10, NA),
      Item = str_replace(label, "_", " "),
      idx = 1:nrow(as.list(dendrogram)$nodes))
  
  # Central node 
  nodes[nodes$height == max(nodes$height), c("Item", "Size")] <- data.frame(Item="Central", Size=15)
  
  # Process Edges
  edges <- as.list(dendrogram)$edges
  edges$linewidth = datawizard::rescale(nodes[edges$from, ]$height, to = c(0.1, 1))
  
  p <- tidygraph::tbl_graph(nodes = nodes, edges = edges) |> 
    ggraph(layout = "dendrogram", circular = TRUE) +
    # geom_edge_diagonal(strength = 0.7, linewidth = 1) +
    geom_edge_elbow2(aes(edge_width=linewidth), color="#212121") +
    geom_node_point(aes(filter=Item %in% c("Central"), size = Size), color="#212121") +
    # geom_node_point(aes(filter=Group %in% c("Visceroception", "Awareness", "Deficit"), color=Group, size = Size)) +
    geom_node_text(aes(
      label = ifelse(Item != "", Item, NA),
      x = x * 1.10,
      y = y * 1.10,
      filter = label != "",
      angle = ifelse(
        x >= 0,
        asin(y) * 360 / 2 / pi,
        360 - asin(y) * 360 / 2 / pi
      ),
      hjust = ifelse(
        x >= 0, 0, 1
      ))) +
    geom_node_point(aes(filter = label != "", color=Item, size=Size), alpha = 1) +
    # geom_node_text(aes(label=idx)) +  # Debug
    scale_edge_width_continuous(range=c(1, 3), guide = "none") +
    scale_size_continuous(range=c(10, 10), guide = "none") +
    scale_color_manual(values = colors,
                       breaks = names(colors)) +
    ggtitle("Hierarchical Clustering Analysis (HCA)", subtitle = "Method = Correlation") +
    coord_equal(clip = "off", xlim = c(-1.25, 1.25), ylim = c(-1.25, 1.25)) +
    theme_void() + 
    # guides(color = guide_legend(override.aes = list(size = c(7.5, 5, 5, 3.5, 7.5, 5, 5, 3.5, 3.5, 3.5, 7.5, 5, 5, 3.5)))) +
    theme(legend.text = ggtext::element_markdown(),
          legend.title = element_blank(),
          legend.position = "none",
          # plot.title = element_blank(), 
          # plot.subtitle = element_blank()
          plot.title = element_text(face="bold"))

  
  list(pvclust = rez_pvclust, p = p)
}


rez_hclust <- make_hclust(data_all)
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
rez_hclust$p

Code
make_ega <- function(data, cols = c("C1" = "#F44336", "C2" = "#FF9800", "C3" = "#795548", "C4" = "#CDDC39", "C5" = "#2196F3", "C6" = "#009688", "C7" = "#9C27B0")) {
  rez <- EGAnet::bootEGA(
    data = data,
    seed = 123,
    model = "glasso",
    algorithm = "leiden",
    EGA.type = "hierEGA",
    type = "resampling",
    plot.itemStability = FALSE,
    verbose = FALSE, 
    allow.singleton = TRUE,
    loading.method = "original"
  )
  
  table <- EGAnet::net.loads(rez$EGA$lower_order)$std |> 
    as.data.frame() |> 
    rownames_to_column("Item") |>
    gt::gt() |>
    gt::tab_header(title = "EGA Loadings") |>
    gt::data_color(
      columns = -Item,
      method = "numeric",
      colors = scales::col_numeric(
        palette = c("red", "white", "green"),
        domain = c(-1, 0, 1)
      )) |> 
    gt::fmt_auto()
  
    
  nodes <- rez$EGA$lower_order$dim.variables |> 
    rename(name = items) |> 
    mutate(dimension = paste0("C", dimension),
           size = apply(EGAnet::net.loads(rez$EGA$lower_order)$std, 1, max)[name])
  
  loadings <- rez$EGA$lower_order$network |> 
    as.data.frame() |> 
    rownames_to_column("to") |>
    pivot_longer(-c(to), names_to = "from", values_to = "Loading") |> 
    filter(Loading > quantile(Loading, 1/3)) 
  
  g <- tidygraph::tbl_graph(nodes = nodes, edges = loadings, directed = FALSE) |> 
    mutate(name = str_replace(name, "_", " "),
           name = ifelse(name == "Sex arousal", "Sex", name)) 
  
  set.seed(1)
  layout <- ggraph::create_layout(g, layout = "fr", weights = abs(loadings$Loading))
  
  xrange <- max(layout$x) - min(layout$x)
  yrange <- max(layout$y) - min(layout$y)
  xmin <- min(layout$x) - xrange * 0.05
  xmax <- max(layout$x) + xrange * 0.05
  ymin <- min(layout$y) - yrange * 0.05
  ymax <- max(layout$y) + yrange * 0.05
  
  p <- layout |> 
    ggraph() +
    geom_edge_bend(aes(edge_width = Loading, edge_alpha = Loading), color = "#212121",
                   strength = 0.3) +
    geom_node_point(aes(size = size, color = dimension), alpha = 0.95) +
    geom_node_text(aes(label = name), size = 3, color = "white", fontface = "bold") +
    scale_size_continuous(range = c(23, 28)) +
    scale_edge_width(range = c(0.3, 4)) +
    scale_edge_alpha(range = c(0.1, 0.9)) +
    scale_color_manual(values = cols) +
    labs(title = "Exploratory Graph Analysis (EGA)", subtitle = "Method = Leiden") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold")) +
    coord_cartesian(xlim =c(xmin, xmax), ylim = c(ymin, ymax)) 
    
  list(rez = rez, table = table, p = p)
}

rez_ega <- make_ega(data_all, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]]))

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.526 0.036 0.008 0.069 0.037 0.019 0.112
Breathing 0.526 0.189 0.044 0.18  0.058 0.002 0.002
Hungry 0.129 0.522 0.074 0.123 0     0.004 0.04 
Thirsty 0.093 0.522 0.124 0.108 0.009 0.013 0.01 
Urinate 0.02  0.185 0.564 0.148 0.038 0.017 0    
Defecate 0.043 0.057 0.564 0.153 0.115 0.041 0    
Pain 0.036 0.107 0.052 0.393 0.045 0     0.131
Temp 0.026 0.091 0.088 0.363 0.125 0.03  0.006
Muscles 0.104 0.025 0.039 0.357 0.058 0.089 0.123
Sex_arousal 0.067 0.069 0.071 0.3   0.076 0.078 0    
Affective_touch 0.036 0.001 0.035 0.279 0.05  0     0.202
Taste 0.07  0.027 0.061 0.229 0.116 0.04  0.054
Sneeze 0.032 0     0.052 0.107 0.564 0.085 0.046
Cough 0.039 0.011 0.011 0.129 0.446 0.161 0.091
Vomit 0.046 0     0.094 0.185 0.293 0.043 0.045
Wind 0.008 0.01  0.057 0.118 0.101 0.571 0.058
Burp 0.019 0.011 0.003 0.098 0.19  0.571 0.107
Bruise 0.031 0     0     0.18  0.052 0.08  0.444
Blood_Sugar 0.067 0.052 0     0.073 0.022 0.025 0.385
Itch 0.019 0     0     0.133 0.078 0.031 0.298
Code
rez_ega$p

Code
make_efa <- function(data, n = 3, cols = c("F1" = "#009688", "F2" = "#795548", "F3" = "red")) {
  n_fac <- n_factors(data)

  rez_efa <- factor_analysis(data, n = n, sort = TRUE, rotation = "oblimin")
  
  # Exctract loadings
  loadings <- attributes(rez_efa)$loadings_long |> 
    select(from = Component, to=Variable, Loading) |> 
    mutate(Type = "Loading",
           to = str_replace(to, "_", " "),
           from = str_replace(from, "MR", "F")) 
  
  # Reorder factors
  fct_order <- c("Pain", "Muscles", "Blood Sugar", "Bruise",
                 "Affective touch", "Sex arousal", "Temp", "Itch", "Tickle", 
                 "Thirsty", "Hungry", "Breathing", "Heart", "Defecate", "Urinate",
                 "Sneeze", "Cough",  "Wind", "Burp", "Vomit", "Taste")
  loadings <- loadings |> 
    mutate(to = fct_relevel(to, fct_order[fct_order %in% unique(loadings$to)])) |> 
    arrange(to) |> 
    mutate(to = as.character(to))

  # Get correlations between factors
  phi <- attributes(rez_efa)$model$Phi
  phi[upper.tri(phi, diag = TRUE)] <- NA
  phi <- phi |> 
    as.data.frame() |> 
    rownames_to_column("from") %>%
    pivot_longer(-from, names_to = "to", values_to = "Loading") |> 
    filter(!is.na(Loading)) |>
    mutate(Type = "Correlation",
           from = str_replace(from, "MR", "F"),
           to = str_replace(to, "MR", "F")) 
  
  # Reverse so that arcs are on the same side
  phi[paste(phi$from, phi$to) %in% c("F3 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F3")
  phi[paste(phi$from, phi$to) %in% c("F2 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F2")
  phi[paste(phi$from, phi$to) %in% c("F3 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F3")
  phi[paste(phi$from, phi$to) %in% c("F4 F1"), c("from", "to")] <- data.frame(from = "F1", to = "F4")
  phi[paste(phi$from, phi$to) %in% c("F4 F2"), c("from", "to")] <- data.frame(from = "F2", to = "F4")
  phi[paste(phi$from, phi$to) %in% c("F4 F3"), c("from", "to")] <- data.frame(from = "F3", to = "F4")
  
  loadings <- rbind(loadings, phi)
  
  # Make layout
  factor_names <- unique(loadings$from[loadings$Type == "Loading"])
  item_names <- unique(loadings$to[loadings$Type == "Loading"])
  
  nodes <- tibble(name = unique(c(loadings$from, loadings$to))) |>
    mutate(type = case_when(
      name %in% factor_names ~ "Factor",
      name %in% item_names ~ "Item",
      TRUE ~ "Other"
    ))
  
  y_position <- seq(0.8, 0.2, length.out = length(factor_names))
  
  layout <- nodes |>
    mutate(x = ifelse(type == "Item", 0, 1),
           y = case_when(
            name == "F1" ~ y_position[1],
            name == "F2" ~ y_position[2],
            name == "F3" ~ y_position[3],
            name == "F4" ~ y_position[4],
            name == "F5" ~ y_position[5],
            name == "F6" ~ y_position[6],
            name == "F7" ~ y_position[7],
            name == "F8" ~ y_position[8],
            .default = NA
          ))
  layout[layout$type == "Item", "y"] <- seq(0, 1, length.out = sum(layout$type == "Item"))
  
  
  g <- tidygraph::as_tbl_graph(loadings, directed = TRUE) |>
    left_join(layout, by = c("name"))
  
  p_efa <- ggraph(g, layout = as.data.frame(layout)[c("x", "y")]) + 
    geom_edge_link(
      aes(
        edge_width = abs(Loading),
        edge_alpha = abs(Loading),
        # edge_color = Loading,
        filter = Type == "Loading"
        # label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
      ),
      color = "#212121",
      arrow = arrow(length = unit(2.5, "mm"), type = "closed"),
      end_cap = circle(5, 'mm')
    ) +
      geom_edge_arc(
      aes(
        edge_width = abs(Loading),
        filter = Type == "Correlation",
        label = insight::format_value(Loading, lead_zero = FALSE, zap_small = TRUE)
      ),
      edge_color = "forestgreen",
      angle_calc = "along",
      label_dodge = unit(2.5, "mm"),
      force_flip = TRUE,
      strength = 0.2  # controls curvature
    ) +
    geom_node_label(data=filter(layout, type == "Item"), 
                    aes(fill = name), 
                    label = paste(rep(" ", 20), collapse = ""),
                    size = 3, 
                    hjust = 0.5, 
                    nudge_x = -0.05,
                    label.r = unit(0.2, "lines"),
                    label.padding = unit(0.5, "lines"),
                    label.size = 0) +
    geom_node_label(data=filter(layout, type == "Item"), 
                    aes(label = name), 
                    fill = NA,
                    size = 3, 
                    hjust = 0.5, 
                    vjust = 0.5,
                    color = "white",
                    fontface = "bold",
                    nudge_x = -0.05,
                    label.size = 0) +
    geom_node_point(data=filter(layout, type == "Factor"), 
                    aes(color = name), 
                    size = 12) +
    geom_node_text(data=filter(layout, type == "Factor"), 
                   aes(label = name), 
                   size = 4, 
                   color = "white", 
                   fontface = "bold") +
    scale_edge_width(range = c(0.3, 2)) +
    scale_edge_alpha(range = c(0.1, 1)) +
    scale_fill_manual(values = colors) +
    scale_color_manual(values = cols) +
    labs(title = "Exploratory Factor Analysis (EFA)", subtitle = "Method = Oblimin") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold")) +
    coord_cartesian(xlim = c(-0.08, 1.08)) 

  list(n_fac=n_fac, efa=rez_efa, p=p_efa)
}


rez_efa <- make_efa(data_all, n = 3, cols = c("F1" = colors[["Cough"]], 
                                              "F2" = colors[["Thirsty"]], 
                                              "F3" = colors[["Bruise"]]))
Loading required namespace: GPArotation
Code
plot(rez_efa$n_fac)

Code
rez_efa$p

Final

We discarded:

  • Taste: Lone item + unstable
  • Affective_touch: Cross-loadings + unstable
  • Vomit: Less strongly associated
  • Itch: Less strongly associated, did not form consistent cluster
  • More controversial: Temp and Sex_arousal: similar cluster yet less reliable
data_allf <- select(data_all, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data1af <- select(data1a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data1bf <- select(data1b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data2f <- select(data2, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data3f <- select(data3, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data4f <- select(data4, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data5f <- select(data5, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data6f <- select(data6, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7af <- select(data7a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7bf <- select(data7b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data7cf <- select(data7c, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data8af <- select(data8a, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data8bf <- select(data8b, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data9f <- select(data9, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data10f <- select(data10, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data11f <- select(data11, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data12f <- select(data12, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data13f <- select(data13, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data14f <- select(data14, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data15f <- select(data15, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data16f <- select(data16, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
data17f <- select(data17, -Taste, -Affective_touch, -Vomit, -Itch,
                   -Sex_arousal, -Temp)
Code
rez_hclust <- make_hclust(data_allf)

plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
rez_hclust$p

Code
rez_ega <- make_ega(data_allf, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.529 0.04  0.017 0.04  0.026 0.046 0.131
Breathing 0.529 0.199 0.062 0.056 0.011 0.169 0.017
Hungry 0.139 0.526 0.087 0.008 0.013 0.096 0.056
Thirsty 0.096 0.526 0.129 0.017 0.02  0.107 0.017
Urinate 0.031 0.197 0.572 0.043 0.027 0.114 0    
Defecate 0.065 0.072 0.572 0.087 0.06  0.113 0    
Sneeze 0.058 0.01  0.088 0.561 0.107 0.091 0.028
Cough 0.054 0.02  0.036 0.561 0.174 0.115 0.084
Wind 0.014 0.025 0.072 0.116 0.575 0.079 0.065
Burp 0.032 0.017 0.016 0.184 0.575 0.095 0.125
Muscles 0.126 0.037 0.067 0.083 0.107 0.493 0.154
Pain 0.053 0.134 0.087 0.065 0.01  0.493 0.19 
Bruise 0.042 0     0     0.055 0.09  0.276 0.489
Blood_Sugar 0.077 0.06  0     0.022 0.035 0.059 0.489
Code
rez_ega$p

Code
rez_efa <- make_efa(data_allf, n = 3, cols = 
                      c("F1" = colors[["Wind"]], "F2" = colors[["Defecate"]], 
                        "F3" = colors[["Heart"]]))
plot(rez_efa$n_fac)

Code
rez_efa$p

Code
make_cfa <- function(data) {
  
  # Main models
  m_g <- '
IAS =~ Burp + Wind + Cough + Sneeze + Hungry + Thirsty  + Urinate + Defecate + Heart + Breathing + Pain + Muscles 
  '
  m_ega <- '
    CoughSneeze =~ Cough + Sneeze
    UrinateDefecate =~ Urinate + Defecate
    BruiseBloodsugar =~ Bruise + Blood_Sugar
    WindBurp =~ Wind + Burp
    HeartBreath =~ Heart + Breathing
    MusclesPain =~ Muscles + Pain
    HungryThirsty =~ Hungry + Thirsty
  '
  
  m_efa <- '
    UrinateDefecate =~ Urinate + Defecate
    Sudden =~ Cough + Sneeze + Wind + Burp
    Homeostasis =~ Heart + Breathing + Muscles + Pain + Hungry + Thirsty + Bruise + Blood_Sugar
  '
    
    
  m_hclust <- '
    HeartBreath =~ Heart + Breathing
    BruiseBloodsugar =~ Bruise + Blood_Sugar
    HungryThirsty =~ Hungry + Thirsty
    Sudden =~ Cough + Sneeze + Wind + Burp
    Unpleasant =~ Urinate + Defecate + Muscles + Pain
  '
  
  # Test bifactor models
  m_egag <- paste0(m_ega, "\nIAS =~ CoughSneeze + WindBurp + BruiseBloodsugar + UrinateDefecate + HeartBreath + MusclesPain")
  m_efag <- paste0(m_efa, "\nIAS =~ UrinateDefecate + Sudden + Homeostasis")
  m_hclustg <- paste0(m_hclust, "\nIAS =~ HeartBreath + BruiseBloodsugar + HungryThirsty + Sudden + Unpleasant")

  fit_g <- lavaan::cfa(m_g, data = data, std.lv = TRUE)
  fit_ega <- lavaan::cfa(m_ega, data = data, std.lv = TRUE)
  fit_egag <- lavaan::cfa(m_egag, data = data, std.lv = TRUE)
  fit_efa <- lavaan::cfa(m_efa, data = data, std.lv = TRUE)
  fit_efag <- lavaan::cfa(m_efag, data = data, std.lv = TRUE)
  fit_hclust <- lavaan::cfa(m_hclust, data = data, std.lv = TRUE)
  fit_hclustg <- lavaan::cfa(m_hclustg, data = data, std.lv = TRUE)
  
  # Test
  bf <- test_bf(fit_ega, fit_g, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg) |> 
    as.data.frame() |> 
    select(Name = Model, BF)
  
  # test <- lavaan::lavTestLRT(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg)
  test <- compare_performance(fit_g, fit_ega, fit_egag, fit_efa, fit_efag, fit_hclust, fit_hclustg,
                      metrics = c("BIC", "RMSEA", "Chi2", "CFI")) |> 
    select(-BIC_wt, -Model) |> 
    merge(bf, by = "Name") |>
    arrange(BIC) |> 
    gt::gt() |> 
    gt::tab_header(title = "CFA Model Comparison") |>
    gt::data_color(
      columns = c(BIC, Chi2, RMSEA),
      method = "numeric",
      palette = c("green", "white")) |>
    gt::data_color(
      columns = c(CFI),
      method = "numeric",
      palette = c("white", "green")) |>
    gt::data_color(
      columns = c(BF),
      method = "numeric",
      palette = c("red", "yellow", "white", "yellow", "green"),
      domain = c(0, 0.3, 1, 3, 10000000)) |>
    gt::fmt_number(columns = c(RMSEA, Chi2, CFI), decimals = 3) 

  
  param  <- parameters::parameters(fit_ega)
  perf <- performance::performance(fit_ega)
  
  # Graph
  nodes <- tidySEM::get_nodes(fit_ega) |> 
    mutate(name = str_replace(name, "_", " "))
  edges <- tidySEM::get_edges(fit_ega) |> 
    mutate(from = str_replace(from, "_", " "),
           to = str_replace(to, "_", " "), 
           est_std = as.numeric(est_std)) |> 
    filter(from != to)
  
  g <- tidygraph::tbl_graph(nodes = nodes, edges = edges, directed = TRUE) 
  set.seed(1)
  layout <- ggraph::create_layout(g, layout = "fr", weights = est_std)
  
  p_cfa <- ggraph(layout) +
    geom_edge_bend(aes(edge_width = abs(est_std), edge_alpha = abs(est_std),
                       label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
                       filter = arrow == "none"), 
                   color = "forestgreen",
                   angle_calc = "along", 
                   check_overlap = TRUE,
                   label_dodge = unit(0, "mm"),
                   label_size = unit(3, "mm"),
                   strength = 0.3) +
    geom_edge_link(aes(edge_width = abs(est_std), edge_color = abs(est_std), 
                       label = format_value(est_std, zap_small = TRUE, lead_zero = FALSE),
                      filter = arrow != "none"), 
                   angle_calc = "along", 
                   label_dodge = unit(3, "mm"),
                   arrow = arrow(length = unit(2.5, "mm"), type = "open", ends = "first"), 
                   lineend = "round",
                   linejoin = "bevel",
                   start_cap = circle(9, 'mm')) +
    geom_node_point(aes(color = name, shape=shape), alpha = 0.97, size = 20) +
    geom_node_text(data = filter(layout, shape == "rect"), aes(label = name), 
                   size = 3, color = "white", fontface = "bold") +
    scale_size_continuous(range = c(10, 20)) +
    scale_edge_width(range = c(0.3, 4)) +
    scale_edge_alpha(range = c(0.1, 0.9)) +
    scale_edge_color_gradient(low = "#EEEEEE", high = "#212121") +
    scale_color_manual(values = c(colors, 
                               "CoughSneeze" = "#CDDC39", 
                               "UrinateDefecate" = "#795548", "BruiseBloodsugar" = "#673AB7",
                               "WindBurp" = "#009688", "HeartBreath" = "#F44336",
                               "MusclesPain" = "#2196F3", "HungryThirsty" = "#FF9800")) +
    scale_shape_manual(values = c("oval" = 16, "rect" = 15)) +
    labs(title = "Confirmatory Factor Analysis (CFA)", subtitle = "Method = Maximum Likelihood") +
    theme_void() +
    theme(legend.position = "none",
          plot.title = element_text(face="bold"))
  p_cfa
  
  list(test = test, param = param, perf = perf, p=p_cfa)
}


rez_cfa0 <- make_cfa(data_allf)
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: Could not access model information.
Warning: When comparing models, please note that probably not all models were fit
  from same data.
When comparing models, please note that probably not all models were fit
  from same data.
Warning: Some values were outside the color scale and will be treated as NA
Code
rez_cfa0$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 940941.8 0.124 27,398.095 0.783 Inf
fit_ega 1110166.4 0.035 2,334.112 0.984 NA
fit_egag 1114563.1 0.054 6,876.395 0.952 0
fit_hclust 1121159.8 0.078 13,441.878 0.906 0
fit_hclustg 1122347.5 0.079 14,681.532 0.897 0
fit_efag 1124348.4 0.083 16,703.253 0.883 0
fit_efa 1124348.4 0.083 16,703.253 0.883 0
Code
rez_cfa0$perf
# Indices of model performance

Chi2(56) | p (Chi2) | Baseline(91) | p (Baseline) |   GFI |  AGFI |   NFI
-------------------------------------------------------------------------
2334.112 |   < .001 |    1.425e+05 |       < .001 | 0.990 | 0.981 | 0.984

Chi2(56) |  NNFI |   CFI | RMSEA |      RMSEA  CI | p (RMSEA) |   RMR |  SRMR
-----------------------------------------------------------------------------
2334.112 | 0.974 | 0.984 | 0.035 | [0.034, 0.036] |    > .999 | 0.019 | 0.019

Chi2(56) |   RFI |  PNFI |   IFI |   RNI | Loglikelihood |     AIC |     BIC | BIC_adjusted
-------------------------------------------------------------------------------------------
2334.112 | 0.973 | 0.605 | 0.984 | 0.984 |    -5.548e+05 | 1.1e+06 | 1.1e+06 |    1.110e+06
Code
rez_cfa0$param
# Loading

Link                            | Coefficient |       SE |       95% CI
-----------------------------------------------------------------------
CoughSneeze =~ Cough            |        0.68 | 4.57e-03 | [0.67, 0.69]
CoughSneeze =~ Sneeze           |        0.60 | 4.32e-03 | [0.59, 0.61]
UrinateDefecate =~ Urinate      |        0.67 | 4.85e-03 | [0.66, 0.68]
UrinateDefecate =~ Defecate     |        0.63 | 4.55e-03 | [0.62, 0.64]
BruiseBloodsugar =~ Bruise      |        0.86 | 8.34e-03 | [0.85, 0.88]
BruiseBloodsugar =~ Blood_Sugar |        0.63 | 7.65e-03 | [0.61, 0.64]
WindBurp =~ Wind                |        0.72 | 5.05e-03 | [0.71, 0.73]
WindBurp =~ Burp                |        0.80 | 5.19e-03 | [0.79, 0.81]
HeartBreath =~ Heart            |        0.57 | 5.82e-03 | [0.56, 0.58]
HeartBreath =~ Breathing        |        0.69 | 5.68e-03 | [0.68, 0.70]
MusclesPain =~ Muscles          |        0.58 | 4.63e-03 | [0.57, 0.59]
MusclesPain =~ Pain             |        0.58 | 4.83e-03 | [0.57, 0.59]
HungryThirsty =~ Hungry         |        0.73 | 6.34e-03 | [0.72, 0.74]
HungryThirsty =~ Thirsty        |        0.73 | 6.19e-03 | [0.72, 0.74]

Link                            |      z |      p
-------------------------------------------------
CoughSneeze =~ Cough            | 148.76 | < .001
CoughSneeze =~ Sneeze           | 138.72 | < .001
UrinateDefecate =~ Urinate      | 139.04 | < .001
UrinateDefecate =~ Defecate     | 138.08 | < .001
BruiseBloodsugar =~ Bruise      | 103.66 | < .001
BruiseBloodsugar =~ Blood_Sugar |  82.24 | < .001
WindBurp =~ Wind                | 143.16 | < .001
WindBurp =~ Burp                | 153.21 | < .001
HeartBreath =~ Heart            |  97.94 | < .001
HeartBreath =~ Breathing        | 121.97 | < .001
MusclesPain =~ Muscles          | 126.04 | < .001
MusclesPain =~ Pain             | 120.36 | < .001
HungryThirsty =~ Hungry         | 114.88 | < .001
HungryThirsty =~ Thirsty        | 118.19 | < .001

# Correlation

Link                                | Coefficient |       SE |       95% CI
---------------------------------------------------------------------------
CoughSneeze ~~ UrinateDefecate      |        0.56 | 5.86e-03 | [0.55, 0.57]
CoughSneeze ~~ BruiseBloodsugar     |        0.55 | 7.21e-03 | [0.54, 0.56]
CoughSneeze ~~ WindBurp             |        0.72 | 4.73e-03 | [0.72, 0.73]
CoughSneeze ~~ HeartBreath          |        0.52 | 6.61e-03 | [0.51, 0.53]
CoughSneeze ~~ MusclesPain          |        0.66 | 5.94e-03 | [0.65, 0.68]
CoughSneeze ~~ HungryThirsty        |        0.46 | 6.89e-03 | [0.45, 0.48]
UrinateDefecate ~~ BruiseBloodsugar |        0.40 | 7.70e-03 | [0.38, 0.41]
UrinateDefecate ~~ WindBurp         |        0.51 | 5.99e-03 | [0.50, 0.52]
UrinateDefecate ~~ HeartBreath      |        0.53 | 6.60e-03 | [0.52, 0.54]
UrinateDefecate ~~ MusclesPain      |        0.64 | 6.07e-03 | [0.63, 0.66]
UrinateDefecate ~~ HungryThirsty    |        0.66 | 5.99e-03 | [0.64, 0.67]
BruiseBloodsugar ~~ WindBurp        |        0.58 | 6.99e-03 | [0.56, 0.59]
BruiseBloodsugar ~~ HeartBreath     |        0.50 | 7.93e-03 | [0.48, 0.52]
BruiseBloodsugar ~~ MusclesPain     |        0.73 | 7.30e-03 | [0.72, 0.75]
BruiseBloodsugar ~~ HungryThirsty   |        0.45 | 8.13e-03 | [0.44, 0.47]
WindBurp ~~ HeartBreath             |        0.46 | 6.73e-03 | [0.44, 0.47]
WindBurp ~~ MusclesPain             |        0.63 | 5.96e-03 | [0.62, 0.65]
WindBurp ~~ HungryThirsty           |        0.44 | 6.86e-03 | [0.43, 0.45]
HeartBreath ~~ MusclesPain          |        0.64 | 6.76e-03 | [0.63, 0.66]
HeartBreath ~~ HungryThirsty        |        0.64 | 6.75e-03 | [0.63, 0.65]
MusclesPain ~~ HungryThirsty        |        0.63 | 6.77e-03 | [0.62, 0.65]

Link                                |      z |      p
-----------------------------------------------------
CoughSneeze ~~ UrinateDefecate      |  95.13 | < .001
CoughSneeze ~~ BruiseBloodsugar     |  76.32 | < .001
CoughSneeze ~~ WindBurp             | 153.13 | < .001
CoughSneeze ~~ HeartBreath          |  78.74 | < .001
CoughSneeze ~~ MusclesPain          | 111.88 | < .001
CoughSneeze ~~ HungryThirsty        |  67.02 | < .001
UrinateDefecate ~~ BruiseBloodsugar |  51.46 | < .001
UrinateDefecate ~~ WindBurp         |  84.78 | < .001
UrinateDefecate ~~ HeartBreath      |  80.33 | < .001
UrinateDefecate ~~ MusclesPain      | 106.28 | < .001
UrinateDefecate ~~ HungryThirsty    | 109.51 | < .001
BruiseBloodsugar ~~ WindBurp        |  82.78 | < .001
BruiseBloodsugar ~~ HeartBreath     |  63.04 | < .001
BruiseBloodsugar ~~ MusclesPain     | 100.58 | < .001
BruiseBloodsugar ~~ HungryThirsty   |  55.88 | < .001
WindBurp ~~ HeartBreath             |  67.88 | < .001
WindBurp ~~ MusclesPain             | 106.23 | < .001
WindBurp ~~ HungryThirsty           |  64.04 | < .001
HeartBreath ~~ MusclesPain          |  95.23 | < .001
HeartBreath ~~ HungryThirsty        |  94.69 | < .001
MusclesPain ~~ HungryThirsty        |  93.63 | < .001
Code
rez_cfa0$p

Figure

The structural analysis seem to converge on the idea of small clusters (of pairs or triplets) that are potentially inter-related parts of larger clusters (although with unstable associations). Best way to assess the values of this new granular structure is to assess its predictive value against convergent measures, and see if it’s superior to a unique score (-> Study 2).

Code
fig2 <- (rez_hclust$p | rez_ega$p) / (rez_efa$p | rez_cfa0$p) 
ggsave("figures/Figure2.png", fig2, width=14, height=14, dpi=300, bg="white")

Structure Invariance

Initial

Code
rez_hclust <- make_hclust(data1a)
rez_ega <- make_ega(data1a, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1a, n = 4, cols = c("F1" = colors[["Itch"]], "F2" = colors[["Burp"]], "F3" = colors[["Defecate"]], "F4" = colors[["Hungry"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7 8
Heart 0.504 0.125  0.046 0     0.093 0.114 0      0    
Breathing 0.504 0.14   0.032 0.098 0.158 0.123 0      0.036
Hungry 0.2   0.511  0.12  0.004 0.028 0.064 0      0.034
Thirsty 0.074 0.511  0.166 0     0.068 0.13  0.005  0.053
Defecate 0.044 0.093  0.632 0.107 0.077 0.036 0.059 −0.045
Urinate 0.026 0.253  0.421 0.035 0.167 0.035 0.046 −0.02 
Taste 0.032 0.018  0.21  0.105 0.18  0.149 0.082  0.028
Sneeze 0     0.006  0.078 0.676 0.094 0.001 0.092  0.044
Cough 0.053 0      0.014 0.377 0.11  0     0.171  0.121
Vomit 0.088 0      0.177 0.299 0.066 0.067 0      0.054
Muscles 0.084 0.022  0.13  0.03  0.454 0.032 0.029  0.093
Pain 0.126 0.032  0.133 0.067 0.438 0.195 0      0.143
Temp 0.086 0.056  0.117 0.125 0.264 0.245 0      0.038
Sex_arousal 0.103 0.042  0.055 0.006 0.069 0.435 0.119  0.007
Affective_touch 0.059 0.086  0.058 0.027 0.205 0.435 0.011  0.152
Wind 0     0.007  0.189 0.051 0     0.261 0.571  0.058
Burp 0     0      0.007 0.203 0.034 0     0.571  0.105
Blood_Sugar 0.013 0     −0.02  0.04  0.033 0.033 0.015  0.455
Bruise 0     0.108 −0.044 0.08  0.193 0.044 0.068  0.401
Itch 0.034 0      0.027 0.076 0.074 0.224 0.07   0.388
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data1af)
rez_ega <- make_ega(data1af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1af, n = 3, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Urinate"]], "F3" = colors[["Bruise"]]))
rez_cfa1a  <- make_cfa(data1af)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.511 0.128  0.036 0.009 0.01  0.108  0    
Breathing 0.511 0.16   0.048 0.068 0.015 0.151  0.045
Hungry 0.202 0.514  0.105 0.015 0     0.042  0.048
Thirsty 0.091 0.514  0.178 0     0.018 0.041  0.089
Urinate 0.029 0.253  0.565 0.038 0.05  0.185 −0.017
Defecate 0.08  0.108  0.565 0.098 0.087 0.09  −0.042
Sneeze 0.033 0.02   0.12  0.574 0.117 0.091  0    
Cough 0.071 0      0.021 0.574 0.171 0.06   0.157
Wind 0.021 0.024  0.136 0.071 0.577 0.028  0.089
Burp 0.012 0      0.008 0.22  0.577 0.041  0.089
Muscles 0.091 0.042  0.118 0.058 0.04  0.523  0.125
Pain 0.183 0.044  0.107 0.061 0.014 0.523  0.18 
Bruise 0.01  0.128 −0.031 0.059 0.093 0.222  0.502
Blood_Sugar 0.033 0     −0.013 0.052 0.032 0.053  0.502
Code
plot(rez_efa$n_fac)

Code
rez_cfa1a$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 13466.62 0.106 325.704 0.784 Inf
fit_ega 15870.94 0.026 72.825 0.988 NA
fit_egag 15892.96 0.059 180.403 0.924 1.654961e-05
fit_hclustg 15911.16 0.065 210.828 0.905 1.846821e-09
fit_hclust 15913.39 0.062 182.492 0.921 6.080789e-10
fit_efa 15941.59 0.073 253.472 0.877 4.576430e-16
fit_efag 15941.59 0.073 253.472 0.877 4.576430e-16

Initial

Code
rez_hclust <- make_hclust(data1b)
rez_ega <- make_ega(data1b, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Defecate"]], "C3" = colors[["Bruise"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Temp"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1b, n = 4, cols = c("F1" = colors[["Urinate"]], "F2" = colors[["Bruise"]], 
                                            "F3" = colors[["Sneeze"]], "F4" = colors[["Wind"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5
Breathing 0.42  0.118 0.065 0.016 0.037
Hungry 0.415 0.078 0.13  0     0.023
Thirsty 0.388 0.121 0.078 0.036 0.053
Temp 0.314 0     0.206 0.158 0.08 
Heart 0.207 0.043 0.169 0.053 0.027
Urinate 0.214 0.574 0.028 0.098 0.016
Defecate 0.121 0.574 0.106 0.041 0.077
Bruise 0.055 0     0.45  0     0.094
Pain 0.154 0.051 0.444 0.049 0.062
Itch 0.054 0     0.353 0.098 0.12 
Affective_touch 0.113 0     0.31  0     0.128
Blood_Sugar 0     0     0.303 0.037 0.049
Muscles 0.178 0.045 0.268 0.047 0.09 
Taste 0.151 0.061 0.189 0.067 0.101
Sneeze 0.048 0.036 0.125 0.621 0.043
Cough 0.107 0.043 0.02  0.502 0.103
Vomit 0.095 0.064 0.116 0.228 0.126
Wind 0.017 0     0.037 0.09  0.651
Burp 0.034 0     0.27  0.113 0.465
Sex_arousal 0.154 0.093 0.24  0.061 0.202
Code
plot(rez_efa$n_fac)

Final

Note: CFA did not converge

Code
rez_hclust <- make_hclust(data1bf)
rez_ega <- make_ega(data1af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data1bf, n = 3, cols = c("F1" = colors[["Urinate"]], "F2" = colors[["Bruise"]], 
                                            "F3" = colors[["Cough"]], "F4" = colors[["Itch"]]))
# rez_cfa1b  <- make_cfa(data1bf)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.511 0.128  0.036 0.009 0.01  0.108  0    
Breathing 0.511 0.16   0.048 0.068 0.015 0.151  0.045
Hungry 0.202 0.514  0.105 0.015 0     0.042  0.048
Thirsty 0.091 0.514  0.178 0     0.018 0.041  0.089
Urinate 0.029 0.253  0.565 0.038 0.05  0.185 −0.017
Defecate 0.08  0.108  0.565 0.098 0.087 0.09  −0.042
Sneeze 0.033 0.02   0.12  0.574 0.117 0.091  0    
Cough 0.071 0      0.021 0.574 0.171 0.06   0.157
Wind 0.021 0.024  0.136 0.071 0.577 0.028  0.089
Burp 0.012 0      0.008 0.22  0.577 0.041  0.089
Muscles 0.091 0.042  0.118 0.058 0.04  0.523  0.125
Pain 0.183 0.044  0.107 0.061 0.014 0.523  0.18 
Bruise 0.01  0.128 −0.031 0.059 0.093 0.222  0.502
Blood_Sugar 0.033 0     −0.013 0.052 0.032 0.053  0.502
Code
plot(rez_efa$n_fac)

Code
# rez_cfa1b$test

Initial

Code
rez_hclust <- make_hclust(data2)
rez_ega <- make_ega(data2, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Urinate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Pain"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data2, n = 4, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Bruise"]], 
                                            "F3" = colors[["Heart"]], "F4" = colors[["Sex arousal"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.539 0.057 0.053 0     0.051 0.067  0.047
Breathing 0.539 0.126 0.067 0.027 0.16  0      0.058
Hungry 0.096 0.524 0.066 0.028 0.18  0.026  0.069
Thirsty 0.074 0.524 0.166 0.035 0.027 0.091  0    
Defecate 0.004 0.057 0.52  0.097 0.061 0.126  0    
Urinate 0     0.175 0.511 0.01  0.065 0.021 −0.013
Sex_arousal 0.04  0.089 0.306 0.04  0.135 0.051  0.086
Vomit 0.027 0.002 0.275 0.171 0.024 0.037  0.062
Taste 0.084 0     0.199 0.074 0.148 0.086  0.056
Sneeze 0.001 0.031 0.233 0.542 0.033 0.138  0.006
Cough 0.026 0.039 0.074 0.542 0.238 0.198  0.103
Pain 0.078 0.086 0.027 0.038 0.47  0.144  0.159
Muscles 0.061 0.027 0.072 0.074 0.337 0.076  0.09 
Temp 0.058 0.095 0.213 0.138 0.299 0.055  0.039
Wind 0     0.026 0.162 0.156 0.142 0.51   0.044
Burp 0.058 0.082 0.051 0.13  0.113 0.51   0.125
Blood_Sugar 0.039 0     0.025 0     0     0      0.469
Bruise 0.009 0.03  0.006 0.023 0.229 0.086  0.375
Itch 0.038 0     0.042 0.1   0.073 0.091  0.332
Affective_touch 0.034 0.055 0.095 0     0.052 0.048  0.259
Code
plot(rez_efa$n_fac)

Final

Note: CFA did not converge

Code
rez_hclust <- make_hclust(data2f)
rez_ega <- make_ega(data2f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data2f, n = 3, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Bruise"]], 
                                            "F3" = colors[["Hungry"]]))
# rez_cfa2  <- make_cfa(data2f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.544 0.062 0.021 0     0.083 0.078 0.037
Breathing 0.544 0.148 0.019 0.066 0     0.165 0.06 
Hungry 0.113 0.528 0.077 0.041 0.043 0.154 0.047
Thirsty 0.081 0.528 0.16  0.039 0.096 0.041 0    
Urinate 0.013 0.21  0.567 0.058 0.06  0.031 0    
Defecate 0.031 0.075 0.567 0.136 0.157 0.097 0    
Sneeze 0.021 0.048 0.139 0.56  0.161 0.024 0    
Cough 0.05  0.045 0.048 0.56  0.233 0.215 0.085
Wind 0.009 0.047 0.113 0.175 0.514 0.159 0.021
Burp 0.063 0.083 0.055 0.141 0.514 0.141 0.107
Muscles 0.077 0.044 0.046 0.087 0.083 0.479 0.097
Pain 0.097 0.108 0.037 0.073 0.166 0.479 0.18 
Bruise 0.017 0.037 0     0.039 0.112 0.314 0.508
Blood_Sugar 0.064 0.005 0     0.027 0.013 0.01  0.508
Code
plot(rez_efa$n_fac)

Code
# rez_cfa2$test

Initial

Code
rez_hclust <- make_hclust(data3)
rez_ega <- make_ega(data3, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Urinate"]], "C3" = colors[["Pain"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Muscles"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data3, n = 4, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Urinate"]], 
                                            "F3" = colors[["Blood Sugar"]], "F4" = colors[["Sex arousal"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing 0.618 0.154 0.01  0.005 0     0.046
Heart 0.349 0.002 0.146 0.043 0.024 0.128
Sex_arousal 0.27  0.122 0.127 0.06  0.177 0.022
Urinate 0.078 0.51  0.127 0.093 0     0    
Thirsty 0     0.454 0.055 0.041 0     0.023
Defecate 0.121 0.425 0.063 0.108 0.007 0.048
Hungry 0.158 0.364 0.076 0.052 0.025 0.019
Muscles 0.122 0.025 0.506 0.082 0.031 0.167
Pain 0.022 0.153 0.356 0.091 0.036 0.151
Taste 0.095 0.033 0.221 0.129 0.028 0.064
Cough 0     0     0.061 0.65  0.118 0.047
Sneeze 0.061 0.049 0.078 0.589 0.066 0.053
Temp 0.019 0.169 0.216 0.246 0.085 0    
Vomit 0.055 0.071 0.095 0.23  0.152 0    
Wind 0.168 0.022 0.041 0.118 0.562 0.041
Burp 0.04  0.003 0.074 0.227 0.562 0.076
Blood_Sugar 0.108 0.008 0.019 0     0     0.48 
Itch 0.053 0     0.068 0.031 0.112 0.395
Bruise 0     0     0.257 0.026 0.019 0.39 
Affective_touch 0.065 0.074 0.178 0.034 0     0.295
Code
plot(rez_efa$n_fac)

Final

Note: CFA did not converge

Code
rez_hclust <- make_hclust(data3f)
rez_ega <- make_ega(data3f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Urinate"]], "C3" = colors[["Sneeze"]], "C4" = colors[["Wind"]], 
  "C5" = colors[["Pain"]], "C6" = colors[["Bruise"]], "C7" = colors[["Muscles"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data3f, n = 3, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Hungry"]], 
                                            "F3" = colors[["Blood Sugar"]]))
# rez_cfa3  <- make_cfa(data3f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.537 0.006 0.058 0.031 0.086 0.134
Breathing 0.537 0.199 0     0.031 0.066 0.015
Urinate 0.1   0.527 0.036 0.006 0.141 0    
Thirsty 0     0.468 0     0     0.084 0.029
Defecate 0.035 0.44  0.067 0.067 0.096 0    
Hungry 0.16  0.378 0.026 0.049 0.086 0    
Sneeze 0.085 0.095 0.623 0.111 0.054 0    
Cough 0     0.037 0.623 0.16  0.099 0.061
Wind 0.016 0.078 0.099 0.581 0.059 0    
Burp 0.06  0.026 0.127 0.581 0.118 0.046
Muscles 0.122 0.064 0.054 0.075 0.515 0.141
Pain 0.014 0.19  0.04  0.055 0.515 0.125
Bruise 0     0     0.04  0.039 0.248 0.544
Blood_Sugar 0.154 0.021 0.003 0     0.058 0.544
Code
plot(rez_efa$n_fac)

Code
# rez_cfa3$test

Initial

Code
rez_hclust <- make_hclust(data4)
rez_ega <- make_ega(data4, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Pain"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Affective touch"]]))
rez_efa <- make_efa(data4, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Itch"]], 
                                            "F3" = colors[["Breathing"]], "F4" = colors[["Urinate"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7 8
Heart 0.617 0      0.021 0.044 0     0.01   0.151 0.043
Breathing 0.617 0.154  0.045 0     0.067 0.028  0.029 0.026
Hungry 0.07  0.56   0.151 0     0.094 0.022  0     0.007
Thirsty 0.049 0.56   0.106 0.009 0.051 0.014  0     0.023
Urinate 0     0.161  0.559 0.041 0.108 0.04  −0.01  0    
Defecate 0.016 0.02   0.524 0.08  0.072 0.093  0     0.034
Taste 0.035 0.072  0.169 0.08  0.131 0.089  0     0.08 
Sneeze 0.003 0.01   0.036 0.587 0.016 0.117  0     0.07 
Cough 0.005 0      0.1   0.462 0.063 0.203  0.065 0.04 
Vomit 0.029 0      0.083 0.291 0.122 0.069  0.097 0.005
Pain 0.036 0.032  0.156 0.086 0.527 0.028  0.14  0.032
Muscles 0.002 0.023  0.069 0.071 0.449 0.136  0.06  0.048
Temp 0.019 0.102  0.117 0.045 0.372 0.148  0     0.051
Wind 0.006 0      0.073 0.107 0.085 0.604  0.039 0.024
Burp 0     0      0     0.178 0.093 0.451  0.02  0.02 
Sex_arousal 0.022 0.035  0.146 0.067 0.102 0.186  0     0.085
Bruise 0.05  0      0     0.084 0.085 0.006  0.448 0.059
Blood_Sugar 0.028 0     −0.006 0     0.018 0.027  0.448 0.143
Tickle 0.014 0      0.003 0.056 0.014 0.055  0.009 0.62 
Itch 0.026 0      0.04  0.042 0.053 0.01   0.257 0.542
Affective_touch 0.023 0.035  0.093 0.027 0.076 0.091  0.159 0.275
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data4f)
rez_ega <- make_ega(data4f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data4f, n = 3, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Breathing"]], 
                                            "F3" = colors[["Urinate"]]))
rez_cf4  <- make_cfa(data4f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.619 0     0.031 0.028 0.004 0     0.197
Breathing 0.619 0.173 0     0.003 0.027 0.086 0.038
Hungry 0.083 0.565 0.126 0     0     0.069 0.006
Thirsty 0.053 0.565 0.09  0.022 0     0.069 0.003
Urinate 0     0.185 0.582 0.061 0.017 0.112 0    
Defecate 0.026 0.048 0.582 0.076 0.065 0.074 0.025
Sneeze 0.012 0.023 0.054 0.57  0.09  0.052 0.111
Cough 0.013 0     0.076 0.57  0.189 0.093 0.122
Wind 0.022 0     0.077 0.129 0.569 0.088 0.094
Burp 0.004 0     0     0.149 0.569 0.11  0.067
Muscles 0.011 0.058 0.041 0.068 0.131 0.556 0.152
Pain 0.054 0.074 0.124 0.067 0.055 0.556 0.137
Blood_Sugar 0.036 0     0     0     0.04  0.038 0.434
Bruise 0.058 0     0     0.051 0.011 0.121 0.351
Tickle 0.045 0.007 0.017 0.12  0.068 0.067 0.238
Code
plot(rez_efa$n_fac)

Code
rez_cf4$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 24092.56 0.131 802.775 0.734 Inf
fit_ega 28320.01 0.025 83.537 0.991 NA
fit_egag 28365.84 0.052 222.983 0.950 1.118915e-10
fit_hclust 28472.07 0.067 309.152 0.920 9.577515e-34
fit_hclustg 28503.33 0.072 373.847 0.901 1.560998e-40
fit_efag 28732.35 0.096 616.239 0.822 2.902305e-90
fit_efa 28732.35 0.096 616.239 0.822 2.902305e-90

Initial

Code
rez_hclust <- make_hclust(data5)
rez_ega <- make_ega(data5, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Sneeze"]], "C3" = colors[["Vomit"]], "C4" = colors[["Cough"]], 
  "C5" = colors[["Defecate"]], "C6" = colors[["Muscles"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data5, n = 8, cols = c("F1" = colors[["Urinate"]], "F2" = colors[["Itch"]], 
                                            "F3" = colors[["Burp"]], "F4" = colors[["Hungry"]], 
                                           "F5" = colors[["Affective touch"]], "F6" = colors[["Breathing"]], "F7" = colors[["Sex arousal"]], "F8" = colors[["Taste"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.435 0.01 0.075 0.039 0.089 0.164  0.046
Thirsty 0.435 0.154 0     0.094 0.084 0.16   0.037
Burp 0     0.454 0     0.071 0.082 0.117  0.042
Sneeze 0.084 0.323 0.069 0.131 0.024 0.052  0.112
Hungry 0.153 0.226 0.054 0.01  0.062 0.085 −0.051
Breathing 0.057 0 0.517 0     0.087 0.115  0.049
Vomit 0.06  0.134 0.517 0.13  0.018 0.183  0    
Urinate 0.089 0.014 0.129 0.534 0.143 0.006  0    
Cough 0.114 0.202 0     0.404 0.059 0.252  0.085
Pain 0.074 0.09 0.044 0.404 0.023 0.246  0.055
Taste 0.152 0.026 0.094 0     0.343 0.266  0    
Affective_touch 0     0.073 0     0.017 0.337 0.074  0.005
Defecate 0.086 0.068 0.012 0.154 0.3   0.069  0.015
Blood_Sugar 0.036 0.017 0     0     0.133 0.011  0    
Wind 0.21  2.563 × 10−4 0.19  0.128 0.172 0.368  0.006
Sex_arousal 0.201 0.146 0     0.089 0.1   0.346  0    
Muscles 0     0.059 0.082 0.056 0.131 0.293  0.096
Temp 0.162 0.105 0.065 0.155 0.065 0.245  0    
Bruise 0.153 0.118 0     0.068 0.024 0.006  0.55 
Itch 0     0.012 0.057 0.055 0     0.1    0.55 
Code
plot(rez_efa$n_fac)

Final

Note: EGA errors. Note: CFA did not converge

Code
rez_hclust <- make_hclust(data5f)
# rez_ega <- make_ega(data5f, cols = c(
#   "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
#   "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data5f, n = 5, cols = c("F1" = colors[["Urinate"]], "F2" = colors[["Burp"]], 
                                            "F3" = colors[["Wind"]], "F4" = colors[["Thirsty"]],
                                            "F5" = colors[["Sneeze"]]))
# rez_cfa5  <- make_cfa(data5f)

rez_hclust$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)

Code
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

plot(rez_efa$n_fac)

Code
# rez_cfa5$test

Initial

Code
rez_hclust <- make_hclust(data6)
rez_ega <- make_ega(data6, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Itch"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data6, n = 3, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Hungry"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.524 0.034 0     0.086 0     7.366 × 10−4 0.172
Breathing 0.524 0.135 0     0.127 0.028 0.078 0.116
Hungry 0.124 0.554 0.023 0.124 0     0.086 0.028
Thirsty 0.072 0.554 0.095 0.045 0     0.045 0.007
Urinate 0 0.082 0.612 0.037 0.052 0.015 0    
Defecate 0 0.071 0.612 0.096 0.053 0.205 0    
Wind 0.025 0.014 0.059 0.545 0.057 0.044 0.043
Burp 0.105 0     0     0.537 0.13  0.004 0.178
Sex_arousal 0 0.099 0.031 0.288 0.048 0.185 0    
Muscles 0.126 0.033 0     0.259 0.052 0.216 0.116
Taste 0.044 0.06  0.036 0.2   0.075 0.106 0    
Cough 0.025 0     0     0.154 0.602 0.122 0.061
Sneeze 0 0     0.099 0.051 0.588 0.152 0    
Vomit 0.014 0     0     0.156 0.253 0.099 0    
Itch 0 0     0.06  0.102 0.087 0.509 0    
Pain 9.275 × 10−4 0.077 0.015 0.1   0.048 0.431 0.24 
Affective_touch 0 0.065 0.013 0.161 0.03  0.275 0.085
Temp 0.098 0     0.098 0.132 0.168 0.243 0    
Bruise 0.135 0     0     0.087 0.016 0.147 0.46 
Blood_Sugar 0.071 0.021 0     0.083 0.015 0.037 0.46 
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data6f)
rez_ega <- make_ega(data6f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Pain"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data6f, n = 3, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Thirsty"]]))
rez_cfa6  <- make_cfa(data6f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.525 0.036 0     0     0.022 0.174
Breathing 0.525 0.148 0.007 0.038 0.079 0.14 
Hungry 0.141 0.559 0.039 0     0.039 0.116
Thirsty 0.076 0.559 0.1   0     0     0.079
Urinate 0     0.098 0.618 0.055 0     0.006
Defecate 0.011 0.082 0.618 0.086 0.085 0.104
Sneeze 0.005 0     0.138 0.627 0.013 0.101
Cough 0.055 0     0.008 0.627 0.148 0.146
Wind 0.036 0.047 0.08  0.02  0.604 0.11 
Burp 0.109 0     0     0.126 0.604 0.225
Bruise 0.136 0     0     0.013 0.064 0.437
Pain 0.013 0.118 0.066 0.073 0.015 0.399
Muscles 0.152 0.047 0.018 0.08  0.145 0.284
Blood_Sugar 0.071 0.031 0     0.017 0.051 0.249
Code
plot(rez_efa$n_fac)

Code
rez_cfa6$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 17201.04 0.144 745.632 0.776 Inf
fit_egag 20221.10 0.058 217.353 0.957 5.020192e+00
fit_ega 20224.33 0.046 130.586 0.978 NA
fit_efag 20516.73 0.101 538.689 0.865 3.212108e-64
fit_efa 20516.73 0.101 538.689 0.865 3.212108e-64
fit_hclust 20579.19 0.109 556.159 0.857 8.756386e-78
fit_hclustg 20587.13 0.108 596.235 0.847 1.656411e-79

Initial

Code
rez_hclust <- make_hclust(data7a)
rez_ega <- make_ega(data7a, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Defecate"]], "C3" = colors[["Itch"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Pain"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data7a, n = 4, cols = c("F1" = colors[["Pain"]], "F2" = colors[["Sneeze"]], 
                                            "F3" = colors[["Burp"]], "F4" = colors[["Urinate"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing  0.493 0     7.190 × 10−4  0.083 0     0.146
Thirsty  0.349 0.119 0.053  0     0.018 0.061
Hungry  0.201 0.035 0.105 −0.009 0.098 0.131
Heart  0.193 0     0  0     0.047 0.145
Urinate  0.095 0.572 0.119  0     0.078 0.049
Defecate  0.085 0.572 0.066  0.097 0.067 0.098
Taste  0.057 0.04  0.468  0.095 0.095 0.082
Temp  0.053 0.07  0.468  0.07  0.026 0.231
Sneeze −0.007 0     0  0.639 0.002 0.089
Cough −0.005 0     0.247  0.597 0.152 0.052
Vomit  0.111 0.11  0.068  0.223 0.18  0.021
Wind  0.043 0.037 0.048  0.074 0.656 0.056
Burp  0.012 0.003 0.077  0.169 0.463 0.158
Sex_arousal  0.134 0.104 0.076  0.048 0.193 0.126
Muscles  0.113 0.061 0.094  0.008 0.051 0.411
Pain  0.229 0.021 0.093  0.032 0.028 0.385
Itch  0.092 0.07  0.288  0.085 0.033 0.384
Bruise  0.054 0.006 0  0     0.13  0.352
Affective_touch  0.054 0.002 0.137  0.003 0.085 0.346
Blood_Sugar  0.115 0.011 0  0.039 0.072 0.232
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data7af)
rez_ega <- make_ega(data7af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Pain"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Burp"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data7af, n = 4, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Pain"]], 
                                            "F3" = colors[["Sneeze"]], "F4" = colors[["Urinate"]]))
rez_cf7a  <- make_cfa(data7af)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing 0.546 0.183 0     0.006 0.012 0.154
Heart 0.302 0.035 0     0     0.057 0.158
Thirsty 0.243 0.155 0.135 0.013 5.718 × 10−4 0.045
Pain 0.179 0.487 0.07  0.063 0 0.168
Muscles 0.065 0.441 0.079 0.023 0.041 0.321
Hungry 0.196 0.315 0.069 0     0 0    
Urinate 0.087 0.082 0.585 0.015 0.068 0    
Defecate 0.094 0.166 0.585 0.029 0.03 0.063
Sneeze 0     0.113 0.019 0.636 0.022 0    
Cough 0.032 0.009 0.036 0.636 0.176 0.099
Wind 0.071 0.049 0.071 0.038 0.594 0.119
Burp 0.027 0     0.032 0.128 0.594 0.325
Bruise 0.085 0.204 0.014 0     0.134 0.441
Blood_Sugar 0.146 0.064 0.017 0.038 0.071 0.441
Code
plot(rez_efa$n_fac)

Code
rez_cf7a$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 14844.13 0.129 522.733 0.646 Inf
fit_ega 17647.14 0.038 98.571 0.973 NA
fit_egag 17657.35 0.059 196.386 0.920 6.073383e-03
fit_efa 17782.60 0.084 346.672 0.827 3.842713e-30
fit_efag 17782.60 0.084 346.672 0.827 3.842713e-30
fit_hclust 17812.14 0.087 332.401 0.831 1.484630e-36
fit_hclustg 17819.37 0.089 370.921 0.810 3.993394e-38

Initial

Code
rez_hclust <- make_hclust(data7b)
rez_ega <- make_ega(data7b, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Pain"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data7b, n = 4, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Urinate"]], 
                                            "F3" = colors[["Heart"]], "F4" = colors[["Sneeze"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.56  0     0.026 0.019 0.052 0.03  0.208
Breathing 0.56  0.194 0.02  0.041 0.127 0.01  0.008
Hungry 0.047 0.477 0.07  0.007 0.089 0     0.091
Thirsty 0.064 0.376 0.124 0.013 0.013 0.023 0.046
Sex_arousal 0.044 0.203 0.09  0.092 0.174 0.128 0    
Defecate 0.009 0.122 0.601 0.11 0.067 0     0    
Urinate 0.013 0.184 0.537 5.716 × 10−4 0.133 0     0    
Taste 0.028 0.077 0.201 0.038 0.199 0.071 0.14 
Sneeze 0     0.01  0.037 0.579 0.143 0.08  0    
Cough 0.033 0.054 0.045 0.474 0.122 0.163 0    
Vomit 0.033 0.089 0.067 0.296 0.045 0.078 0.007
Itch 0.033 0.016 0.113 0.123 0.505 0.003 0    
Pain 0.054 0.102 0.077 0.045 0.399 0.035 0.194
Muscles 0.064 0.075 0.059 0.035 0.32  0.105 0.146
Temp 0.025 0.102 0.118 0.11 0.318 0.036 0    
Affective_touch 0.031 0.107 0.064 0.018 0.244 0.05  0.106
Wind 0.011 0.097 0.027 0.081 0.116 0.548 0.203
Burp 0.026 0.081 0.036 0.199 0.069 0.548 0.167
Bruise 0.011 0.062 0.008 0 0.188 0.172 0.449
Blood_Sugar 0.111 0.036 0.066 0.004 0.029 0.049 0.449
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data7bf)
rez_ega <- make_ega(data7bf, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data7bf, n = 3, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Urinate"]], 
                                            "F3" = colors[["Breathing"]]))
rez_cfa7b  <- make_cfa(data7bf)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.616 0.002 0.012 0.018 0.04  0.043
Breathing 0.44  0.168 0.032 0.039 0.018 0.156
Blood_Sugar 0.187 0.054 0.011 0 0.071 0.153
Hungry 0.097 0.495 0.059 0 0.021 0.162
Thirsty 0.069 0.495 0.107 0.02 0.033 0.041
Urinate 0.023 0.179 0.604 0.027 0.017 0.108
Defecate 0.045 0.098 0.604 0.102 0.029 0.073
Sneeze 0.01  0     0.059 0.578 0.112 0.128
Cough 0.054 0.03  0.056 0.578 0.185 0.093
Wind 0.081 0.018 0.031 0.076 0.557 0.189
Burp 0.048 0.056 0.006 0.195 0.557 0.206
Pain 0.089 0.147 0.068 0.11 0.058 0.437
Muscles 0.088 0.036 0.064 0.074 0.138 0.387
Bruise 0.149 0.072 0.003 7.260 × 10−4 0.17  0.347
Code
plot(rez_efa$n_fac)

Code
rez_cfa7b$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 57540.78 0.106 1,258.337 0.789 Inf
fit_ega 68203.78 0.033 178.485 0.982 NA
fit_egag 68259.45 0.044 340.521 0.959 8.144885e-13
fit_efa 68668.80 0.069 780.259 0.893 1.051831e-101
fit_efag 68668.80 0.069 780.259 0.893 1.051831e-101
fit_hclustg 68700.47 0.071 796.735 0.891 1.395370e-108
fit_hclust 68709.50 0.072 767.779 0.894 1.526592e-110

Initial

Code
rez_hclust <- make_hclust(data7c)
rez_ega <- make_ega(data7c, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Pain"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Sex arousal"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data7c, n = 3, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Urinate"]], 
                                            "F3" = colors[["Bruise"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7 8
Heart 0.521 0     0.045 0.039 0.009 0     0.053 0.117
Breathing 0.521 0.159 0.059 0.201 0.036 0.007 0     0.025
Hungry 0.074 0.54  0.09  0.09  0     0.124 0     0    
Thirsty 0.101 0.54  0.082 0.04  0     0.056 0.046 0.022
Urinate 0.111 0.125 0.555 0.076 0.015 0.065 0.026 0    
Defecate 0.011 0.061 0.555 0.127 0.087 0.122 0.067 0    
Pain 0.068 0.078 0.065 0.513 0     0.144 0.015 0.095
Temp 0     0.034 0.061 0.353 0.09  0.138 0.046 0.029
Taste 0.101 0.032 0.083 0.281 0.063 0.074 0     0    
Muscles 0.125 0     0     0.241 0.001 0.102 0.055 0.146
Sneeze 0.045 0     0.017 0.011 0.648 0.034 0.121 0    
Cough 0.014 0     0     0.15  0.484 0     0.136 0.138
Vomit 0.005 0     0.107 0.021 0.318 0.007 0.04  0.009
Sex_arousal 0     0.116 0.138 0.142 0.004 0.494 0.097 0    
Affective_touch 0.006 0.027 0     0.186 0.02  0.494 0     0.124
Wind 0.069 0.029 0.069 0.062 0.147 0.052 0.577 0.001
Burp 0     0.026 0.034 0.062 0.123 0.093 0.577 0.138
Bruise 0.048 0.022 0     0.211 0.031 0     0.069 0.455
Blood_Sugar 0.1   0     0     0.031 0.005 0     0     0.374
Itch 0.011 0     0     0.003 0.076 0.157 0.049 0.351
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data7cf)
rez_ega <- make_ega(data7cf, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data7cf, n = 4, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Hungry"]], 
                                            "F3" = colors[["Bruise"]], "F4" = colors[["Breathing"]]))
rez_cfa7c  <- make_cfa(data7cf)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.522 0     0.047 0.01  0.054 0.032 0.15 
Breathing 0.522 0.164 0.066 0.047 0     0.233 0.038
Hungry 0.081 0.544 0.114 0     0.013 0.116 0    
Thirsty 0.102 0.544 0.085 0     0.058 0.037 0.041
Urinate 0.124 0.142 0.568 0.011 0.038 0.021 0    
Defecate 0.015 0.08  0.568 0.056 0.088 0.146 0    
Sneeze 0.054 0     0.053 0.603 0.137 0.004 0    
Cough 0.03  0     0.024 0.603 0.16  0.048 0.098
Wind 0.071 0.039 0.077 0.155 0.581 0.073 0    
Burp 0     0.045 0.058 0.115 0.581 0.097 0.127
Muscles 0.133 0     0.007 0     0.076 0.489 0.206
Pain 0.091 0.116 0.107 0.03  0.033 0.489 0.148
Bruise 0.05  0.028 0     0.056 0.08  0.298 0.487
Blood_Sugar 0.107 0.003 0     0     0     0.051 0.487
Code
plot(rez_efa$n_fac)

Code
rez_cfa7c$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 19927.16 0.115 514.981 0.728 Inf
fit_ega 23336.33 0.036 102.630 0.976 NA
fit_egag 23351.38 0.055 208.177 0.928 5.413702e-04
fit_hclustg 23444.76 0.072 314.491 0.873 2.853466e-24
fit_hclust 23462.50 0.074 299.902 0.878 4.020609e-28
fit_efa 23478.41 0.078 361.067 0.850 1.409458e-31
fit_efag 23478.41 0.078 361.067 0.850 1.409458e-31

Initial

Code
rez_hclust <- make_hclust(data8a)
rez_ega <- make_ega(data8a, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Itch"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data8a, n = 5, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Thirsty"]], 
                                            "F3" = colors[["Bruise"]], "F4" = colors[["Heart"]], "F5" = colors[["Urinate"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing 0.472 0.137 0.019 0.066 0.037 0.064
Heart 0.467 0.061 0     0.054 0     0.082
Sex_arousal 0.238 0.134 0.064 0.044 0.082 0.097
Hungry 0.155 0.454 0.095 0.027 0.049 0.016
Thirsty 0.01  0.405 0.118 0.041 0     0.102
Temp 0.12  0.337 0.028 0.067 0.008 0.137
Taste 0.072 0.158 0.025 0.113 0.09  0.137
Urinate 0.076 0.115 0.544 0.063 0     0.103
Defecate 0.008 0.136 0.544 0.13  0.08  0.066
Sneeze 0.03  0.084 0.068 0.545 0.208 0.07 
Cough 0.053 0.147 0.052 0.4   0.285 0.013
Vomit 0.098 0.019 0.088 0.316 0.057 0.094
Wind 0.059 0.094 0.011 0.184 0.503 0.08 
Burp 0.039 0.02  0.054 0.233 0.503 0.148
Bruise 0.075 0.004 0.01  0.029 0.041 0.479
Itch 0.107 0.081 0.007 0.05  0.046 0.396
Affective_touch 0.018 0.037 0     0.036 0.109 0.373
Muscles 0.031 0.109 0.009 0.007 0.095 0.319
Blood_Sugar 0.049 0.009 0.083 0.091 0.049 0.304
Pain 0.051 0.256 0.117 0.009 0.035 0.278
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data8af)
rez_ega <- make_ega(data8af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Pain"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data8af, n = 4, cols = c("F1" = colors[["Cough"]], "F2" = colors[["Thirsty"]], 
                                            "F3" = colors[["Bruise"]], "F4" = colors[["Heart"]]))
rez_cfa8a  <- make_cfa(data8af)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.531 0.072 0     0.045 0     0.128
Breathing 0.531 0.086 0.041 0.095 0.06  0.051
Hungry 0.146 0.533 0.107 0.052 0.059 0.044
Thirsty 0.014 0.533 0.122 0.063 0     0.138
Urinate 0.031 0.09  0.546 0.065 0.003 0.152
Defecate 0.013 0.154 0.546 0.102 0.089 0.09 
Sneeze 0.063 0.032 0.096 0.533 0.238 0.099
Cough 0.078 0.082 0.06  0.533 0.301 0.02 
Wind 0.054 0.053 0.018 0.199 0.509 0.095
Burp 0     0     0.059 0.28  0.509 0.162
Bruise 0.031 0     0.036 0     0.094 0.464
Muscles 0.043 0.013 0.015 0.026 0.121 0.401
Blood_Sugar 0.062 0     0.092 0.087 0.066 0.316
Pain 0.082 0.207 0.129 0.031 0.067 0.28 
Code
plot(rez_efa$n_fac)

Code
rez_cfa8a$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 34374.91 0.091 570.623 0.795 Inf
fit_ega 40816.35 0.033 128.727 0.976 NA
fit_hclust 40867.57 0.049 257.618 0.936 7.558221e-12
fit_egag 40871.15 0.051 282.380 0.929 1.263321e-12
fit_hclustg 40889.24 0.054 314.597 0.919 1.487098e-16
fit_efag 41018.42 0.067 457.898 0.872 1.322953e-44
fit_efa 41018.42 0.067 457.898 0.872 1.322953e-44

Initial

Code
rez_hclust <- make_hclust(data8b)
rez_ega <- make_ega(data8b, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Cough"]], 
  "C5" = colors[["Pain"]], "C6" = colors[["Bruise"]], "C7" = colors[["Itch"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data8b, n = 5, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Thirsty"]], 
                                            "F3" = colors[["Breathing"]], "F4" = colors[["Bruise"]], "F5" = colors[["Urinate"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing 0.575  0.032 0 0.077 0      0.025
Heart 0.422  0.037 0 0.056 0.061  0.024
Sex_arousal 0.233  0.165 0.021 0.083 0      0.12 
Hungry 0.124  0.469 0.039 0.099 0.058  0.02 
Thirsty 0.013  0.392 0.124 0.023 0.15  −0.066
Temp 0.08   0.377 2.195 × 10−4 0.05  0.093  0.033
Taste 0.038  0.22  0 0.113 0.188  0.083
Urinate 0.022  0.088 0.563 0.059 0.015  0.115
Defecate 0      0.07  0.563 0.197 0.125  0.009
Cough 0.062  0.124 0.06 0.465 0      0.03 
Burp 0.035  0.052 0.092 0.444 0.018  0.124
Sneeze 0.07   0.081 0.09 0.404 0.052  0.091
Wind 0.051  0.02  0 0.38  0.155  0.018
Vomit 0.055  0.052 0.065 0.21  0.029  0.139
Muscles 0.025  0.091 0.017 0.06  0.468  0.192
Pain 0.015  0.199 0.07 0.07  0.468  0.07 
Affective_touch 0      0.023 0.011 0.084 0.083  0.408
Bruise 0.025 −0.006 0.027 0.041 0.17   0.4  
Blood_Sugar 0.026 −0.008 0.095 0.089 0.066  0.352
Itch 0.137  0.063 0 0.143 0.134  0.338
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data8bf)
rez_ega <- make_ega(data8bf, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Cough"]], 
  "C5" = colors[["Pain"]], "C6" = colors[["Bruise"]], "C7" = colors[["Itch"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data8bf, n = 4, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Heart"]], 
                                            "F3" = colors[["Thirsty"]], "F4" = colors[["Bruise"]]))
rez_cfa8b  <- make_cfa(data8bf)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.551  0.054 0     0.053 0.078  0.051
Breathing 0.551  0.046 0     0.09  0.012  0    
Hungry 0.075  0.538 0.045 0.115 0.101  0.058
Thirsty 0.019  0.538 0.13  0.035 0.161 −0.107
Urinate 0      0.111 0.566 0.051 0.019  0.189
Defecate 0      0.088 0.566 0.173 0.127  0.022
Burp 0      0.036 0.094 0.499 0.05   0.116
Cough 0.083  0.097 0.067 0.473 0      0.051
Wind 0      0.014 0     0.409 0.18   0    
Sneeze 0.102  0.06  0.11  0.359 0.088  0.017
Muscles 0.042  0.013 0.024 0.083 0.48   0.243
Pain 0.021  0.182 0.071 0.088 0.48   0.061
Bruise 0     −0.034 0.044 0.031 0.241  0.499
Blood_Sugar 0.039 −0.006 0.108 0.079 0.095  0.499
Code
plot(rez_efa$n_fac)

Code
rez_cfa8b$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 15186.05 0.097 308.676 0.751 Inf
fit_egag 18012.76 0.048 151.514 0.932 1.159343e+04
fit_ega 18031.48 0.031 83.225 0.977 NA
fit_hclustg 18052.96 0.061 204.139 0.890 2.166613e-05
fit_hclust 18058.18 0.058 178.287 0.907 1.592270e-06
fit_efa 18123.87 0.076 287.481 0.822 8.654529e-21
fit_efag 18123.87 0.076 287.481 0.822 8.654529e-21

Initial

Code
rez_hclust <- make_hclust(data9)
rez_ega <- make_ega(data9, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data9, n = 3, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Thirsty"]], 
                                            "F3" = colors[["Bruise"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.556 0.036 0.007 0.051 0.037 0.007 0.102
Breathing 0.556 0.237 0.041 0.161 0.02  0     0    
Hungry 0.15 0.536 0.068 0.097 0     0     0.036
Thirsty 0.098 0.536 0.128 0.102 0     0     0    
Urinate 0.015 0.223 0.607 0.114 0.018 0.011 0    
Defecate 0.045 0.047 0.607 0.129 0.098 0.04  0    
Pain 0.026 0.087 0.037 0.426 0.032 0     0.164
Temp 0.012 0.073 0.078 0.405 0.13  0.018 0    
Muscles 0.09 0.03  0.018 0.384 0.042 0.089 0.123
Sex_arousal 0.054 0.047 0.048 0.316 0.065 0.064 0    
Affective_touch 0.018 0.015 0.01 0.285 0.034 0     0.2  
Taste 0.06 0.016 0.047 0.219 0.112 0.016 0.059
Sneeze 0.03 0     0.028 0.081 0.643 0.063 0.04 
Cough 0.024 0     7.041 × 10−4 0.14  0.466 0.133 0.062
Vomit 0.014 0     0.081 0.182 0.316 0.032 0.017
Wind 2.528 × 10−4 0     0.052 0.097 0.084 0.609 0.036
Burp 0.008 0     0 0.096 0.159 0.609 0.076
Bruise 0.01 0     0 0.183 0.023 0.05  0.444
Blood_Sugar 0.055 0.037 0 0.077 0.007 0.006 0.391
Itch 0.03 0     0 0.156 0.064 0.026 0.342
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data9f)
rez_ega <- make_ega(data9f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data9f, n = 3, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Urinate"]], 
                                            "F3" = colors[["Pain"]]))
rez_cfa9  <- make_cfa(data9f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.559 0.039 0.012 0.041 0.015 0.092
Breathing 0.559 0.244 0.051 0.033 0     0.104
Hungry 0.158 0.539 0.078 0.003 0     0.087
Thirsty 0.1   0.539 0.133 0.007 0     0.075
Urinate 0.02  0.232 0.613 0.03  0.02  0.06 
Defecate 0.059 0.063 0.613 0.067 0.052 0.061
Sneeze 0.049 0.004 0.07  0.598 0.085 0.059
Cough 0.038 0.009 0.02  0.598 0.145 0.118
Wind 0.005 0     0.07  0.098 0.612 0.063
Burp 0.014 0     0.002 0.146 0.612 0.135
Bruise 0.017 0     0     0.036 0.059 0.477
Pain 0.039 0.115 0.07  0.048 0     0.461
Muscles 0.107 0.044 0.04  0.07  0.107 0.393
Blood_Sugar 0.063 0.046 0     0.019 0.014 0.256
Code
plot(rez_efa$n_fac)

Code
rez_cfa9$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 611904.5 0.137 22,155.189 0.754 Inf
fit_ega 720615.2 0.041 2,068.422 0.980 NA
fit_egag 724475.5 0.063 6,068.587 0.941 0
fit_hclust 730031.8 0.089 11,594.850 0.886 0
fit_hclustg 731202.9 0.090 12,815.915 0.874 0
fit_efag 731887.4 0.091 13,520.449 0.867 0
fit_efa 731887.4 0.091 13,520.449 0.867 0

Initial

Code
rez_hclust <- make_hclust(data10)
rez_ega <- make_ega(data10, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Bruise"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data10, n = 4, cols = c("F1" = colors[["Hungry"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Itch"]], "F4" = colors[["Breathing"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5
Breathing 0.525 0.096 0.014 0.167 0.095
Heart 0.459 0.199 0.01  0.102 0.06 
Muscles 0.357 0.185 0.043 0.108 0.088
Hungry 0.148 0.49  0.018 0.02  0.025
Sex_arousal 0.031 0.418 0.141 0     0.022
Pain 0.085 0.351 0.144 0     0    
Thirsty 0.084 0.349 0.18  0.02  0.03 
Temp 0.113 0.329 0.047 0.16  0.105
Affective_touch 0.099 0.234 0.014 0.031 0.154
Defecate 0.032 0.085 0.457 0.356 0.046
Urinate 0.006 0.258 0.42  0.093 0.032
Vomit 0.024 0.084 0.37  0.251 0.125
Sneeze 0.136 0.039 0.209 0.462 0.039
Burp 0.076 0.074 0.222 0.462 0.119
Itch 0.091 0.058 0.052 0.087 0.474
Bruise 0.039 0.055 0.072 0.012 0.363
Wind 0.052 0.103 0.043 0.111 0.23 
Code
plot(rez_efa$n_fac)

Final

Note: No CFA (missing items)

Code
rez_hclust <- make_hclust(data10f)
rez_ega <- make_ega(data10f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Wind"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data10f, n = 3, cols = c("F1" = colors[["Breathing"]], "F2" = colors[["Defecate"]], 
                                            "F3" = colors[["Hungry"]]))
# rez_cfa10  <- make_cfa(data10f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4
Breathing 0.524 0.081 0.108 0.243
Heart 0.47  0.172 0.083 0.147
Muscles 0.372 0.175 0.108 0.096
Hungry 0.168 0.616 0.056 0.019
Thirsty 0.105 0.377 0.182 0.155
Pain 0.124 0.297 0.138 0.073
Defecate 0.049 0.089 0.528 0.099
Sneeze 0.149 0     0.391 0.153
Burp 0.103 0.06  0.364 0.184
Urinate 0.015 0.278 0.337 0.118
Wind 0.079 0.042 0.123 0.378
Bruise 0.079 0.045 0.049 0.378
Code
plot(rez_efa$n_fac)

Code
# rez_cfa10$test

Initial

Code
rez_hclust <- make_hclust(data11)
rez_ega <- make_ega(data11, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Bruise"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data11, n = 4, cols = c("F1" = colors[["Hungry"]], "F2" = colors[["Wind"]], 
                                            "F3" = colors[["Bruise"]], "F4" = colors[["Breathing"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.491 0.061 0.08  0.104 0.107 0.117
Breathing 0.491 0.101 0.021 0.178 0.067 0.194
Hungry 0     0.492 0.122 0.141 0.012 0.21
Thirsty 0.162 0.492 0.14  0.094 0.071 0.198
Defecate 0.036 0     0.477 0.044 0.113 0.066
Wind 0     0.012 0.404 0.095 0.257 0.072
Urinate 0.009 0.212 0.285 0.075 0.033 0.098
Sex_arousal 0.046 0.171 0.274 0.098 0.02  0.197
Vomit 0.02  0.037 0.274 0.093 0.08  0.084
Burp 0.055 0     0.265 0.188 0.248 0.017
Bruise 0.074 0.052 0.101 0.389 0.052 6.507 × 10−4
Itch 0.034 0.044 0.101 0.308 0.226 0.037
Blood_Sugar 0.104 0.072 0.003 0.26  0.016 0.019
Affective_touch 0.069 0.156 0.166 0.252 0.026 0.079
Taste 0.107 0     0.122 0.223 0.186 0.191
Sneeze 0.081 0.062 0.198 0.138 0.441 0.006
Cough 0.051 0     0.145 0.14  0.441 0.172
Pain 0.07  0.218 0.099 0.098 0.051 0.404
Muscles 0.163 0.02  0.143 0.114 0.207 0.355
Temp 0.136 0.245 0.141 0.068 0.022 0.346
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data11f)
rez_ega <- make_ega(data11f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Bruise"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Pain"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data11f, n = 3, cols = c("F1" = colors[["Wind"]], "F2" = colors[["Hungry"]], 
                                            "F3" = colors[["Breathing"]]))
rez_cfa11  <- make_cfa(data11f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4
Breathing 0.453 0.147 0.07  0.169
Heart 0.44  0.077 0.119 0.082
Muscles 0.257 0.173 0.251 0.047
Hungry 0.007 0.494 0.049 0.14 
Thirsty 0.221 0.485 0.058 0.081
Pain 0.257 0.318 0.078 0.02 
Urinate 0.032 0.312 0.211 0.015
Wind 0.084 0.087 0.476 0.086
Burp 0.109 0.021 0.386 0.106
Sneeze 0.099 0.063 0.373 0.044
Cough 0.187 0.065 0.32  0.068
Defecate 0.114 0.174 0.299 0    
Bruise 0.077 0.077 0.14  0.463
Blood_Sugar 0.127 0.057 0.016 0.463
Code
plot(rez_efa$n_fac)

Code
rez_cfa11$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 24319.43 0.088 401.264 0.843 Inf
fit_hclustg 29002.67 0.057 268.750 0.919 3.113615e+04
fit_egag 29008.58 0.057 261.203 0.922 1.621090e+03
fit_hclust 29010.84 0.056 243.273 0.928 5.246595e+02
fit_ega 29023.36 0.052 181.784 0.948 NA
fit_efa 29055.99 0.065 335.520 0.893 8.250859e-08
fit_efag 29055.99 0.065 335.520 0.893 8.250859e-08

Initial

Code
rez_hclust <- make_hclust(data12)
rez_ega <- make_ega(data12, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Bruise"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data12, n = 2, cols = c("F1" = colors[["Urinate"]], "F2" = colors[["Bruise"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)

Code
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing 0.416 0.161 0.053 0.152 0.003 0.122
Muscles 0.362 0.05  0.155 0     0.088 0.115
Heart 0.264 0.016 0.089 0.02  0.046 0.052
Affective_touch 0.148 0.091 0.092 0.057 0.043 0    
Hungry 0.025 0.367 0.177 0.015 0.014 0.197
Burp 0.09  0.347 0.066 0     0.144 0    
Temp 0.205 0.314 0.156 0.021 0.108 0.119
Blood_Sugar 0.01  0.212 0     0.029 0.038 0    
Urinate 0.118 0.161 0.552 0.01  0.188 0.015
Wind 0.126 0.089 0.262 0     0.117 0    
Defecate 0.018 0.065 0.246 0.093 0.137 0.154
Thirsty 0.159 0.1   0.236 0.077 0.027 0.099
Itch 0.066 0     0.01  0.535 0.1   0.035
Taste 0.142 0.021 0.117 0.352 0.025 0.15 
Bruise 0.03  0.044 0.047 0.286 0.044 0.125
Sneeze 0.09  0.138 0.295 0     0.554 0.019
Cough 0.083 0.095 0.016 0.134 0.348 0.03 
Vomit 0     0.049 0.109 0.023 0.206 0.082
Sex_arousal 0.217 0     0.17  0.116 0.09  0.516
Pain 0.051 0.281 0.06  0.16  0.035 0.516
Code
plot(rez_efa$n_fac)

Final

Note: CFA errors.

Code
rez_hclust <- make_hclust(data12f)
rez_ega <- make_ega(data8af, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Bruise"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data12f, n = 5, cols = c("F1" = colors[["Urinate"]], "F2" = colors[["Burp"]], 
                                             "F3" = colors[["Heart"]], "F4" = colors[["Pain"]], "F5" = colors[["Sneeze"]]))
# rez_cfa12  <- make_cfa(data12f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)

Code
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.531 0.072 0     0.045 0     0.128
Breathing 0.531 0.086 0.041 0.095 0.06  0.051
Hungry 0.146 0.533 0.107 0.052 0.059 0.044
Thirsty 0.014 0.533 0.122 0.063 0     0.138
Urinate 0.031 0.09  0.546 0.065 0.003 0.152
Defecate 0.013 0.154 0.546 0.102 0.089 0.09 
Sneeze 0.063 0.032 0.096 0.533 0.238 0.099
Cough 0.078 0.082 0.06  0.533 0.301 0.02 
Wind 0.054 0.053 0.018 0.199 0.509 0.095
Burp 0     0     0.059 0.28  0.509 0.162
Bruise 0.031 0     0.036 0     0.094 0.464
Muscles 0.043 0.013 0.015 0.026 0.121 0.401
Blood_Sugar 0.062 0     0.092 0.087 0.066 0.316
Pain 0.082 0.207 0.129 0.031 0.067 0.28 
Code
plot(rez_efa$n_fac)

Code
# rez_cfa12$test

Initial

Code
rez_hclust <- make_hclust(data13)
rez_ega <- make_ega(data13, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Sex arousal"]], "C6" = colors[["Sneeze"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data13, n = 4, cols = c("F1" = colors[["Thirsty"]], "F2" = colors[["Heart"]], 
                                            "F3" = colors[["Burp"]], "F4" = colors[["Bruise"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.621 0.033 0     0.154 0.009 0.029 0.084
Breathing 0.621 0     0.038 0.177 0.015 0.039 0.044
Hungry 0     0.514 0.151 0.089 0.123 0.072 0.046
Thirsty 0.02  0.514 0.145 0.228 0.083 0.015 0.023
Urinate 0     0.342 0.544 0.108 0.053 0.069 0    
Defecate 0.027 0     0.544 0.1   0.285 0.129 0    
Taste 0.061 0.009 0.07  0.391 0.11  0.084 0.092
Pain 0.05  0.15  0.003 0.358 0.209 0.029 0    
Temp 0.07  0.147 0.136 0.335 0.145 0.033 0.001
Muscles 0.069 0.085 0.014 0.274 0.146 0.064 0.05 
Affective_touch 0.014 0.085 0.027 0.186 0.41  0.001 0.042
Sex_arousal 0     0.118 0.164 0.219 0.397 0.035 0    
Vomit 0     0     0.097 0.083 0.24  0.14  0.038
Burp 0.016 0.069 0.017 0.033 0.102 0.608 0.015
Sneeze 0.031 0     0.101 0.115 0.082 0.461 0.011
Wind 0     0.05  0.143 0.037 0.036 0.363 0.012
Cough 0.019 0.018 0.007 0.081 0.059 0.321 0.119
Bruise 0     0.055 0     0.002 0.06  0.059 0.581
Itch 0.006 0     0     0.13  0     0.058 0.356
Blood_Sugar 0.085 0.024 0     0.002 0.034 0     0.253
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data13f)
rez_ega <- make_ega(data13f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Bruise"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data13f, n = 4, cols = c("F1" = colors[["Thirsty"]], "F2" = colors[["Heart"]], 
                                            "F3" = colors[["Burp"]], "F4" = colors[["Bruise"]]))
rez_cfa13  <- make_cfa(data13f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5
Heart 0.625 0.162 0.005 0.031 0.119
Breathing 0.625 0.089 0.057 0.06  0.068
Thirsty 0.034 0.536 0.154 0.024 0.051
Hungry 0     0.367 0.176 0.085 0.088
Pain 0.079 0.359 0.074 0.055 0    
Muscles 0.091 0.235 0.068 0.106 0.017
Urinate 0.004 0.31  0.562 0.088 0    
Defecate 0.043 0.13  0.562 0.166 0    
Burp 0.025 0.064 0.043 0.622 0    
Sneeze 0.036 0.09  0.115 0.472 0.042
Wind 0     0.096 0.138 0.361 0    
Cough 0.026 0.068 0.025 0.336 0.138
Bruise 0.008 0.061 0     0.094 0.479
Blood_Sugar 0.086 0.035 0     0     0.479
Code
plot(rez_efa$n_fac)

Code
rez_cfa13$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 21720.12 0.120 623.029 0.770 Inf
fit_hclust 25894.83 0.049 184.972 0.955 8.984228e+00
fit_ega 25899.23 0.038 116.734 0.977 NA
fit_hclustg 25903.37 0.054 226.518 0.942 1.260075e-01
fit_egag 25965.74 0.063 275.684 0.922 3.602767e-15
fit_efag 26137.94 0.086 474.294 0.849 1.458637e-52
fit_efa 26137.94 0.086 474.294 0.849 1.458637e-52

Initial

Code
rez_hclust <- make_hclust(data14)
rez_ega <- make_ega(data14, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data14, n = 7, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Thirsty"]], 
                                            "F3" = colors[["Bruise"]], "F4" = colors[["Temp"]], 
                                            "F5" = colors[["Pain"]], "F6" = colors[["Sneeze"]],
                                            "F7" = colors[["Heart"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7 8 9
Heart 0.557 0     0 0     0.142 0     0      0     0    
Breathing 0.557 0.052 0.026 0.093 0     0.062 0.12   0     0.138
Thirsty 0     0.632 0.114 0.106 0     0.027 0.026  0.044 0.113
Hungry 0.034 0.351 0.011 0.031 0.08  0.016 0.009  0     0    
Taste 0.019 0.28  0.043 0     0.208 0.041 0.089  0     0.172
Urinate 0     0.149 0.6 0.167 0     0.028 0      0     0    
Defecate 0.032 0.052 0.6 0.001 0     0.072 0.152 −0.043 0    
Sneeze 0     0     0 0.565 0.298 0     0.046  0     0.215
Vomit 0.022 0.098 9.556 × 10−4 0.322 0.048 0     0      0     0.038
Sex_arousal 0.066 0.029 0.13 0.302 0.049 0.14  0      0     0    
Cough 0.103 0     0 0.303 0.491 0.158 0.113  0.102 0.082
Temp 0     0.205 0 0     0.491 0     0      0     0.11 
Wind 0     0.081 0.026 0.161 0.137 0.574 0.065  0     0.089
Burp 0.067 0.008 0.063 0     0.1   0.574 0      0.03  0    
Muscles 0.123 0     0.037 0.053 0.168 0.03  0.573  0     0    
Pain 0.006 0.131 0.097 0     0     0.034 0.573  0.054 0.069
Bruise 0     0.027 −0.022 0     0.052 0     0.031  0.463 0.092
Itch 0     0     0 0     0.035 0.017 0      0.463 0    
Blood_Sugar 0     0.15  0 0.025 0.078 0.052 0      0.094 0.465
Affective_touch 0.087 0.026 0 0.143 0.089 0     0.041  0     0.465
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data14f)
rez_ega <- make_ega(data14f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Wind"]], "C6" = colors[["Pain"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data14f, n = 5, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Defecate"]], 
                                             "F3" = colors[["Thirsty"]], "F4" = colors[["Heart"]], 
                                             "F5" = colors[["Muscles"]]))
rez_cfa14  <- make_cfa(data14f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Heart 0.568  0      0     0.11 0     0    
Breathing 0.568  0.055  0.045 0.017 0.078 0.132
Thirsty 0      0.57   0.132 0 0.05  0.063
Hungry 0.052  0.364  0.018 0.024 0.025 0.011
Blood_Sugar 0      0.257  0     0.094 0.07  0    
Bruise 0      0.164 −0.051 0.054 0     0.054
Urinate 0.015  0.192  0.613 0.013 0.052 0    
Defecate 0.04  −0.065  0.613 −5.199 × 10−4 0.086 0.165
Sneeze 0      0.025  0.031 0.568 0.011 0.072
Cough 0.126  0.155 −0.021 0.568 0.188 0.116
Wind 0.014  0.163  0.051 0.137 0.583 0.066
Burp 0.069  0      0.071 0.077 0.583 0    
Muscles 0.124  0      0.039 0.192 0.025 0.58 
Pain 0.016  0.143  0.104 0.007 0.04  0.58 
Code
plot(rez_efa$n_fac)

Code
rez_cfa14$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 3379.581 0.152 187.096 0.637 2.344492e+129
fit_egag 3941.802 0.061 97.490 0.929 1.928811e+07
fit_hclustg 3972.347 0.092 137.380 0.831 4.494018e+00
fit_ega 3975.352 0.040 65.620 0.975 NA
fit_hclust 3983.183 0.090 124.852 0.851 1.993355e-02
fit_efa 3985.451 0.104 159.830 0.779 6.413983e-03
fit_efag 3985.451 0.104 159.830 0.779 6.413983e-03

Initial

Code
rez_hclust <- make_hclust(data15)
rez_ega <- make_ega(data15, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Sex arousal"]], "C4" = colors[["Sneeze"]], 
  "C5" = colors[["Muscles"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data15, n = 3, cols = c("F1" = colors[["Sneeze"]], "F2" = colors[["Hungry"]], 
                                            "F3" = colors[["Bruise"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)

Code
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Breathing 0.604  0.184 0      0     0.029 0.037  0.106
Hungry 0.425  0.236 0.028  0     0.004 0      0.064
Heart 0.352  0.024 0.219  0.055 0     0      0.081
Pain 0.066  0.491 0.389  0     0.186 0      0.177
Urinate 0.105  0.425 0.031  0     0.156 0     −0.06 
Defecate 0.069  0.328 0.142  0.145 0     0     −0.049
Thirsty 0.207  0.314 0     −0.067 0.056 0.06   0.042
Taste 0.028  0.122 0.438  0.174 0     0      0.063
Sex_arousal 0.089  0.142 0.438  0.057 0.055 0.176  0.047
Sneeze 0      0.15  0.124  0.716 0.076 0      0    
Vomit 0.019  0     0.196  0.313 0.121 0.126  0.055
Cough 0.038  0     0.13   0.3   0.09  0.211  0    
Itch 0     −0.069 0.059  0.288 0.014 0.21   0.102
Temp 0      0.196 0.096  0.272 0.543 0.052  0    
Muscles 0.027  0.13  0     −0.033 0.543 0.116  0.291
Wind 0      0.02  0.293  0.198 0.02  0.533  0    
Burp 0.03   0.026 0      0.216 0.14  0.533  0.176
Blood_Sugar 0     −0.081 0      0     0.139 0      0.479
Bruise 0.126  0.04  0      0.019 0.115 0.134  0.369
Affective_touch 0.06   0.122 0.172  0.094 0.008 0.032  0.242
Code
plot(rez_efa$n_fac)

Final

Note: CFA errors.

Code
rez_hclust <- make_hclust(data15f)
rez_ega <- make_ega(data15f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Sneeze"]], "C4" = colors[["Wind"]], 
  "C5" = colors[["Bruise"]], "C6" = colors[["Wind"]], "C7" = colors[["Pain"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data15f, n = 3, cols = c("F1" = colors[["Hungry"]], "F2" = colors[["Burp"]], 
                                            "F3" = colors[["Blood Sugar"]]))
# rez_cfa15  <- make_cfa(data15f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)

Code
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5
Breathing 0.606  0.181 0     0.036  0.119
Hungry 0.428  0.231 0     0      0.058
Heart 0.368  0.067 0.075 0.042  0    
Pain 0.104  0.579 0     0.079  0.173
Urinate 0.119  0.453 0.027 0      0.032
Defecate 0.078  0.345 0.223 0.009 −0.043
Thirsty 0.203  0.292 0     0.035  0.074
Sneeze 0      0.22  0.571 0.127  0    
Cough 0.069  0     0.571 0.234  0.061
Wind 0.037  0.072 0.104 0.557  0.035
Burp 0.031  0.03  0.236 0.557  0.233
Blood_Sugar 0.013 −0.077 0     0      0.492
Bruise 0.125  0.055 0     0.135  0.423
Muscles 0.023  0.226 0.06  0.146  0.382
Code
plot(rez_efa$n_fac)

Code
# rez_cfa15$test

Initial

Code
rez_hclust <- make_hclust(data16)
rez_ega <- make_ega(data16, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Bruise"]], "C7" = colors[["Wind"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data16, n = 4, cols = c("F1" = colors[["Thirsty"]], "F2" = colors[["Sneeze"]], 
                                            "F3" = colors[["Blood Sugar"]], "F4" = colors[["Burp"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Heart 0.478 0.037  0.017 0.094  0.044  0.085 0.035
Breathing 0.478 0.189  0.068 0.136  0.021  0.059 0
Thirsty 0.088 0.571  0.142 0      0      0.09  0.044
Hungry 0.172 0.466  0.125 0.013  0      0.125 7.113 × 10−4
Sex_arousal 0.081 0.226  0.091 0      0.02   0.185 0.036
Urinate 0.038 0.28  0.555 0      0.02   0.066 0.04
Defecate 0.087 0.069  0.555 0.091  0.134  0.045 0.018
Taste 0.041 0  0.073 0.512  0.089  0.021 0.085
Muscles 0.173 0.012  0.028 0.395  0.077  0.212 0
Itch 0.105 0 −0.015 0.277  0.074  0.164 0.011
Sneeze 0.003 0  0.02  0.133  0.585  0.13  0
Cough 0.099 0  0.039 0.113  0.47  −0.011 0.145
Vomit 0.004 0.022  0.112 0.036  0.292  0.042 0.053
Pain 0.236 0.086  0     0.146  0.063  0.462 0
Bruise 0     0  0     0.097  0.103  0.421 0.069
Blood_Sugar 0     0.04  0     0.029 −0.029  0.349 0
Affective_touch 0     0.223  0     0.079  0      0.311 0
Temp 0.015 0.113  0.131 0.147  0.036  0.232 0.086
Wind 0.03  0.098  0.046 0.102  0.053  0.106 0.604
Burp 0.034 8.632 × 10−4  0.026 0.024  0.168  0.056 0.604
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data16f)
rez_ega <- make_ega(data16f, cols = c(
  "C1" = colors[["Heart"]], "C2" = colors[["Hungry"]], "C3" = colors[["Defecate"]], "C4" = colors[["Muscles"]], 
  "C5" = colors[["Cough"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data16f, n = 4, cols = c("F1" = colors[["Thirsty"]], "F2" = colors[["Cough"]], 
                                            "F3" = colors[["Bruise"]], "F4" = colors[["Wind"]]))
# rez_cfa16  <- make_cfa(data16f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6
Breathing 0.429 0.145 0.078  0.03  0      0    
Pain 0.416 0.093 0.035  0.03  0.031  0.242
Muscles 0.378 0.037 0.075  0.134 0.013  0.011
Heart 0.342 0.038 0.021  0.048 0.046  0    
Hungry 0.18  0.567 0.133  0     0.008  0.081
Thirsty 0.108 0.567 0.162  0     0.06   0    
Urinate 0.04  0.265 0.568  0.025 0.057  0    
Defecate 0.151 0.03  0.568  0.08  0.042  0    
Sneeze 0.05  0     0.062  0.574 0.027  0.156
Cough 0.179 0     0.045  0.574 0.152 −0.026
Wind 0.076 0.073 0.088  0.075 0.612  0.044
Burp 0.025 0.009 0.033  0.135 0.612  0.068
Bruise 0.193 0     0      0.126 0.077  0.529
Blood_Sugar 0     0.067 0     −0.021 0      0.529
Code
plot(rez_efa$n_fac)

Code
# rez_cfa16$test

Initial

  • the first optimal number of factors was 11, which, although too large for practical higher-order structure description, hints at the scattered nature of the items
Code
rez_hclust <- make_hclust(data17)
rez_ega <- make_ega(data17, cols = c(
  "C1" = colors[["Hungry"]], "C2" = colors[["Defecate"]], "C3" = colors[["Muscles"]], "C4" = colors[["Cough"]], 
  "C5" = colors[["Itch"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data17, n = 5, cols = c("F1" = colors[["Sneeze"]], "F2" = colors[["Hungry"]], 
                                            "F3" = colors[["Defecate"]], "F4" = colors[["Burp"]],
                                            "F5" = colors[["Vomit"]]))

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5 6 7
Hungry 0.597 0.04  0.032 0     0     1.557 × 10−4 0    
Breathing 0.407 0     0.088 0.046 0     0.02 0.096
Thirsty 0.337 0.125 0.065 0     0.043 0.09 0    
Blood_Sugar 0.224 0     0.105 0.044 0     0.037 0    
Heart 0.211 0.057 0.092 0.039 0     0.05 0.12 
Urinate 0.14  0.557 0.144 0.026 0.058 0.012 0.003
Defecate 0.049 0.557 0.126 0.065 0.068 0.021 0    
Sex_arousal 0.034 0.18  0.528 0     0.003 0.109 0    
Affective_touch 0.025 0.043 0.306 0.111 0.067 0 0.114
Muscles 0.131 0     0.297 0.053 0.122 0.127 0.182
Taste 0.14  0.05  0.239 0.109 0.015 0.084 0    
Sneeze 0.047 0.079 0.006 0.699 0.05  0.03 0.056
Cough 0.041 0     0.114 0.504 0.174 0.068 0.029
Vomit 0.038 0.026 0.191 0.195 0.034 0.017 0.151
Temp 0.03  0.05  0.112 0.122 0.514 0.052 0.163
Itch 0     0.052 0.054 0.059 0.514 0.044 0.208
Wind 0.089 0.033 0.215 0.01  0.118 0.557 0    
Burp 0.08  0     0.101 0.09  0     0.557 0.282
Bruise 0.034 0     0     0.055 0.081 0.209 0.496
Pain 0.102 0.002 0.216 0.097 0.258 0 0.496
Code
plot(rez_efa$n_fac)

Final

Code
rez_hclust <- make_hclust(data17f)
rez_ega <- make_ega(data17f, cols = c(
  "C1" = colors[["Hungry"]], "C2" = colors[["Defecate"]], "C3" = colors[["Cough"]], "C4" = colors[["Burp"]], 
  "C5" = colors[["Muscles"]], "C6" = colors[["Wind"]], "C7" = colors[["Bruise"]], "C8" = colors[["Blood Sugar"]]))
rez_efa <- make_efa(data17f, n = 4, cols = c("F1" = colors[["Burp"]], "F2" = colors[["Hungry"]], 
                                            "F3" = colors[["Urinate"]], "F4" = colors[["Sneeze"]]))
rez_cfa17  <- make_cfa(data17f)

rez_hclust$p / rez_ega$p / rez_efa$p

Code
plot(rez_hclust$pvclust, hang = -1, cex = 0.5)
pvclust::pvrect(rez_hclust$pvclust, alpha=0.95, max.only=FALSE)

Code
stab <- EGAnet::itemStability(rez_ega$rez)

Code
rez_ega$table
EGA Loadings
Item 1 2 3 4 5
Hungry 0.599 0.045 0     0.002 0.028
Breathing 0.414 0     0.034 0.015 0.179
Thirsty 0.335 0.128 0.003 0.099 0.058
Blood_Sugar 0.242 0     0.031 0.055 0    
Heart 0.213 0.069 0.039 0.058 0.132
Urinate 0.158 0.575 0.033 0.038 0.051
Defecate 0.063 0.575 0.059 0.061 0.014
Sneeze 0.05  0.113 0.621 0.04  0.093
Cough 0.071 0     0.621 0.096 0.147
Wind 0.104 0.094 0.029 0.564 0.132
Burp 0.095 0     0.078 0.564 0.304
Pain 0.108 0.046 0.098 0.012 0.613
Bruise 0.048 0     0     0.219 0.343
Muscles 0.172 0.013 0.081 0.185 0.28 
Code
plot(rez_efa$n_fac)

Code
rez_cfa17$test
CFA Model Comparison
Name BIC RMSEA Chi2 CFI BF
fit_g 7047.020 0.108 214.312 0.784 Inf
fit_ega 8493.471 0.057 102.289 0.947 NA
fit_efag 8516.023 0.089 224.513 0.828 1.267374e-05
fit_efa 8516.023 0.089 224.513 0.828 1.267374e-05
fit_hclust 8557.833 0.097 227.561 0.817 1.056922e-14
fit_egag NA NA NA NA 3.324503e+10
fit_hclustg NA NA NA NA 3.130020e-10

Scores

Note the usage of descriptive factor names relating directly to the items to avoid abstraction.

Code
# Add empirical variables
add_empirical <- function(data, sample = "Sample 1a") {
  x <- data.frame(
    Original = rowMeans(data),
    HungryThirsty = (data$Hungry + data$Thirsty) / 2,
    MusclesPain = (data$Muscles + data$Pain) / 2,
    WindBurp = (data$Wind + data$Burp) / 2,
    UrinateDefecate = (data$Urinate + data$Defecate) / 2,
    BreathingHeart = (data$Heart + data$Breathing) / 2)
  
  if("Blood_Sugar" %in% names(data)) {
    x$BruiseBlood <- (data$Blood_Sugar + data$Bruise) / 2
  } else {
    x$BruiseBlood <- data$Bruise
  }
  
  if("Cough" %in% names(data)) {
    x$CoughSneeze <- (data$Cough + data$Sneeze) / 2
  } else {
    x$CoughSneeze <- data$Sneeze
  }
  x$Sample <- sample
  x
}

scores <- list(
  sample1a = add_empirical(data1a, "Sample 1a"),
  sample1b = add_empirical(data1b, "Sample 1b"),
  sample2 = add_empirical(data2, "Sample 2"),
  sample3 = add_empirical(data3, "Sample 3"),
  sample4 = add_empirical(data4, "Sample 4"),
  sample5 = add_empirical(data5, "Sample 5"),
  sample6 = add_empirical(data6, "Sample 6"),
  sample7a = add_empirical(data7a, "Sample 7a"),
  sample7b = add_empirical(data7b, "Sample 7b"),
  sample7c = add_empirical(data7c, "Sample 7c"),
  sample8a = add_empirical(data8a, "Sample 8a"),
  sample8b = add_empirical(data8b, "Sample 8b"),
  sample9 = add_empirical(data9, "Sample 9"),
  sample10 = add_empirical(data10, "Sample 10"),
  sample11 = add_empirical(data11 , "Sample 11"),
  sample12 = add_empirical(data12 , "Sample 12"),
  sample13 = add_empirical(data13, "Sample 13"),
  sample14 = add_empirical(data14 , "Sample 14"),
  sample15 = add_empirical(data15 , "Sample 15"),
  sample16 = add_empirical(data16, "Sample 16"),
  sample17 = add_empirical(data17, "Sample 17")
)

save(scores, file = "../data/scores.RData")