require(readxl)
require(easyCODA)
require(car)
require(emmeans)

### Brewers’ grains dataset

BGds <- as.data.frame(read_excel("~/Desktop/Excel datasets.xlsx",sheet = "Brewer's grains", skip = 1))
AGV_prod_BG <- BGds[,1:7]
AGV_prop_BG <- BGds[,c(1,9:14)]

#(1) Subcomposition

Total_C234_BG <- AGV_prod_BG$PAc + AGV_prod_BG$PPr + AGV_prod_BG$PBut
C234_Ac <- AGV_prod_BG$PAc/Total_C234_BG *100
C234_Pr <- AGV_prod_BG$PPr/Total_C234_BG *100
C234_Bu <- AGV_prod_BG$PBut/Total_C234_BG *100

C234_BG <- as.data.frame(cbind(C234_Ac,C234_Pr,C234_Bu))
type_BG <- AGV_prod_BG$Type

options(contrasts = c("contr.sum", "contr.poly"))
p.values_AGV_mol_3C_BG <- data.frame(matrix(,nrow=ncol(C234_BG),ncol = 5))
colnames(p.values_AGV_mol_3C_BG) <- c("Variable","p.values","CRAFT","INDUSTRIAL","SEM")

for (i in 1:3){
  sum <- Anova(model <- lm(C234_BG[,i] ~ type_BG), type = 3)
  SEM <- sqrt((sum$`Sum Sq`[3]/sum$Df[3])/5)
  means <- summary(emmeans(model, ~ type_BG))
  p.values_AGV_mol_3C_BG[i,3:4] <- round(means$emmean,3)
  p.values_AGV_mol_3C_BG$Variable[i] <- colnames(C234_BG)[i]
  p.values_AGV_mol_3C_BG$p.values[i] <- round(sum[2,4],3)
  p.values_AGV_mol_3C_BG$SEM[i] <- round(SEM,4)
}

#(2) Molar proportions

options(contrasts = c("contr.sum", "contr.poly"))
p.values_AGV_prop_BG <- data.frame(matrix(,nrow=ncol(AGV_prop_BG)-1,ncol = 5))
colnames(p.values_AGV_prop_BG) <- c("Variable","p.values","CRAFT","INDUSTRIAL","SEM")

for (i in 2:ncol(AGV_prop_BG)){
  sum <- Anova(model <- lm(AGV_prop_BG[,i] ~ type_BG), type = 3)
  SEM <- sqrt((sum$`Sum Sq`[3]/sum$Df[3])/5)
  means <- summary(emmeans(model, ~ type_BG))
  p.values_AGV_prop_BG[i-1,3:4] <- round(means$emmean,3)
  p.values_AGV_prop_BG$Variable[i-1] <- colnames(AGV_prop_BG)[i]
  p.values_AGV_prop_BG$p.values[i-1] <- round(sum[2,4],3)
  p.values_AGV_prop_BG$SEM[i-1] <- round(SEM,4)
}

#(3) Production

options(contrasts = c("contr.sum", "contr.poly"))
p.values_AGV_prod_BG <- data.frame(matrix(,nrow=ncol(AGV_prod_BG)-1,ncol = 5))
colnames(p.values_AGV_prod_BG) <- c("Variable","p.values","CRAFT","INDUSTRIAL","SEM")

for (i in 2:ncol(AGV_prod_BG)){
  sum <- Anova(model <- lm(AGV_prod_BG[,i] ~ type_BG), type = 3)
  SEM <- sqrt((sum$`Sum Sq`[3]/sum$Df[3])/5)
  means <- summary(emmeans(model, ~ type_BG))
  p.values_AGV_prod_BG[i-1,3:4] <- round(means$emmean,3)
  p.values_AGV_prod_BG$Variable[i-1] <- colnames(AGV_prod_BG)[i]
  p.values_AGV_prod_BG$p.values[i-1] <- round(sum[2,4],3)
  p.values_AGV_prod_BG$SEM[i-1] <- round(SEM,4)
}

#(4) CLR

AGV_closed_BG_1 <-CLOSE(AGV_prod_BG[,2:7])
AGV_CLR_1 <- CLR(AGV_closed_BG_1,weight=FALSE)$LR; AGV_CLR_1[1,]
AGV_closed_BG <-CLOSE(AGV_prop_BG[,2:7])
AGV_CLR <- CLR(AGV_closed_BG,weight=FALSE)$LR; AGV_CLR[1,] #equal

options(contrasts = c("contr.sum", "contr.poly"))
p.values_AGV_CLR_BG <- data.frame(matrix(,nrow=ncol(AGV_CLR),ncol = 5)) 
colnames(p.values_AGV_CLR_BG) <- c("Variable","p.values","CRAFT","INDUSTRIAL","SEM")
for (i in 1:ncol(AGV_CLR)){
  sum <- Anova(model <- lm(AGV_CLR[,i] ~ type_BG), type = 3)
  SEM <- sqrt((sum$`Sum Sq`[3]/sum$Df[3])/5)
  means <- summary(emmeans(model, ~ type_BG))
  p.values_AGV_CLR_BG[i,3:4] <- round(means$emmean,2)
  p.values_AGV_CLR_BG$Variable[i] <- colnames(AGV_CLR)[i]
  p.values_AGV_CLR_BG$p.values[i] <- round(sum[2,4],3)
  p.values_AGV_CLR_BG$SEM[i] <- round(SEM,4)
}

