###################################################################################################################
####“MIT License Copyright (c) 2026 Silvia Volpe, Rita Sorrentino, Martina De Cola, Davide Mameli, ###############
########Teresa Nicolosi, Carla Figus, Annalisa Pietrobelli, Maria Giovanna Belcastro” ############################
###################################################################################################################
#################          Log file and full codes for analyses contained in the paper              ###############
#################          Volpe, Sorrentino, De Cola, Mameli, Nicolosi, Figus, Pietrobelli         ##############
#################                and Belcastro                                                      ###############
#################               'Sexual dimorphism and allometric variation in the human            ###############
################                     coxal bone from two coeval Italian skeletal samples'           ##############         
#################                                                                                   ###############
#################                              submitted to Humans                                  ###############
#################                                                                                   ###############
#################                                  October 2026                                     ###############
#################                                                                                   ###############
#################                                   authors:                                        ###############
#################             Silvia Volpe, Rita Sorrentino, Martina De Cola, Davide Mameli,        ###############
#################                    Teresa Nicolosi, Carla Figus, Annalisa Pietrobelli             ##############
################                               and Maria Giovanna Belcastro                         ###############
#################                                                                                   ###############
#################                              silvia.volpe5@unibo.it                               ###############
#################                            rita.sorrentino01@unipa.it                             ###############
#################                                decola14@libero.it                                 ###############
#################                             davide.mameli3@unibo.it                               ###############
#################                          teresa.nicolosi2@unibo.it                                ###############
#################                             carla.figus3@unibo.it                                 ###############
#################                        annalisa_pietrobelli@eva.mpg.de                            ###############
#################                       mariagiovanna.belcastro@unibo.it                            ###############
#################                                                                                   ###############
###################################################################################################################
###################################################################################################################

# using R version 4.5.0

## 1) ICC (Intraclass Correlation Coefficient)

install.packages("irr")
library(psych)
library("irr")

#https://www.r-bloggers.com/2021/06/intraclass-correlation-coefficient-in-r-quick-guide/

## 1.1) interobserver-ICC  
femur.sample.inter <- read.table(file.choose(),header=T,row.names = 1,sep="\t")

Femur_sx<-femur.sample.inter[,c(4,8)]
Femur_sx<-Femur_sx[complete.cases(Femur_sx),]
icc(Femur_sx, model= "twoway", type="agreement", unit="single") 

Femur_dx<-femur.sample.inter[,c(5,9)]
Femur_dx<-Femur_dx[complete.cases(Femur_dx),]
icc(Femur_dx, model= "twoway", type="agreement", unit="single")  

FemurHead_sx<-femur.sample.inter[,c(6,10)]
FemurHead_sx<-FemurHead_sx[complete.cases(FemurHead_sx),]
icc(FemurHead_sx, model= "twoway", type="agreement", unit="single") 

FemurHead_dx<-femur.sample.inter[,c(7,11)]
FemurHead_dx<-FemurHead_dx[complete.cases(FemurHead_dx),]
icc(FemurHead_dx, model= "twoway", type="agreement", unit="single")  

## 1.2) Intraobserver-ICC 
femur.sample.intra <- read.table(file.choose(),header=T,row.names = 1,sep="\t")

Femur_sx<-femur.sample.intra[,c(4,8)]
Femur_sx<-Femur_sx[complete.cases(Femur_sx),]
icc(Femur_sx, model= "twoway", type="agreement", unit="single")  

Femur_dx<-femur.sample.intra[,c(5,9)]
Femur_dx<-Femur_dx[complete.cases(Femur_dx),]
icc(Femur_dx, model= "twoway", type="agreement", unit="single")  

FemurHead_sx<-femur.sample.intra[,c(6,10)]
FemurHead_sx<-FemurHead_sx[complete.cases(FemurHead_sx),]
icc(FemurHead_sx, model= "twoway", type="agreement", unit="single")  

FemurHead_dx<-femur.sample.intra[,c(7,11)]
FemurHead_dx<-FemurHead_dx[complete.cases(FemurHead_dx),]
icc(FemurHead_dx, model= "twoway", type="agreement", unit="single") 

#######################################################################################################
## 2 ASYMMETRY test between the measurements performed on right and left femora for FEM and FHSI
femur.sample <- read.table(file.choose(),header=T,row.names = 1,sep="\t")

shapiro.test(femur.sample$Femur_sx)
shapiro.test(femur.sample$Femur_dx)
shapiro.test(femur.sample$FemurHead_sx)
shapiro.test(femur.sample$FemurHead_dx) 

str(femur.sample)
Femur_length<-femur.sample[,c(5,6)]
str(Femur_length)
Femur_length<-Femur_length[complete.cases(Femur_length),]
t.test(Femur_length$Femur_sx, y = Femur_length$Femur_dx)

Femur_head<-femur.sample[,c(7,8)]
str(Femur_head)
Femur_head<-Femur_head[complete.cases(Femur_head),]
wilcox.test(Femur_head$FemurHead_sx, y = Femur_head$FemurHead_dx)

#####################################################################################################
## 3) Descriptive statistics
sample <- read.table(file.choose(),header=T,row.names = 1,sep="\t")

