#                                                                              #         
#              The Importance of Using an Optimal Cutoff Value for             # 
#                  the 10-Item Autism-Spectrum Quotient (AQ10)                 #           
#                                                                              #
#                        Psychological Assessment, 2024                        #
#                                                                              #
################################################################################        


###############################**Load Packages**################################
packages <- c("boot", "caret", "DescTools", "dplyr", "ggplot2", "lme4", "psych", "report", "Rmisc", "rstatix", "sjPlot")
for (package in packages){
  if(!is.element(package, .packages(all.available = TRUE))){
    install.packages(package)}
  library(package, character.only = TRUE)}

data <- read.csv(Sys.glob("C:\\Users\\*\\Downloads\\Analysis_Data.csv"))        #Data is available in open repositories or from the corresponding authors upon request

str(data)                                                                       #See READ_ME.txt for data structure summary
###############################**Pre-Processing**###############################                   
##**Classify By Cut-Off**##
cut_offs <- c("Cut6", "Cut7")                                                   #Classify whole data set above/below the 6+ and 7+ cut-off

for(x in 1:2){
  value = 4+x                                                                   
  data[[cut_offs[x]]] <- ifelse(data$AQ10 > value, 1, 0)}                       #1 = above cut-off, 0 = below cut-off          

##**Remove Missing Data**##
dat <- filter(data, Study != "Thabtah_2017")                                    #No DV
dat <- filter(dat, Study != "Rahman_2021")                                      #No ppts with an AQ10 score 7+

##**Inspect Data**##                                                            
describe(data)                                                                  #See READ_ME.txt for data summary statistics

length(unique(data$ppt))                                                        #13692 ppts across 25 studies
length(unique(dat$ppt))                                                         #13070 ppts included in Effect Size Analyses

############################**AQ10 Group Analyses**############################# 
##**Descriptives, Sensitivity, and Specificity**##                  
Classifications <- data.frame(Study = NA,  L6_N = NA, L6_M = NA, L6_SD = NA,    #Create empty data frame to store AQ10 group descriptives 
                                           H6_N = NA, H6_M = NA, H6_SD = NA,    #N = no. ppts, M = Mean AQ10, SD = SD of Mean
                                           L7_N = NA, L7_M = NA, L7_SD = NA,    #L6/H6 = low traits/high traits classified using the 6+ cut-off    
                                           H7_N = NA, H7_M = NA, H7_SD = NA,    #L7/H7 = low traits/high traits classified using the 7+ cut-off   
                                           ND_N = NA, ND_M = NA, ND_SD = NA,    #ND = No Diagnosis
                                           D_N  = NA, D_M  = NA, D_SD  = NA)    #D = Diagnosis

SenSpec <- data.frame(Study = NA, Cut = NA, Sensitivity = NA, Specificty = NA)  #Create empty data frame to store sensitivity and specificity information

class_data <- data                                                              #Adjust labels to represent separate data sets within 1 study
class_data$Study[class_data$Variable_Name == "Environmental_Behaviour"] <- "Taylor_2021_2"     
class_data$Study[class_data$Variable_Name == "Climate Belief"] <- "Taylor_2021_3"

cols <- list(L6 = c(2:4), H6 = c(5:7), L7 = c(8:10), H7 = c(11:13))             #List of column indices to store low/ligh 6+/7+ data

