#------Packages------
library(PLNmodels)
library(ggplot2)
library(corrplot)
library(tidyr)
library(dplyr)
library(factoextra)
library(ade4)
library(FactoMineR)
library(Metrics)
library(stringr)
library(rlang)
library(ggpubr)
library(pals)
library(ggrepel)
library(igraph)
library(ggraph)

#Put the data into a folder named data

#------INFERRING ASSOCIATIONS WITHOUT SPECIES TRAITS (8 species community)------
#Loading data------
terrain8species_eco<-read.csv2("data/FF8_abund_eco.csv",dec=".")
nom_mouches<-read.csv2("data/fly_labels.csv")

#Famd on ecological covariates------
dat_eco=terrain8species_eco[,c(13,14,22,25:31)]
dat_eco$mois=as.factor(dat_eco$mois)
dat_eco$annee=as.factor(dat_eco$annee)

res.famd=FAMD(dat_eco, graph=F,ncp = 10)
ind<-res.famd$ind$coord
ind=ind[,c(1:10)]
colnames(ind)=c("FAMD.1","FAMD.2","FAMD.3","FAMD.4","FAMD.5","FAMD.6","FAMD.7","FAMD.8","FAMD.9","FAMD.10")

#Visualisation of FAMD results (used in Fig 4 and S2)
fviz_screeplot(res.famd)
fviz_famd_var(res.famd, "quanti.var", col.var = "cos2",
							gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
							repel = TRUE)
 
fviz_famd_var(res.famd, "quali.var", col.var = "cos2",
							gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
							repel = TRUE)
 
fviz_contrib (res.famd, "var", axes = 1)
fviz_contrib (res.famd, "var", axes = 2)
fviz_contrib (res.famd, "var", axes = c(3:10))

#Preparing data for PLN------
donnees=cbind(terrain8species_eco[,c(1,2,4:12)],ind)
rn_abondance<-as.matrix(donnees[,4:11])
colnames(rn_abondance)=nom_mouches$latin[match(names(donnees)[4:11], nom_mouches$terrain)]
rownames(rn_abondance)=donnees$numech
rn_covariate<-donnees[,c(2,3,12:21)]
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance

#PLN modelling------
#Model inference
myPLN0_base <- PLN(Abundance ~ 1, data = rn_pln)
myPLN0_poids <- PLN(Abundance ~ 1 + offset(log(poids)), data = rn_pln)

myPLN_plante <- PLN(Abundance ~ 1 + factor(plantesc) + offset(log(poids)), data = rn_pln)
myPLN_plante_diag <- PLN(Abundance ~ 1 + factor(plantesc) + offset(log(poids)), data = rn_pln, control = list(covariance = "diagonal"))

myPLN_famd <- PLN(Abundance ~ 1 +  FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)
myPLN_famd_diag <- PLN(Abundance ~ 1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10+ offset(log(poids)), data = rn_pln, control = list(covariance = "diagonal"))

myPLN_pl_famd <- PLN(Abundance ~ 1 + factor(plantesc) + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)
myPLN_pl_famd_diag <- PLN(Abundance ~ 1 + factor(plantesc) + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln, control = list(covariance = "diagonal"))

#Residual variance covariance matrices (Fig 2)
col<-brewer.rdgy(10) 
par(mfrow=c(2,2))

cp0_poids=corrplot(sigma(myPLN0_poids), is.corr = F, method="circle", addgrid.col = "lightgray", col=col,tl.col="white", tl.cex=0.0000001, cl.lim=c(-21,21),cl.pos="b")

cp_eco=corrplot(sigma(myPLN_famd),  is.corr = F, method="circle", addgrid.col = "lightgray", col=col,tl.col="white", tl.cex=0.0000001, cl.lim=c(-21,21),cl.pos="b")

cp_pl=corrplot(sigma(myPLN_plante),  is.corr = F, method="circle", addgrid.col = "lightgray", col=col,tl.col="white", tl.cex=0.0000001, cl.lim=c(-21,21),cl.pos="b")

cp=corrplot(sigma(myPLN_pl_famd), is.corr = F, method="circle", addgrid.col = "lightgray", col=col,tl.col="white", tl.cex=0.0000001, cl.lim=c(-21,21),cl.pos="b")

par(mfrow=c(1,1))

#Model selection (Table 1)
noms=cbind(Model=c("Model 1-0", "Model 1-1", "Model 1-2", "Model 1-3","Model 1-4","Model 1-5","Model 1-6"), Covariates=c("None", "Plant", "Plant", "Eco", "Eco", "Plant + Eco","Plant + Eco"), Residual_Matrix=c("Full","Full","Diagonal","Full","Diagonal","Full","Diagonal"))