## 3.1) FEM: maximum length of the femur##
#BOLOGNA FEMALES#
summary(sample$Femur_dx[1:31],na.rm = TRUE)
sd(sample$Femur_dx[1:31],na.rm = TRUE)
#YA (Young Adult)
summary(sample$Femur_dx[1:16],na.rm = TRUE)
sd(sample$Femur_dx[1:16],na.rm = TRUE)
#MA (Middle Adult)
summary(sample$Femur_dx[17:25],na.rm = TRUE)
sd(sample$Femur_dx[17:25],na.rm = TRUE)
#OA (Old Adult)
summary(sample$Femur_dx[26:31],na.rm = TRUE)
sd(sample$Femur_dx[26:31],na.rm = TRUE)

#BOLOGNA MALES#
summary(sample$Femur_dx[32:64],na.rm = TRUE)
sd(sample$Femur_dx[32:64],na.rm = TRUE)
#YA
summary(sample$Femur_dx[32:42],na.rm = TRUE)
sd(sample$Femur_dx[32:42],na.rm = TRUE)
#MA
summary(sample$Femur_dx[43:50],na.rm = TRUE)
sd(sample$Femur_dx[43:50],na.rm = TRUE)
#OA
summary(sample$Femur_dx[51:64],na.rm = TRUE)
sd(sample$Femur_dx[51:64],na.rm = TRUE)

#SASSARI FEMALES#
summary(sample$Femur_dx[65:96],na.rm = TRUE)
sd(sample$Femur_dx[65:96],na.rm = TRUE)
#YA
summary(sample$Femur_dx[65:79],na.rm = TRUE)
sd(sample$Femur_dx[65:79],na.rm = TRUE)
#MA
summary(sample$Femur_dx[80:89],na.rm = TRUE)
sd(sample$Femur_dx[80:89],na.rm = TRUE)
#OA
summary(sample$Femur_dx[90:96],na.rm = TRUE)
sd(sample$Femur_dx[90:96],na.rm = TRUE)

#SASSARI MALES#
summary(sample$Femur_dx[97:132],na.rm = TRUE)
sd(sample$Femur_dx[97:132],na.rm = TRUE)
#YA
summary(sample$Femur_dx[97:104],na.rm = TRUE)
sd(sample$Femur_dx[97:104],na.rm = TRUE)
#MA
summary(sample$Femur_dx[105:118],na.rm = TRUE)
sd(sample$Femur_dx[105:118],na.rm = TRUE)
#OA
summary(sample$Femur_dx[119:132],na.rm = TRUE)
sd(sample$Femur_dx[119:132],na.rm = TRUE)

## 3.2) FHSI: femoral head superoinferior breadth##
#BOLOGNA FEMALES#
summary(sample$FemurHead_dx[1:31],na.rm = TRUE)
sd(sample$FemurHead_dx[1:31],na.rm = TRUE)
#YA
summary(sample$FemurHead_dx[1:16],na.rm = TRUE)
sd(sample$FemurHead_dx[1:16],na.rm = TRUE)
#MA
summary(sample$FemurHead_dx[17:25],na.rm = TRUE)
sd(sample$FemurHead_dx[17:25],na.rm = TRUE)
#OA
summary(sample$FemurHead_dx[26:31],na.rm = TRUE)
sd(sample$FemurHead_dx[26:31],na.rm = TRUE)

#BOLOGNA MALES#
summary(sample$FemurHead_dx[32:64],na.rm = TRUE)
sd(sample$FemurHead_dx[32:64],na.rm = TRUE)
#YA
summary(sample$FemurHead_dx[32:42],na.rm = TRUE)
sd(sample$FemurHead_dx[32:42],na.rm = TRUE)
#MA
summary(sample$FemurHead_dx[43:50],na.rm = TRUE)
sd(sample$FemurHead_dx[43:50],na.rm = TRUE)
#OA
summary(sample$FemurHead_dx[51:64],na.rm = TRUE)
sd(sample$FemurHead_dx[51:64],na.rm = TRUE)

#SASSARI FEMALES#
summary(sample$FemurHead_dx[65:96],na.rm = TRUE)
sd(sample$FemurHead_dx[65:96],na.rm = TRUE)
#YA
summary(sample$FemurHead_dx[65:79],na.rm = TRUE)
sd(sample$FemurHead_dx[65:79],na.rm = TRUE)
#MA
summary(sample$FemurHead_dx[80:89],na.rm = TRUE)
sd(sample$FemurHead_dx[80:89],na.rm = TRUE)
#OA
summary(sample$FemurHead_dx[90:96],na.rm = TRUE)
sd(sample$FemurHead_dx[90:96],na.rm = TRUE)

#SASSARI MALES#
summary(sample$FemurHead_dx[97:132],na.rm = TRUE)
sd(sample$FemurHead_dx[97:132],na.rm = TRUE)
#YA
summary(sample$FemurHead_dx[97:104],na.rm = TRUE)
sd(sample$FemurHead_dx[97:104],na.rm = TRUE)
#MA
summary(sample$FemurHead_dx[105:118],na.rm = TRUE)
sd(sample$FemurHead_dx[105:118],na.rm = TRUE)
#OA
summary(sample$FemurHead_dx[119:132],na.rm = TRUE)
sd(sample$FemurHead_dx[119:132],na.rm = TRUE)

## 3.3) Stature (Trotter_Gleser), Body mass (Ruff) and Geometric_mean for each group
#BOLOGNA FEMALES#
summary(sample$Trotter_Gleser[1:31],na.rm = TRUE)
sd(sample$Trotter_Gleser[1:31],na.rm = TRUE)
summary(sample$Ruff[1:31],na.rm = TRUE)
sd(sample$Ruff[1:31],na.rm = TRUE)
summary(sample$Geometric_mean[1:31],na.rm = TRUE)
sd(sample$Geometric_mean[1:31],na.rm = TRUE)