for(study in unique(class_data$Study)){                                         
  s <- filter(class_data, Study == study)                                       #Filter to the data unique to each study data set
  s <- s[!(duplicated(s$ppt)),]                                                 #Remove duplicated information for studies with multiple variables
  #Sensitivity & Specificity 
  if(is.na(s$Diagnosis[1]) == F){                                               #Can only calculate sensitivity and specificity on data sets with Diagnosis information
    for(cut in cut_offs){                                                         #For each cut-off value:
      cut_grps <- table(s[[cut]], s$Diagnosis)                                      #Create a cross-tabs table original cut-off based groups vs Diagnosis
      SenSpec[nrow(SenSpec) +1, 1] <- study                                         #Import Study name
      SenSpec[nrow(SenSpec), 2] <- cut                                              #Specify which cut-off value was used to define groups
      SenSpec[nrow(SenSpec), 3] <- sensitivity(cut_grps)                            #Sensitivity of cut-off value
      SenSpec[nrow(SenSpec), 4] <- specificity(cut_grps)}}                          #Specificity of cut-off value
  #Group Descriptives
  Classifications[nrow(Classifications) +1, 1] <- s$Study[1]                    #Import Study name
  NoDiag <- filter(s, Diagnosis == 0)                                           #Import descriptives (N, Mean AQ10, SD AQ10) for Diagnosis and No Diagnosis groups
  Diag <- filter(s, Diagnosis == 1)
  Classifications[nrow(Classifications), 14:16] <- c(nrow(NoDiag), mean(NoDiag$AQ10), sd(NoDiag$AQ10)) 
  Classifications[nrow(Classifications), 17:19] <- c(nrow(Diag), mean(Diag$AQ10), sd(Diag$AQ10))
    for(cut in cut_offs){                                                       #Repeat for high/low groups based on each cut-off
      y <- ifelse(cut == "Cut6", list(cols[1:2]), list(cols[3:4]))              #Use 'cols' list to import into correct columns
      High <- filter(s, s[[cut]] == 1)
      Low  <- filter(s, s[[cut]] == 0)
      Classifications[nrow(Classifications), unlist(y[[1]][1])] <- c(nrow(Low), mean(Low$AQ10), sd(Low$AQ10))
      Classifications[nrow(Classifications), unlist(y[[1]][2])] <- c(nrow(High), mean(High$AQ10), sd(High$AQ10))}}

Classifications <- Classifications[-c(1), ]                                     #Remove empty NA row

##**Percentage Ratios**##
Classifications$TotalN <- Classifications$L6_N + Classifications$H6_N           #Calculate Total N for each data set 

for(col in c("L6", "H6", "L7", "H7")){                                          #Calculate percentage of total sample for each high/low group
  Classifications[[paste0(col, "_P")]] <- Classifications[[paste0(col, "_N")]]/ Classifications$TotalN}

Classifications$Reclassify <-  (Classifications$H7_N -  Classifications$H6_N) / Classifications$H6_N    #Percentage difference between the 6+ and 7+ high group N

describe(Classifications)
SenSpec %>% filter(Cut == "Cut6") %>% describe()         #Due to format of data, sensitivity and specificity are reversed 
SenSpec %>% filter(Cut == "Cut7") %>% describe()         #View more decimal places

##**T-Tests**##
tdata <- data.frame(Low = c(Classifications$L6_M, Classifications$L7_M),        #Re-structure data for t-test analyses
                    High = c(Classifications$H6_M, Classifications$H7_M),       
                    Cut = c(rep("Cut6", 25), rep("Cut7", 25)), 
                    Study = rep(Classifications$Study, 2))

tdata <- tdata %>% filter(Study != "Rahman_2021")                               #Remove data set with no. ppts in the 7+ high trait group

set.seed(1234)  
t.test(tdata$Low[tdata$Cut == "Cut6"], tdata$Low[tdata$Cut == "Cut7"], paired = T)                                            #10,000 bootstrapped 95% CI of Cohen's d                                                
tdata %>% cohens_d(Low ~ Cut, paired = T, ci = TRUE, conf.level = 0.95, nboot = 10000) 

set.seed(1234)                                                                  #Repeat for high AQ10 Mean
t.test(tdata$High[tdata$Cut == "Cut6"], tdata$High[tdata$Cut == "Cut7"], paired = T)                                           
tdata %>% cohens_d(High ~ Cut, paired = T, ci = TRUE, conf.level = 0.95, nboot = 10000)

##########################**Sex Comparison Analyses**###########################
#Create data frame to store total n by sex, no. participants who meet 6+ and 7+ cut-off by sex, AQ10 mean by sex, and AQ10 mean comparison test results
sex_stats <- data.frame(Data = NA, Female = NA, Male= NA, Female_6= NA, Male_6= NA, Female_7= NA, Male_7= NA, Female_Mean = NA, Male_Mean = NA, t= NA, df= NA, p = NA)