tab_res=rbind(myPLN0_poids$criteria, 
							myPLN_plante$criteria,myPLN_plante_diag$criteria, 
							myPLN_famd$criteria,myPLN_famd_diag$criteria,
							myPLN_pl_famd$criteria,myPLN_pl_famd_diag$criteria)

tab_res=cbind(noms, tab_res)[,-7]

tab_res=tab_res[order(tab_res$BIC, decreasing = T),] 

tab_res$BIC=-2*tab_res$BIC

delta=sapply(1:length(tab_res$BIC),function(i){
	best=min(tab_res$BIC)
	sortie=tab_res$BIC[[i]]-best
})
names(delta)="Delta"

tab_res=cbind(tab_res,delta)

tab_res%>% knitr::kable()

#Best model : fitted vs observed (Fig S3)
tmp=data.frame(
	fitted   = as.vector(fitted(myPLN_pl_famd)+1),
	observed = as.vector(rn_pln$Abundance+1)
)
ggplot(tmp,aes(x = observed, y = fitted)) +
	geom_point(size = .5, alpha =.25 ) +
	scale_x_log10() +
	scale_y_log10() +
	theme_bw() + annotation_logticks()

#Best model : Species abundances'responses to ecological covariates (in fig 4)
coefs=data.frame(
	Species = rownames(coef(myPLN_pl_famd)),
	FAMD1=coef(myPLN_pl_famd)[,22],
	FAMD2=coef(myPLN_pl_famd)[,23],
	SD1=standard_error(myPLN_pl_famd)[,22],
	SD2=standard_error(myPLN_pl_famd)[,23], 
	spe=c("gen","gen","gen","gen","spe","spe","spe","spe") 
)

ggplot(coefs, aes(x = FAMD1, y=FAMD2))+
	geom_hline(yintercept=0)+
	geom_vline(xintercept=0)+
	geom_point(aes(color=spe))+
	geom_pointrange(aes(ymin=FAMD2-1.96*SD2, ymax=FAMD2+1.96*SD2, color=spe))+
	geom_pointrange(aes(xmin=FAMD1-1.96*SD1, xmax=FAMD1+1.96*SD1, color=spe))+
	scale_color_brewer(palette="Dark2")+
	theme_minimal()+
	geom_text_repel(aes(label = Species,  color = spe,fontface=3), size = 5, box.padding = 1.1)
	