#BOLOGNA MALES#
summary(sample$Trotter_Gleser[32:64],na.rm = TRUE)
sd(sample$Trotter_Gleser[32:64],na.rm = TRUE)
summary(sample$Ruff[32:64],na.rm = TRUE)
sd(sample$Ruff[32:64],na.rm = TRUE)
summary(sample$Geometric_mean[32:64],na.rm = TRUE)
sd(sample$Geometric_mean[32:64],na.rm = TRUE)

#SASSARI FEMALES#
summary(sample$Trotter_Gleser[65:96],na.rm = TRUE)
sd(sample$Trotter_Gleser[65:96],na.rm = TRUE)
summary(sample$Ruff[65:96],na.rm = TRUE)
sd(sample$Ruff[65:96],na.rm = TRUE)
summary(sample$Geometric_mean[65:96],na.rm = TRUE)
sd(sample$Geometric_mean[65:96],na.rm = TRUE)

#SASSARI MALES#
summary(sample$Trotter_Gleser[97:132],na.rm = TRUE)
sd(sample$Trotter_Gleser[97:132],na.rm = TRUE)
summary(sample$Ruff[97:132],na.rm = TRUE)
sd(sample$Ruff[97:132],na.rm = TRUE)
summary(sample$Geometric_mean[97:132],na.rm = TRUE)
sd(sample$Geometric_mean[97:132],na.rm = TRUE)

##########################################################################################################
### 4) #Shapiro-Wilk (normality test) and consequent Anova/Kruskal-Wallis for each variable and each group
sample <- read.table(file.choose(),header=T,row.names = 1,sep="\t")

## 4.1) FEM: maximum length of the femur##
#BOLOGNA FEMALES#
shapiro.test(sample$Femur_dx[1:31])
#YA
shapiro.test(sample$Femur_dx[1:16]) 
#MA
shapiro.test(sample$Femur_dx[17:25])
#OA
shapiro.test(sample$Femur_dx[26:31])
#BOLOGNA MALES#
shapiro.test(sample$Femur_dx[32:64])
#YA
shapiro.test(sample$Femur_dx[32:42])
#MA
shapiro.test(sample$Femur_dx[43:50])
#OA
shapiro.test(sample$Femur_dx[51:64])
#SASSARI FEMALES#
shapiro.test(sample$Femur_dx[65:96])
#YA
shapiro.test(sample$Femur_dx[65:79])
#MA
shapiro.test(sample$Femur_dx[80:89])
#OA
shapiro.test(sample$Femur_dx[90:96])
#SASSARI MALES#
shapiro.test(sample$Femur_dx[97:132])
#YA
shapiro.test(sample$Femur_dx[97:104])
#MA
shapiro.test(sample$Femur_dx[105:118])
#OA
shapiro.test(sample$Femur_dx[119:132])

##ANOVA/KRUSKAL
#BOLOGNA FEMALES
filtered_BOF <- subset(sample, Group_Sex == "BO_F")
table(filtered_BOF$Class)
kruskal.test(Femur_dx ~ Class, data = filtered_BOF)
pairwise.wilcox.test(filtered_BOF$Ruff,
                     filtered_BOF$Class,
                     p.adjust.method = "bonferroni")

#BOLOGNA MALES
filtered_BOM <- subset(sample, Group_Sex == "BO_M")
table(filtered_BOM$Class)
anova_BOM_age <- aov(Femur_dx ~ Class, data = filtered_BOM)
summary(anova_BOM_age)
TukeyHSD(anova_BOM_age)

#SASSARI FEMALES
filtered_SASF <- subset(sample, Group_Sex == "SAS_F")
table(filtered_SASF$Class)
anova_SASF_age <- aov(Femur_dx ~ Class, data = filtered_SASF)
summary(anova_SASF_age)
TukeyHSD(anova_SASF_age)

#SASSARI MALES
filtered_SASM <- subset(sample, Group_Sex == "SAS_M")
table(filtered_SASM$Class)
anova_SASM_age <- aov(Femur_dx ~ Class, data = filtered_SASM)
summary(anova_SASM_age)
TukeyHSD(anova_SASM_age)

## 4.2) FHSI: femoral head superoinferior breadth##
#BOLOGNA FEMALES#
shapiro.test(sample$FemurHead_dx[1:31]) 
#YA
shapiro.test(sample$FemurHead_dx[1:16]) 
#MA
shapiro.test(sample$FemurHead_dx[17:25])
#OA
shapiro.test(sample$FemurHead_dx[26:31])
#BOLOGNA MALES#
shapiro.test(sample$FemurHead_dx[32:64])
#YA
shapiro.test(sample$FemurHead_dx[32:42])
#MA
shapiro.test(sample$FemurHead_dx[43:50])
#OA
shapiro.test(sample$FemurHead_dx[51:64])
#SASSARI FEMALES#
shapiro.test(sample$FemurHead_dx[65:96]) 
#YA
shapiro.test(sample$FemurHead_dx[65:79])
#MA
shapiro.test(sample$FemurHead_dx[80:89])
#OA
shapiro.test(sample$FemurHead_dx[90:96]) 
#SASSARI MALES#
shapiro.test(sample$FemurHead_dx[97:132])
#YA
shapiro.test(sample$FemurHead_dx[97:104])
#MA
shapiro.test(sample$FemurHead_dx[105:118])
#OA
shapiro.test(sample$FemurHead_dx[119:132])