#### Almond hulls silages dataset

AHSds <- as.data.frame(read_excel("~/Desktop/Excel datasets.xlsx",sheet = "Almond hulls silages", skip = 1))

trat_AHS <- factor(AHSds$Treatment,levels = c("Control","AC1","AC2","BAC"))
REM <- 100 - apply(AHSds[,3:7],1,sum)
HEM <- AHSds$NDF - AHSds$ADF
CEL <- AHSds$ADF - AHSds$ADL

AHSds_100 <- AHSds[,3:6]
AHSds_100$HEM <- HEM
AHSds_100$CEL <- CEL
AHSds_100$ADL <- AHSds$ADL
AHSds_100$REM <- REM
AHSds_100 <- AHSds_100*10

#(1) DM basis

options(contrasts = c("contr.sum", "contr.poly"))
p.values_AH <- data.frame(matrix(,nrow=ncol(AHSds_100),ncol = 7))
colnames(p.values_AH) <- c("Variable","p.values","CON","AC1","AC2","BAC","SEM")

for (i in 1:ncol(AHSds_100)){
  sum <- Anova(model <- lm(AHSds_100[,i] ~ trat_AHS), type = 3)
  SEM <- sqrt((sum$`Sum Sq`[3]/sum$Df[3])/5)
  means <- summary(emmeans(model, ~ trat_AHS))
  p.values_AH[i,3:6] <- round(means$emmean,2)
  p.values_AH$Variable[i] <- colnames(AHSds_100)[i]
  p.values_AH$p.values[i] <- round(sum[2,4],3)
  p.values_AH$SEM[i] <- round(SEM,4)
}

#(2) ALR

AH.FINDALR <- FINDALR(CLOSE(AHSds_100), weight=FALSE)
AH.FINDALR$procrust.cor
AH.FINDALR$var.log
colnames(AHSds_100)

options(contrasts = c("contr.sum", "contr.poly"))
p.values_AH_ALR <- data.frame(matrix(,nrow=ncol(AHSds_100)-1,ncol = 7)) 
colnames(p.values_AH_ALR) <- c("Variable","p.values","CON","AC1","AC2","BAC","SEM")

for (i in 2:ncol(AHSds_100)){
  log.ASH <- log(AHSds_100[,i]/AHSds_100$Ashes)
  sum <- Anova(model <- lm(log.ASH ~ trat_AHS), type = 3)
  SEM <- sqrt((sum$`Sum Sq`[3]/sum$Df[3])/5)
  means <- summary(emmeans(model, ~ trat_AHS))
  p.values_AH_ALR[i-1,3:6] <- round(means$emmean,2)
  p.values_AH_ALR$Variable[i-1] <- colnames(AHSds_100)[i]
  p.values_AH_ALR$p.values[i-1] <- round(sum[2,4],3)
  p.values_AH_ALR$SEM[i-1] <- round(SEM,4)
}

#### Agroindustrial by-products dataset

ABYPds <- as.data.frame(read_excel("~/Desktop/Excel datasets.xlsx",sheet = "Agroindustrial by-products", skip = 1))

ABYPds$Ashes <- 100 - ABYPds$OM
HEM <- ABYPds$NDF - ABYPds$ADF
CEL <- ABYPds$ADF - ABYPds$ADL
ABYPds$HEM <- HEM
ABYPds$CEL <- CEL

ABYP_QC <- ABYPds[,c(5:8,10:12)]
ABYP_AGV <- ABYPds[,9]

#(1) DM basis
cor(ABYP_QC,ABYP_AGV)

#(2) LR
ABYPds.PLR <- LR(CLOSE(ABYP_QC), weight=FALSE)$LR 
cor(ABYPds.PLR,ABYP_AGV)

#### Fatty acids dataset

FAds <- as.data.frame(read_excel("~/Desktop/Excel datasets.xlsx",sheet = "Fatty acids", skip = 1))

SC <- FAds[,2:24]
IM <- FAds[,25:48]

SC_close <- CLOSE(SC) 
SC_close <- na.omit(SC_close)
IM_close <- CLOSE(IM)

df_SC <- data.frame(matrix(,nrow = length(colnames(SC_close)), ncol = 3))
colnames(df_SC) <- c("Variable","Corr.","Var.")

AGV.FINDALR <- FINDALR(SC_close, weight=FALSE)
df_SC[,2] <- AGV.FINDALR$procrust.cor
df_SC[,3] <- AGV.FINDALR$var.log
df_SC[,1] <- colnames(SC_close)

df_IM <- data.frame(matrix(,nrow = length(colnames(IM_close)), ncol = 3))
colnames(df_IM) <- c("Variable","Corr.","Var.")

AGV.FINDALR <- FINDALR(IM_close, weight=FALSE)
df_IM[,2] <- AGV.FINDALR$procrust.cor
df_IM[,3] <- AGV.FINDALR$var.log
df_IM[,1] <- colnames(IM_close)