#Network inference------
network_models <- PLNnetwork(Abundance ~ 1 + factor(plantesc) + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

model_BIC <- getBestModel(network_models, "BIC") 
net<-model_BIC$plot_network(type = "partial_cor",plot = F,edge.color = c("dodgerblue","firebrick"))
plot(net)

# Selection frequency in the bootstrap subsamples of the StARS procedure (from Pauvert et al. 2020 BioRxiv)
## need for a stratified sampling to be able to keep enough representatives of each host plant in each network
df <- rn_pln %>%
	mutate(i = 1:nrow(.)) %>% #create row number if you dont have one
	select(i, everything()) # put 'i' at the front of the dataset

subsmp<-replicate(100,df %>% # 100 subsamples
										group_by(plantesc) %>% #any number of variables you wish to partition by proportionally
										sample_frac(0.6) %>% pull(i),simplify = F) # without replacement 
# Reconstruction of subnetworks based on the 100 subsamples for each penalty
network_models$stability_selection(subsamples = subsmp,mc.cores = 1)
#Selection of the Stars network
model_StARS<-network_models$getBestModel("StARS")
net_stars<-model_StARS$plot_network(type = "partial_cor",plot = T,edge.color = c("dodgerblue","firebrick")) 
# extraction of the Selection frequency
probs<-extract_probs(network_models,penalty = model_BIC$penalty, format = "vector")
probs

#------INFERRING INTERACTIONS WITH SPECIES TRAITS (7 species community)------
#Loading data------
terrain7species_lab_eco<-read.csv2("data/FF7_abund_traits_eco.csv",dec=".")
nom_mouches<-read.csv2("data/fly_labels.csv")

#Famd on ecological covariates------
tmp=pivot_wider(data=terrain7species_lab_eco, id_cols = names(terrain7species_lab_eco)[c(1,6,7,25,28:34)], names_from="espece", values_from = "abondance")

dat_eco=tmp[,c(2:11)]
dat_eco$mois=as.factor(dat_eco$mois)
dat_eco$annee=as.factor(dat_eco$annee)

res.famd=FAMD(dat_eco, graph=F,ncp = 10)
ind<-res.famd$ind$coord
colnames(ind)=c("FAMD.1","FAMD.2","FAMD.3","FAMD.4","FAMD.5","FAMD.6","FAMD.7","FAMD.8","FAMD.9","FAMD.10")
ind=as.data.frame(ind)
ind$numech <-tmp$numech

donnees_long=left_join(terrain7species_lab_eco[,c(1,2,4,5,15,16,22,23)],ind, by="numech")

#Preparing data------
donnees=pivot_wider(data=donnees_long, id_cols = names(donnees_long)[c(1:4,9:18)], names_from="espece", values_from = "abondance")
rn_abondance<-as.matrix(donnees[,15:21])
colnames(rn_abondance)=nom_mouches$latin[match(names(donnees)[15:21], nom_mouches$correct)]
rownames(rn_abondance)=donnees$numech
rn_covariate<-donnees[,c(2:14)]
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#Model inference (for models without traits)------
myPLN0_base <- PLN(Abundance ~ 1, data = rn_pln)
myPLN0_poids <- PLN(Abundance ~ 1 + offset(log(poids)), data = rn_pln)

myPLN_plante <- PLN(Abundance ~ 1 + factor(plantesc) + offset(log(poids)), data = rn_pln)
myPLN_plante_diag <- PLN(Abundance ~ 1 + factor(plantesc) + offset(log(poids)), data = rn_pln, control = list(covariance = "diagonal"))

myPLN_famd <- PLN(Abundance ~ 1 +  FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)
myPLN_famd_diag <- PLN(Abundance ~ 1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10+ offset(log(poids)), data = rn_pln, control = list(covariance = "diagonal"))

myPLN_pl_famd <- PLN(Abundance ~ 0 + factor(plantesc) + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)
myPLN_pl_famd_diag <- PLN(Abundance ~ 0 + factor(plantesc) + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln, control = list(covariance = "diagonal"))

noms=cbind(Model=c("Model 2-0", "Model 2-1", "Model 2-2", "Model 2-3","Model 2-4","Model 2-5","Model 2-6"), Covariates=c("None", "Plant", "Plant", "Eco", "Eco", "Plant + Eco","Plant + Eco"), Residual_Matrix=c("Full","Full","Diagonal","Full","Diagonal","Full","Diagonal"))

tab_res=rbind(myPLN0_poids$criteria, 
							myPLN_plante$criteria,myPLN_plante_diag$criteria, 
							myPLN_famd$criteria,myPLN_famd_diag$criteria,
							myPLN_pl_famd$criteria,myPLN_pl_famd_diag$criteria)

tab_res=cbind(noms, tab_res)[,1:6]

tab_res%>% knitr::kable()

tab_res_base<-tab_res


#PLN models with traits - 7 species------
tcl=split(donnees_long, donnees_long$espece, drop = T)

#****Bacterocera_zonata------
#Abondances
donnees=tcl[["Bactrocera_zonata"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Bz")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_bz_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_bz_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_bz_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_bz_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_bz_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_bz_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_bz_logfec$criteria,
							myPLN_bz_logfit1$criteria,
							myPLN_bz_logfecfit1$criteria,
							myPLN_bz_logfec_eco$criteria,
							myPLN_bz_logfit1_eco$criteria,
							myPLN_bz_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_BZ<-tab_res

#****Zeugodacus_cucurbitae------
#Abondances
donnees=tcl[["Zeugodacus_cucurbitae"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Zc")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_zc_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_zc_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_zc_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_zc_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_zc_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_zc_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_zc_logfec$criteria,
							myPLN_zc_logfit1$criteria,
							myPLN_zc_logfecfit1$criteria,
							myPLN_zc_logfec_eco$criteria,
							myPLN_zc_logfit1_eco$criteria,
							myPLN_zc_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_ZC<-tab_res

#****Dacus_demmerezi------
#Abondances

donnees=tcl[["Dacus_demmerezi"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Dd")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_dd_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_dd_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_dd_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_dd_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_dd_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_dd_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_dd_logfec$criteria,
							myPLN_dd_logfit1$criteria,
							myPLN_dd_logfecfit1$criteria,
							myPLN_dd_logfec_eco$criteria,
							myPLN_dd_logfit1_eco$criteria,
							myPLN_dd_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_DD<-tab_res

#****Neoceratitis cyanescens------
#Abondances

donnees=tcl[["Neoceratitis_cyanescens"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Nc")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_nc_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_nc_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_nc_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_nc_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_nc_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_nc_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_nc_logfec$criteria,
							myPLN_nc_logfit1$criteria,
							myPLN_nc_logfecfit1$criteria,
							myPLN_nc_logfec_eco$criteria,
							myPLN_nc_logfit1_eco$criteria,
							myPLN_nc_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_NC<-tab_res

#****Ceratitis capitata------
#Abondances
donnees=tcl[["Ceratitis_capitata"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Cc")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_cc_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_cc_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_cc_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_cc_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_cc_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_cc_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_cc_logfec$criteria,
							myPLN_cc_logfit1$criteria,
							myPLN_cc_logfecfit1$criteria,
							myPLN_cc_logfec_eco$criteria,
							myPLN_cc_logfit1_eco$criteria,
							myPLN_cc_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_CC<-tab_res

#****Ceratitis catoirii------
#Abondances
donnees=tcl[["Ceratitis_catoirii"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Ct")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_ct_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_ct_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_ct_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_ct_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_ct_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_ct_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_ct_logfec$criteria,
							myPLN_ct_logfit1$criteria,
							myPLN_ct_logfecfit1$criteria,
							myPLN_ct_logfec_eco$criteria,
							myPLN_ct_logfit1_eco$criteria,
							myPLN_ct_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_CT<-tab_res

#****Ceratitis quilicii------
#Abondances
donnees=tcl[["Ceratitis_quilicii"]]
names(donnees)
rn_abondance<-as.matrix(donnees[,6])
colnames(rn_abondance)=c("Cq")
rownames(rn_abondance)=donnees$numech
str(rn_abondance)

#Covariables
names(donnees)
rn_covariate<-donnees[,c(2:4,5:18)]
names(rn_covariate)

#Mise en forme
rn_pln<-data.frame(rn_covariate)
rn_pln$Abundance=rn_abondance
str(rn_pln)

#pln traits
myPLN_cq_logfec <- PLN(Abundance ~ 1 + log_fec + offset(log(poids)), data = rn_pln)

myPLN_cq_logfit1 <- PLN(Abundance ~ 1 + fitness1 + offset(log(poids)), data = rn_pln)

myPLN_cq_logfecfit1 <- PLN(Abundance ~ 1 + log_fec + fitness1 + offset(log(poids)), data = rn_pln)

#pln traits + eco
myPLN_cq_logfec_eco <- PLN(Abundance ~ 1 + log_fec + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_cq_logfit1_eco <- PLN(Abundance ~ 1 + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

myPLN_cq_logfecfit1_eco <- PLN(Abundance ~ 1 + log_fec + fitness1 + FAMD.1 + FAMD.2 + FAMD.3 + FAMD.4 + FAMD.5 + FAMD.6 + FAMD.7 + FAMD.8 + FAMD.9 + FAMD.10 + offset(log(poids)), data = rn_pln)

noms=data.frame(Model=c("log(fec)", "log(fitness1)", "log(fec)+log(fitness1)",
												"log(fec) FAMD", "log(fitness1) FAMD", "log(fec)+log(fitness1) FAMD"))

tab_res=rbind(myPLN_cq_logfec$criteria,
							myPLN_cq_logfit1$criteria,
							myPLN_cq_logfecfit1$criteria,
							myPLN_cq_logfec_eco$criteria,
							myPLN_cq_logfit1_eco$criteria,
							myPLN_cq_logfecfit1_eco$criteria)

tab_res=cbind(noms$Model, tab_res)

tab_res[order(tab_res$BIC),] %>% knitr::kable()

tab_res_poids_CQ<-tab_res

#****Assembling models------
tab_res=rbind(tab_res_poids_BZ,tab_res_poids_CQ,tab_res_poids_CC,tab_res_poids_CT,tab_res_poids_ZC,tab_res_poids_DD,tab_res_poids_NC)

tab_mod= split(tab_res, tab_res$`noms$Model`, drop = T)

crit_summary=sapply(names(tab_mod),function(modele){
	crit=apply(tab_mod[[modele]][,2:3], 2,sum)
})

crit_summary=as.data.frame(t(crit_summary))

crit_summary$BIC=crit_summary$loglik-crit_summary$nb_param*log(dim(donnees)[1])/2

noms=cbind(
	Model=c("Model 2-7", "Model 2-8", "Model 2-9", "Model 2-10","Model 2-11","Model 2-12"),
	Covariates=c("Preference", "Preference + Eco", "Preference + Performance", "Preference + Performance + Eco", "Performance", "Performance + Eco"), Residual_Matrix=c("Diagonal","Diagonal","Diagonal","Diagonal","Diagonal","Diagonal"))

tab_res=cbind(noms, crit_summary)
rownames(tab_res)<-NULL
tab_res=rbind(tab_res,tab_res_base[,-7])
tab_res=tab_res[order(tab_res$BIC, decreasing = T),] 
tab_res$BIC=-2*tab_res$BIC

delta=sapply(1:length(tab_res$BIC),function(i){
	best=min(tab_res$BIC)
	sortie=tab_res$BIC[[i]]-best
})
names(delta)="Delta"

tab_res=cbind(tab_res,delta)
tab_res%>% knitr::kable()

#Species abundances' responses to host plants (fig 3)------
labels_plantes=read.csv2(file="data/plant_labels.csv", sep=";", header=F)
mouches_bon_ordre=c("Bactrocera zonata","Ceratitis quilicii","Ceratitis capitata","Ceratitis catoirii","Neoceratitis cyanescens","Dacus demmerezi","Zeugodacus cucurbitae")
plantes_bon_ordre=c("Mangifera indica","Annona reticulata","Terminalia catappa","Psidium cattleyanum","Psidium guajava","Syzygium jambos","Syzygium samarangense","Averrhoa carambola","Eriobotrya japonica","Prunus persica","Coffea arabica","Citrus reticulata","Citrullus lanatus","Cucumis melo","Cucumis sativus","Cucurbita maxima","Cucurbita pepo","Sechium edule","Capsicum annuum","Solanum lycopersicum","Solanum mauritianum")
niche_labo<-read.csv2(file="data/species_traits.csv", dec=".") #Laboratory-measured fitnesses

gfit=ggplot(data = niche_labo, aes(espece, plantesc, fill = fitness))+
	geom_tile(color = "lightgrey")+
	theme_minimal()+ 
	theme(axis.text.x = element_text(angle = 90, size=8,vjust = 1, hjust = 1, face="italic"),
				axis.text.y = element_text(face="italic"))+
	coord_fixed()+
	scale_y_discrete(name ="Plant species", labels=labels_plantes[seq(21,1),2],
									 limits=plantes_bon_ordre[seq(21,1)])+
	scale_x_discrete(name ="Fly species", 
									 limits=mouches_bon_ordre[c(6,7,5,4,3,2,1)])+
	scale_fill_gradient(low = "white", high = "black", space = "Lab", name="Log\nFitness")
gfit

#Inferred species abundances's responses to host plants
tmp=as.data.frame(coef(myPLN_pl_famd_diag)[,1:21])
tmp=cbind(tmp,rownames(tmp))
colnames(tmp)=c(sapply(colnames(tmp)[1:21], function(nom){str_sub(nom, start = 17)}),"espece")
beta_df=pivot_longer(tmp, cols=names(tmp)[1:21],names_to = "plantesc", values_to = "beta")

gbeta=ggplot(data = beta_df, aes(x=espece, y=plantesc, fill = beta))+
	geom_tile(color = "lightgrey")+
	theme_minimal()+ 
	theme(axis.text.x = element_text(angle = 90, size=8,vjust = 1, hjust = 1, face="italic"),
				axis.text.y = element_text(face="italic"))+
	coord_fixed()+
	scale_y_discrete(name ="Plant species", labels=labels_plantes[seq(21,1),2],
									 limits=plantes_bon_ordre[seq(21,1)])+
	scale_x_discrete(name ="Fly species", 
									 limits=mouches_bon_ordre[c(6,7,5,4,3,2,1)])+
	scale_fill_gradient(low = "white", high = "black", space = "Lab", name="Reponse")

gbeta

#Relationship between lab fitness and abundances'responses
beta_df=left_join(beta_df,niche_labo)
beta_df$spe2=sapply(beta_df$espece, function(esp){if(esp%in%c("Bactrocera zonata","Ceratitis quilicii","Ceratitis capitata","Ceratitis catoirii")) "Generalist" else "Specialist"})
beta_df$inferred_host=sapply(beta_df$beta, function(b){if(b <= -50) "low" else "Moderate"})

tmp=split(beta_df, beta_df$inferred_host, drop = T)

formula<-y~x
p <- ggplot(data = tmp[[2]], aes(x=fitness, y = beta, color=spe2))+
	geom_point()+
	stat_smooth(aes(fill = spe2, color = spe2), method = "lm", formula = formula)+
	stat_regline_equation(
		aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~")),
		formula = formula)+
	theme_minimal()

gcorr_coefs_fit=ggpar(p, palette = "Dark2",legend.title = "Host use strategy")
gcorr_coefs_fit