#ANOVA/KRUSKAL
#BOLOGNA FEMALES
filtered_BOF <- subset(sample, Group_Sex == "BO_F")
table(filtered_BOF$Class)

kruskal.test(FemurHead_dx ~ Class, data = filtered_BOF)
pairwise.wilcox.test(filtered_BOF$Ruff,
                     filtered_BOF$Class,
                     p.adjust.method = "bonferroni")

#BOLOGNA MALES
filtered_BOM <- subset(sample, Group_Sex == "BO_M")
table(filtered_BOM$Class)

kruskal.test(FemurHead_dx ~ Class, data = filtered_BOM)
pairwise.wilcox.test(filtered_BOM$Ruff,
                     filtered_BOM$Class,
                     p.adjust.method = "bonferroni")

#SASSARI FEMALES
filtered_SASF <- subset(sample, Group_Sex == "SAS_F")
table(filtered_SASF$Class)

kruskal.test(FemurHead_dx ~ Class, data = filtered_SASF)
pairwise.wilcox.test(filtered_SASF$Ruff,
                     filtered_SASF$Class,
                     p.adjust.method = "bonferroni")

#SASSARI MALES
filtered_SASM <- subset(sample, Group_Sex == "SAS_M")
table(filtered_SASM$Class)

anova_SASM_age <- aov(FemurHead_dx ~ Class, data = filtered_SASM)
summary(anova_SASM_age)
TukeyHSD(anova_SASM_age)

## 4.3) Stature (Trotter_Gleser), Body mass (Ruff) and Geometric_mean for each group
#BOLOGNA FEMALES
shapiro.test(sample$Trotter_Gleser[1:31]) 
shapiro.test(sample$Ruff[1:31]) 
shapiro.test(sample$Geometric_mean[1:31]) 
#BOLOGNA MALES
shapiro.test(sample$Trotter_Gleser[32:64]) 
shapiro.test(sample$Ruff[32:64]) 
shapiro.test(sample$Geometric_mean[32:64]) 
#SASSARI FEMALES
shapiro.test(sample$Trotter_Gleser[65:96]) 
shapiro.test(sample$Ruff[65:96]) 
shapiro.test(sample$Geometric_mean[65:96]) 
#SASSARI MALES
shapiro.test(sample$Trotter_Gleser[97:132]) 
shapiro.test(sample$Ruff[97:132]) 
shapiro.test(sample$Geometric_mean[97:132])

#ANOVA/KRUSKAL
#STATURE (Trotter_Gleser)#
filtered_sampleBO <- subset(sample, Group_Sex %in% c("BO_F", "BO_M"))
filtered_sampleBO$Group_Sex <- as.factor(filtered_sampleBO$Group_Sex)
anova_BOfemminevsmaschi <- aov(Trotter_Gleser ~ Group_Sex, data = filtered_sampleBO)
summary(anova_BOfemminevsmaschi)

filtered_sampleSAS <- subset(sample, Group_Sex %in% c("SAS_F", "SAS_M"))
filtered_sampleSAS$Group_Sex <- as.factor(filtered_sampleSAS$Group_Sex)
anova_SASfemminevsmaschi <- aov(Trotter_Gleser ~ Group_Sex, data = filtered_sampleSAS)
summary(anova_SASfemminevsmaschi)

filtered_sampleF <- subset(sample, Group_Sex %in% c("BO_F", "SAS_F"))
filtered_sampleF$Group_Sex <- as.factor(filtered_sampleF$Group_Sex)
anova_femminevsfemmine <- aov(Trotter_Gleser ~ Group_Sex, data = filtered_sampleF)
summary(anova_femminevsfemmine)

filtered_sampleM <- subset(sample, Group_Sex %in% c("BO_M", "SAS_M"))
filtered_sampleM$Group_Sex <- as.factor(filtered_sampleM$Group_Sex)
anova_maschivsmaschi <- aov(Trotter_Gleser ~ Group_Sex, data = filtered_sampleM)
summary(anova_maschivsmaschi)

#BODY MASS (Ruff)#
filtered_sampleBO <- subset(sample, Group_Sex %in% c("BO_M", "BO_F"))
filtered_sampleBO$Group_Sex <- as.factor(filtered_sampleBO$Group_Sex)
anova_BOfemminevsmaschi <- aov(Ruff ~ Group_Sex, data = filtered_sampleBO)
summary(anova_BOfemminevsmaschi)

filtered_data <- subset(sample, Group_Sex %in% c("SAS_M", "SAS_F"))
kruskal_test <- kruskal.test(Ruff ~ Group_Sex, data = filtered_data)
print(kruskal_test)

filtered_data <- subset(sample, Group_Sex %in% c("BO_F", "SAS_F"))
kruskal_test <- kruskal.test(Ruff ~ Group_Sex, data = filtered_data)
print(kruskal_test)

filtered_sampleM <- subset(sample, Group_Sex %in% c("BO_M", "SAS_M"))
filtered_sampleM$Group_Sex <- as.factor(filtered_sampleM$Group_Sex)
anova_maschivsmaschi <- aov(Ruff ~ Group_Sex, data = filtered_sampleM)
summary(anova_maschivsmaschi)