#Function to extract above information for each data set 
get_stats <- function(data, study, var){
  if(study == "Thabtah_2017"){d <- data %>% filter(Study == study)} else{d <- data %>% filter(Study == study) %>% filter(Variable_Name == var)}
  if(study == "Talipski_2020"){ ttest <- list(statistic = NA, parameter = NA, p.value = NA)} else{ttest <- t.test(d$AQ10~ d$Sex)}
  sex_stats <- data.frame(Data = study,
                          Female = table(d$Sex)[[1]],
                          Male = ifelse(length(table(d$Sex)) < 2 , 0, table(d$Sex)[[2]]),
                          Female_6 = table(d$Sex[d$Cut6 == 1])[[1]],
                          Male_6 = ifelse(length(table(d$Sex[d$Cut6 == 1])) < 2 , 0, table(d$Sex[d$Cut6 == 1])[[2]]),
                          Female_7 = ifelse(length(table(d$Sex[d$Cut7 ==1])) < 1, 0, table(d$Sex[d$Cut7 == 1])[[1]]),
                          Male_7 = ifelse(length(table(d$Sex[d$Cut7 == 1])) < 2 , 0, table(d$Sex[d$Cut7 == 1])[[2]]),
                          Female_Mean = d %>% filter(Sex == 0) %>% summarise(Female_Mean = mean(AQ10)),
                          Female_Mean = d %>% filter(Sex == 1) %>% summarise(Male_Mean = mean(AQ10)),
                          t = ttest$statistic,
                          df = ttest$parameter,
                          p = ttest$p.value)
  return(sex_stats)}

#Run function for each data set with sex information (isolating 1 variable to avoid duplication) and append to sex_stats
studies <- c("Brewer_2017", "Brewer_2022", "Clutterbuck_2021", "Clutterbuck_2022", "Hochard_2020", "Huang_2021", "Kember_2021", "Livingston_2020", "Lorenz_2021", "Nijhof_2017", 
             "Rahman_2021", "Rogge_2021", "Rogge_2022", "Ryan_2021", "Shand_2020", "Talipski_2020", "Taylor_2021", "Taylor_2021", "Taylor_2021", "Thabtah_2017")
variables <- c("AToM_Physical", "Physical_AToM", "FIMI", "Anthropomorphism", "Sleep", "Anxiety", "RAADS", "Compensation", "BAPQ", "PP", "KOR_Known", "SunkBias", "Fear_Max", "Alexithymia", 
               "Attitudes", "Angry_GCE", "Environmental_Behaviour", "Environmental_Attitudes", "Climate Belief", "" )

for(df in 1:length(studies)){
  stats <- get_stats(data, studies[df], variables[df])
  sex_stats <- rbind(sex_stats, stats)}

#Append information combined data set results
all <- data[!(duplicated(data$ppt)),]
ttest <- t.test(all$AQ10 ~ all$Sex)

sex_stats[nrow(sex_stats) + 1, c(1:12)] <- c("All_Data", table(all$Sex)[[1]], table(all$Sex)[[2]], table(all$Sex[all$Cut6 == 1])[[1]], table(all$Sex[all$Cut6 == 1])[[2]], 
                                             table(all$Sex[all$Cut7 == 1])[[1]], table(all$Sex[all$Cut7 == 1])[[2]], all %>% filter(Sex == 0) %>% summarise(Female_Mean = mean(AQ10)), 
                                             all %>% filter(Sex == 1) %>% summarise(Male_Mean = mean(AQ10)), ttest$statistic, ttest$parameter, ttest$p.value)
#Remove NA row
sex_stats <- sex_stats[-c(1), ]

#Add % proportion information
sex_stats$F6 <- (sex_stats$Female_6 / sex_stats$Female) * 100; sex_stats$M6 <- (sex_stats$Male_6 / sex_stats$Male) * 100

sex_stats$F7 <- (sex_stats$Female_7 / sex_stats$Female) * 100; sex_stats$M7 <- (sex_stats$Male_7 / sex_stats$Male) * 100

#Add N-1 Chi Square Statistics, to assess whether the proportion of participants excluded from the high 
#autistic traits group between the ≥6 and ≥7 cut-off was significantly different between males and females. 
sex_stats$X2 <- NA; sex_stats$df2 <- NA; sex_stats$X2_p <- NA

for(row in 1:nrow(sex_stats)){
  diff <- c((sex_stats$Female_6[row] - sex_stats$Female_7[row]), (sex_stats$Male_6[row] - sex_stats$Male_7[row]))
  dx2 <- data.frame(Female = c(diff[1], (sex_stats$Female[row] - diff[1])), Male = c(diff[2], (sex_stats$Male[row] - diff[2])))
  row.names(dx2) <- c("Excluded","Not_Excluded")
  chi <- chisq.test(dx2, correct = F) 
  sex_stats$X2[row] <- chi$statistic
  sex_stats$X2_p[row] <- chi$p.value
  sex_stats$df2[row] <- chi$parameter
}