#GEOMETRIC MEAN#
filtered_data <- subset(sample, Group_Sex %in% c("BO_M", "BO_F"))
kruskal_test <- kruskal.test(Geometric_mean ~ Group_Sex, data = filtered_data)
print(kruskal_test)

filtered_sampleSAS <- subset(sample, Group_Sex %in% c("SAS_M", "SAS_F"))
filtered_sampleSAS$Group_Sex <- as.factor(filtered_sampleSAS$Group_Sex)
anova_SASfemminevsmaschi <- aov(Geometric_mean ~ Group_Sex, data = filtered_sampleSAS)
summary(anova_SASfemminevsmaschi)

filtered_data <- subset(sample, Group_Sex %in% c("BO_F", "SAS_F"))
kruskal_test <- kruskal.test(Geometric_mean ~ Group_Sex, data = filtered_data)
print(kruskal_test)

filtered_sampleM <- subset(sample, Group_Sex %in% c("BO_M", "SAS_M"))
filtered_sampleM$Group_Sex <- as.factor(filtered_sampleM$Group_Sex)
anova_maschivsmaschi <- aov(Geometric_mean ~ Group_Sex, data = filtered_sampleM)
summary(anova_maschivsmaschi)

#SHAPIRO TEST -  BOLOGNA AGE CLASS#
shapiro.test(sample$Trotter_Gleser[1:16]) 
shapiro.test(sample$Ruff[1:16]) 
shapiro.test(sample$Geometric_mean[1:16]) 
shapiro.test(sample$Trotter_Gleser[17:25]) 
shapiro.test(sample$Ruff[17:25]) 
shapiro.test(sample$Geometric_mean[17:25]) 
shapiro.test(sample$Trotter_Gleser[26:31]) 
shapiro.test(sample$Ruff[26:31]) 
shapiro.test(sample$Geometric_mean[26:31]) 
shapiro.test(sample$Trotter_Gleser[32:42]) 
shapiro.test(sample$Ruff[32:42]) 
shapiro.test(sample$Geometric_mean[32:42]) 
shapiro.test(sample$Trotter_Gleser[43:50]) 
shapiro.test(sample$Ruff[43:50]) 
shapiro.test(sample$Geometric_mean[43:50]) 
shapiro.test(sample$Trotter_Gleser[51:64]) 
shapiro.test(sample$Ruff[51:64]) 
shapiro.test(sample$Geometric_mean[51:64]) 

#BOLOGNA FEMALES
filtered_BOF <- subset(sample, Group_Sex == "BO_F")

table(filtered_BOF$Class)
anova_BOF_age <- aov(Trotter_Gleser ~ Class, data = filtered_BOF)
summary(anova_BOF_age)
TukeyHSD(anova_BOF_age)

kruskal.test(Ruff ~ Class, data = filtered_BOF)
pairwise.wilcox.test(filtered_BOF$Ruff,
                     filtered_BOF$Class,
                     p.adjust.method = "bonferroni")

kruskal.test(Geometric_mean ~ Class, data = filtered_BOF)
pairwise.wilcox.test(filtered_BOF$Geometric_mean,
                     filtered_BOF$Class,
                     p.adjust.method = "bonferroni")

filtered_BOF <- subset(sample, Group_Sex == "BO_F")
table(filtered_BOF$Class)
anova_BOF_age <- aov(Trotter_Gleser ~ Class, data = filtered_BOF)
summary(anova_BOF_age)
TukeyHSD(anova_BOF_age)

#BOLOGNA MALES
filtered_BOM <- subset(sample, Group_Sex == "BO_M")
table(filtered_BOM$Class)

anova_BOM_age <- aov(Trotter_Gleser ~ Class, data = filtered_BOM)
summary(anova_BOM_age)
TukeyHSD(anova_BOM_age)

anova_BOM_age <- aov(Ruff ~ Class, data = filtered_BOM)
summary(anova_BOM_age)
TukeyHSD(anova_BOM_age)

anova_BOM_age <- aov(Geometric_mean ~ Class, data = filtered_BOM)
summary(anova_BOM_age)
TukeyHSD(anova_BOM_age)

#SHAPIRO TEST -  SASSARI AGE CLASS#
shapiro.test(sample$Trotter_Gleser[65:79]) 
shapiro.test(sample$Ruff[65:79]) 
shapiro.test(sample$Geometric_mean[65:79]) 
shapiro.test(sample$Trotter_Gleser[80:89]) 
shapiro.test(sample$Ruff[80:89]) 
shapiro.test(sample$Geometric_mean[80:89]) 
shapiro.test(sample$Trotter_Gleser[90:96]) 
shapiro.test(sample$Ruff[90:96]) 
shapiro.test(sample$Geometric_mean[90:96]) 
shapiro.test(sample$Trotter_Gleser[97:104]) 
shapiro.test(sample$Ruff[97:104]) 
shapiro.test(sample$Geometric_mean[97:104]) 
shapiro.test(sample$Trotter_Gleser[105:118]) 
shapiro.test(sample$Ruff[105:118]) 
shapiro.test(sample$Geometric_mean[105:118]) 
shapiro.test(sample$Trotter_Gleser[119:132]) 
shapiro.test(sample$Ruff[119:132]) 
shapiro.test(sample$Geometric_mean[119:132]) 

#SASSARI FEMALES
filtered_SASF <- subset(sample, Group_Sex == "SAS_F")
table(filtered_SASF$Class)

anova_SASF_age <- aov(Trotter_Gleser ~ Class, data = filtered_SASF)
summary(anova_SASF_age)
TukeyHSD(anova_SASF_age)

anova_SASF_age <- aov(Ruff ~ Class, data = filtered_SASF)
summary(anova_SASF_age)
TukeyHSD(anova_SASF_age)

anova_SASF_age <- aov(Geometric_mean ~ Class, data = filtered_SASF)
summary(anova_SASF_age)
TukeyHSD(anova_SASF_age)

#SASSARI MALES
filtered_SASM <- subset(sample, Group_Sex == "SAS_M")
table(filtered_SASM$Class)

anova_SASM_age <- aov(Trotter_Gleser ~ Class, data = filtered_SASM)
summary(anova_SASM_age)
TukeyHSD(anova_SASM_age)

kruskal.test(Ruff ~ Class, data = filtered_SASM)
pairwise.wilcox.test(filtered_SASM$Ruff,
                     filtered_SASM$Class,
                     p.adjust.method = "bonferroni")

anova_SASM_age <- aov(Geometric_mean ~ Class, data = filtered_SASM)
summary(anova_SASM_age)
TukeyHSD(anova_SASM_age)

####################################################################################
## 5) STANDARDIZED MAJOR AXIS ANALYSIS - SMA and PAIRWISE TESTS ##
library(smatr)
library(dplyr)
sample <- read.table(file.choose(),header=T,row.names = 1,sep="\t")

#M01#
vars <- c("M01_dx", "Geometric_mean", "M01_log", "Group_Sex")
sampleM01_comp <- sample[complete.cases(sample[, vars]), ]