print(sex_stats)
############################**Effect Size Analysis**############################
##**Main Analysis**##
for(variable in unique(dat$Variable_Name)){                                     #Filter to each variable (some Study data sets have more than one variable)
  var <- filter(data, Variable_Name == variable)
  ts <- list(t6 = t(unlist(t.test(DV_Z ~ Cut6, var, var.equal = F))),           #T-test of standardized DV between the high/low cut-off based groups
             t7 = t(unlist(t.test(DV_Z ~ Cut7, var, var.equal = F))))           #Compute separately for the 6+ and 7+ cut-offs
  if(exists("Comparisons") == FALSE){Comparisons <- as.data.frame(ts$t6)}          
  for (cut in 1:2){                                                             #Collect t-test outputs in new data frame
    Comparisons[nrow(Comparisons) +1, 1:12] <- ts[[cut]]
    Comparisons$CutOff[nrow(Comparisons)] <- cut_offs[cut]                      #Specify cut-off used to define groups
    Comparisons$Study[nrow(Comparisons)] <- var$Study[1]                        #Name of Study data set
    Comparisons$Variable_Name[nrow(Comparisons)] <- var$Variable_Name[1]}}      #And name of variable within Study
 
Comparisons <- Comparisons[-c(1), ]                                             #Remove duplicate row

TTest <- Comparisons %>% select(Study, Variable_Name, CutOff, t = statistic.t,  #Select & rename columns fo interest
                                df = parameter.df, p = p.value,
                                Low_Mean = `estimate.mean in group 0`, 
                                High_Mean = `estimate.mean in group 1`,SE = stderr)

TTest$p <- as.numeric(TTest$p)
TTest$Sig <- ifelse(TTest$p == 0.05 | TTest$p < 0.05, "Sig", "NonSig")          #Create columns to easily compare changes in significance

for (col in c('Low_Mean', 'High_Mean', 't', 'df', 'SE')){
  TTest[[col]] <- unlist(as.numeric(TTest[[col]]))}

TTest$Diff <- TTest$High_Mean - TTest$Low_Mean                                  #Calculate difference in High/Low group DV Mean for each cut-off 

TTest$DiffDiff <- NA
TTest$Direction <- NA
TTest$SameSig <- NA

row6 <- seq(1,(nrow(TTest) - 1), by = 2)                                        #Vectors to separately index the 6+ and 7+ rows
row7 <- seq(2, nrow(TTest), by = 2)

for (x in 1:length(row6)){
  TTest$DiffDiff[row6[x]] <- TTest$Diff[row6[x]] - TTest$Diff[row7[x]]          #Calculate Difference of Differences (6+ Mean difference - 7+ Mean Difference)
  TTest$Direction[row6[x]] <- sign(TTest$High_Mean[row6[x]]) == sign(TTest$High_Mean[row7[x]]) #Does the direction of effect change as a function of cut-off (FALSE = change in effect direction)
  TTest$SameSig[row6[x]] <- TTest$Sig[row6[x]] == TTest$Sig[row7[x]]}           #Does t-test significance change as a function of cut-off 

TTest$DiffDiff <- abs(TTest$DiffDiff)                                           #Calculate absolute difference of differences

describe(TTest)                                                                 
   
mainboot <- boot(na.omit(TTest$DiffDiff),function(u,i) mean(u[i]), R=10000)     #10,000 bootstrapped 95% CI of mean Difference of Differences
boot.ci(mainboot, type=c("norm","basic","perc"))

  ##########**Mixed Model and Intraclass Correlation Robustness Check**###########
model_data <- na.omit(TTest)                                                    #Only Interested in DiffDiff which is only stored in "Cut6" rows

model_data$Variable_Name[model_data$Variable_Name == "AToM_Physical"] <- "Physical_AToM"   #Keep labels consistent across variables where the same measure was used across studies
model_data$Variable_Name[model_data$Variable_Name == "AToM_Social"] <-  "Social_AToM"              

model_data$Study[model_data$Variable_Name == "Environmental_Behaviour"] <- "Taylor_2021_2" #Adjust labels to represent separate data sets within 1 study
model_data$Study[model_data$Variable_Name == "Climate Belief"] <- "Taylor_2021_3"               

model <- lmer(DiffDiff ~ 1 + (1|Study) + (1|Variable_Name), model_data, control=lmerControl(optimizer="bobyqa", optCtrl=list(maxfun=2e5)))
tab_model(model)       
summary(model)

ci <- confint.merMod(model, level = 0.95, method = "boot", nsim = 10000, boot.type = "perc") 
ci
####################################**END**#####################################
session <- sessionInfo()
report_system(session)
as.data.frame(report(session))