sampleM01_clean <- sampleM01_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M01_log, na.rm = TRUE),
    sd_val = sd(M01_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M01_log < lower_bound | M01_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM01_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M01_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM01_clean %>%
  filter(!is_outlier)

sma_M01 <- sma(M01_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M01$groupsummary
summary(sma_M01)
sma_M01_shift <- sma(M01_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M01_shift

plot(sma_M01, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M01")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
#Common slope and elevation
  sma_pair <- sma(
    M01_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
#Shift along the common axis
  sma_pair_shift <- sma(
    M01_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M04#
vars <- c("M04_dx", "Geometric_mean", "M04_log", "Group_Sex")
sampleM04_comp <- sample[complete.cases(sample[, vars]), ]

sampleM04_clean <- sampleM04_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M04_log, na.rm = TRUE),
    sd_val = sd(M04_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M04_log < lower_bound | M04_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM04_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M04_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM04_clean %>%
  filter(!is_outlier)

sma_M04 <- sma(M04_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M04$groupsummary
summary(sma_M04)
sma_M04_shift <- sma(M04_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M04_shift

plot(sma_M04, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M04")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M04_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M04_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M12#
vars <- c("M12_dx", "Geometric_mean", "M12_log", "Group_Sex")
sampleM12_comp <- sample[complete.cases(sample[, vars]), ]

sampleM12_clean <- sampleM12_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M12_log, na.rm = TRUE),
    sd_val = sd(M12_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M12_log < lower_bound | M12_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM12_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M12_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM12_clean %>%
  filter(!is_outlier)

sma_M12 <- sma(M12_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M12$groupsummary
summary(sma_M12)
sma_M12_shift <- sma(M12_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M12_shift

plot(sma_M12, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M12")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M12_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M12_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M13#
vars <- c("M13_dx", "Geometric_mean", "M13_log", "Group_Sex")
sampleM13_comp <- sample[complete.cases(sample[, vars]), ]

sampleM13_clean <- sampleM13_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M13_log, na.rm = TRUE),
    sd_val = sd(M13_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M13_log < lower_bound | M13_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM13_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M13_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM13_clean %>%
  filter(!is_outlier)

sma_M13 <- sma(M13_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M13$groupsummary
summary(sma_M13)
sma_M13_shift <- sma(M13_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M13_shift

plot(sma_M13, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M13")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M13_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M13_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M14#
vars <- c("M14_dx", "Geometric_mean", "M14_log", "Group_Sex")
sampleM14_comp <- sample[complete.cases(sample[, vars]), ]

sampleM14_clean <- sampleM14_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M14_log, na.rm = TRUE),
    sd_val = sd(M14_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M14_log < lower_bound | M14_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM14_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M14_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM14_clean %>%
  filter(!is_outlier)

sma_M14 <- sma(M14_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M14$groupsummary
summary(sma_M14)
sma_M14_shift <- sma(M14_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M14_shift

plot(sma_M14, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M14")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M14_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M14_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M14_1#
vars <- c("M14_1_dx", "Geometric_mean", "M14_1_log", "Group_Sex")
sampleM14_1_comp <- sample[complete.cases(sample[, vars]), ]

sampleM14_1_clean <- sampleM14_1_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M14_1_log, na.rm = TRUE),
    sd_val = sd(M14_1_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M14_1_log < lower_bound | M14_1_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM14_1_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M14_1_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM14_1_clean %>%
  filter(!is_outlier)

sma_M14_1 <- sma(M14_1_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M14_1$groupsummary
summary(sma_M14_1)
sma_M14_1_shift <- sma(M14_1_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M14_1_shift

plot(sma_M14_1, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M14_1")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M14_1_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M14_1_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M15_1#
vars <- c("M15_1._dx", "Geometric_mean", "M15_1._log", "Group_Sex")
sampleM15_1._comp <- sample[complete.cases(sample[, vars]), ]

sampleM15_1._clean <- sampleM15_1._comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M15_1._log, na.rm = TRUE),
    sd_val = sd(M15_1._log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M15_1._log < lower_bound | M15_1._log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM15_1._clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M15_1._log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM15_1._clean %>%
  filter(!is_outlier)

sma_M15_1. <- sma(M15_1._dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M15_1.$groupsummary
summary(sma_M15_1.)
sma_M15_1._shift <- sma(M15_1._dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M15_1._shift

plot(sma_M15_1., col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M15_1.")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M15_1._dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M15_1._dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M15_a#
vars <- c("M15_a._dx", "Geometric_mean", "M15_a._log", "Group_Sex")
sampleM15_a._comp <- sample[complete.cases(sample[, vars]), ]

sampleM15_a._clean <- sampleM15_a._comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M15_a._log, na.rm = TRUE),
    sd_val = sd(M15_a._log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M15_a._log < lower_bound | M15_a._log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM15_a._clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M15_a._log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM15_a._clean %>%
  filter(!is_outlier)

sma_M15_a. <- sma(M15_a._dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M15_a.$groupsummary
summary(sma_M15_a.)
sma_M15_a._shift <- sma(M15_a._dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M15_a._shift

plot(sma_M15_a., col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M15_a.")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M15_a._dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M15_a._dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M17_a#
vars <- c("M17_a_dx", "Geometric_mean", "M17_a_log", "Group_Sex")
sampleM17_a_comp <- sample[complete.cases(sample[, vars]), ]

sampleM17_a_clean <- sampleM17_a_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M17_a_log, na.rm = TRUE),
    sd_val = sd(M17_a_log, na.rm = TRUE),
    lower_bound = mean_val - 1.96 * sd_val,
    upper_bound = mean_val + 1.96 * sd_val,
    is_outlier = M17_a_log < lower_bound | M17_a_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM17_a_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M17_a_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM17_a_clean %>%
  filter(!is_outlier)

sma_M17_a <- sma(M17_a_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M17_a$groupsummary
summary(sma_M17_a)
sma_M17_a_shift <- sma(M17_a_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M17_a_shift

plot(sma_M17_a, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M17_a")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M17_a_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M17_a_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M18#
vars <- c("M18_dx", "Geometric_mean", "M18_log", "Group_Sex")
sampleM18_comp <- sample[complete.cases(sample[, vars]), ]

sampleM18_clean <- sampleM18_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M18_log, na.rm = TRUE),
    sd_val = sd(M18_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M18_log < lower_bound | M18_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM18_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M18_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM18_clean %>%
  filter(!is_outlier)

sma_M18 <- sma(M18_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M18$groupsummary
summary(sma_M18)
sma_M18_shift <- sma(M18_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M18_shift

plot(sma_M18, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M18")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M18_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M18_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M20#
vars <- c("M20_dx", "Geometric_mean", "M20_log", "Group_Sex")
sampleM20_comp <- sample[complete.cases(sample[, vars]), ]

sampleM20_clean <- sampleM20_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M20_log, na.rm = TRUE),
    sd_val = sd(M20_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M20_log < lower_bound | M20_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM20_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M20_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM20_clean %>%
  filter(!is_outlier)

sma_M20 <- sma(M20_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M20$groupsummary
summary(sma_M20)
sma_M20_shift <- sma(M20_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M20_shift

plot(sma_M20, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M20")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M20_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M20_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M21#
vars <- c("M21_dx", "Geometric_mean", "M21_log", "Group_Sex")
sampleM21_comp <- sample[complete.cases(sample[, vars]), ]

sampleM21_clean <- sampleM21_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M21_log, na.rm = TRUE),
    sd_val = sd(M21_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M21_log < lower_bound | M21_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM21_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M21_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM21_clean %>%
  filter(!is_outlier)

sma_M21 <- sma(M21_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M21$groupsummary
summary(sma_M21)
sma_M21_shift <- sma(M21_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M21_shift

plot(sma_M21, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M21")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M21_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M21_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#M22#
vars <- c("M22_dx", "Geometric_mean", "M22_log", "Group_Sex")
sampleM22_comp <- sample[complete.cases(sample[, vars]), ]

sampleM22_clean <- sampleM22_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(M22_log, na.rm = TRUE),
    sd_val = sd(M22_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = M22_log < lower_bound | M22_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleM22_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, M22_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleM22_clean %>%
  filter(!is_outlier)

sma_M22 <- sma(M22_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_M22$groupsummary
summary(sma_M22)
sma_M22_shift <- sma(M22_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_M22_shift

plot(sma_M22, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "M22")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    M22_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    M22_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#PUM#
vars <- c("PUM_dx", "Geometric_mean", "PUM_log", "Group_Sex")
samplePUM_comp <- sample[complete.cases(sample[, vars]), ]

samplePUM_clean <- samplePUM_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(PUM_log, na.rm = TRUE),
    sd_val = sd(PUM_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = PUM_log < lower_bound | PUM_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- samplePUM_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, PUM_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- samplePUM_clean %>%
  filter(!is_outlier)

sma_PUM <- sma(PUM_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_PUM$groupsummary
summary(sma_PUM)
sma_PUM_shift <- sma(PUM_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_PUM_shift

plot(sma_PUM, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "PUM")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    PUM_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    PUM_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#ISM#
vars <- c("ISM_dx", "Geometric_mean", "ISM_log", "Group_Sex")
sampleISM_comp <- sample[complete.cases(sample[, vars]), ]

sampleISM_clean <- sampleISM_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(ISM_log, na.rm = TRUE),
    sd_val = sd(ISM_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = ISM_log < lower_bound | ISM_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleISM_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, ISM_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleISM_clean %>%
  filter(!is_outlier)

sma_ISM <- sma(ISM_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_ISM$groupsummary
summary(sma_ISM)
sma_ISM_shift <- sma(ISM_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_ISM_shift

plot(sma_ISM, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "ISM")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    ISM_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    ISM_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#SPU#
vars <- c("SPU_dx", "Geometric_mean", "SPU_log", "Group_Sex")
sampleSPU_comp <- sample[complete.cases(sample[, vars]), ]

sampleSPU_clean <- sampleSPU_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(SPU_log, na.rm = TRUE),
    sd_val = sd(SPU_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = SPU_log < lower_bound | SPU_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleSPU_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, SPU_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleSPU_clean %>%
  filter(!is_outlier)

sma_SPU <- sma(SPU_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_SPU$groupsummary
summary(sma_SPU)
sma_SPU_shift <- sma(SPU_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_SPU_shift

plot(sma_SPU, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "SPU")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    SPU_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    SPU_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#ISMM#
vars <- c("ISMM_dx", "Geometric_mean", "ISMM_log", "Group_Sex")
sampleISMM_comp <- sample[complete.cases(sample[, vars]), ]

sampleISMM_clean <- sampleISMM_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(ISMM_log, na.rm = TRUE),
    sd_val = sd(ISMM_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = ISMM_log < lower_bound | ISMM_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleISMM_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, ISMM_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleISMM_clean %>%
  filter(!is_outlier)

sma_ISMM <- sma(ISMM_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_ISMM$groupsummary
summary(sma_ISMM)
sma_ISMM_shift <- sma(ISMM_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_ISMM_shift

plot(sma_ISMM, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "ISMM")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    ISMM_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    ISMM_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#SS#
vars <- c("SS_dx", "Geometric_mean", "SS_log", "Group_Sex")
sampleSS_comp <- sample[complete.cases(sample[, vars]), ]

sampleSS_clean <- sampleSS_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(SS_log, na.rm = TRUE),
    sd_val = sd(SS_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = SS_log < lower_bound | SS_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleSS_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, SS_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleSS_clean %>%
  filter(!is_outlier)

sma_SS <- sma(SS_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_SS$groupsummary
summary(sma_SS)
sma_SS_shift <- sma(SS_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_SS_shift

plot(sma_SS, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "SS")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    SS_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    SS_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}

#SA#
vars <- c("SA_dx", "Geometric_mean", "SA_log", "Group_Sex")
sampleSA_comp <- sample[complete.cases(sample[, vars]), ]

sampleSA_clean <- sampleSA_comp %>%
  group_by(Group_Sex) %>%
  mutate(
    mean_val = mean(SA_log, na.rm = TRUE),
    sd_val = sd(SA_log, na.rm = TRUE),
    lower_bound = mean_val - 2.5 * sd_val,
    upper_bound = mean_val + 2.5 * sd_val,
    is_outlier = SA_log < lower_bound | SA_log > upper_bound
  ) %>%
  ungroup()
outliers_removed <- sampleSA_clean %>%
  filter(is_outlier) %>%
  select(Group_Sex, SA_log)
print(outliers_removed)
cat("Total outliers removed:", nrow(outliers_removed), "\n")
cleaned_data <- sampleSA_clean %>%
  filter(!is_outlier)

sma_SA <- sma(SA_dx ~ Geometric_mean + Group_Sex, log = "xy", data = cleaned_data, slope.test = 1)
sma_SA$groupsummary
summary(sma_SA)
sma_SA_shift <- sma(SA_dx ~ Geometric_mean + Group_Sex, log = "xy", type = "shift", data = cleaned_data)
sma_SA_shift

plot(sma_SA, col = c("red", "blue", "orange", "green"), pch = 16, cex = 2, xlab = "Geometric Mean", ylab = "SA")

#Pairwise
comparisons <- list(
  BO_F_vs_BO_M = c("BO_F", "BO_M"),
  SAS_F_vs_SAS_M = c("SAS_F", "SAS_M"),
  BO_F_vs_SAS_F = c("BO_F", "SAS_F"),
  BO_M_vs_SAS_M = c("BO_M", "SAS_M")
)

for (comparison_name in names(comparisons)) {
  
  groups <- comparisons[[comparison_name]]
  
  pair_data <- cleaned_data %>%
    filter(Group_Sex %in% groups) %>%
    droplevels()
  
  cat("\n\n========================================\n")
  cat(comparison_name, "\n")
  cat("========================================\n\n")
  
  #Common slope and elevation
  sma_pair <- sma(
    SA_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    data = pair_data
  )
  
  print(summary(sma_pair))
  
  #Shift along the common axis
  sma_pair_shift <- sma(
    SA_dx ~ Geometric_mean + Group_Sex,
    log = "xy",
    type = "shift",
    data = pair_data
  )
  
  cat("\n--- Shift along common axis ---\n")
  print(sma_pair_shift)
}
