################################################################################
# 												R-script of the statistical analyses related to the article entitled : 												 #
# 														Vegetation structure of plantain-based agrosystems 															     #
#											determines numerical dominance in community of grounddwelling ants 										 #
# 																											by 																														 #
# 											Anicet Dassou, Philippe Tixier, Sylvain Dépigny and Dominique Carval 										 #
#																											(2017) 																												 #
################################################################################

#------------------------------------------------------------------------------#
{#Rawdata legend
#see also the Methods section and the table 2 of related manuscript
#------------------------------------------------------------------------------#
# sampling											number of the sampling date (1 = rainy season & 2 = dry season)
# plot														identification number of the plot (fields 1 to 20)
# subplot												identification number of the subplot (1 to 25 per plot)
# x																x -axis coordinate of the subplot
# y																y -axis coordinate of the subplot
# siteid													unique subplot identification number
# plantain												density of Musa AAB (plants/m²)
# oil_palm												density of Elais guineensis (plants/m²)
# papaya												density of Carica papaya (plants/m²)
# groundnut											density of Arachis hypogaea (plants/m²)
# taro	coffee										density of Colocasia esculenta (plants/m²)	
# yam														density of Dioscorea spp. (plants/m²)	
# maize													density of Zea mays (plants/m²)	
# macabo												density of Xanthosoma sagittifolium (plants/m²)	
# hot_pepper										density of Capsicum anuum (plants/m²)	
# cacao													density of Theobroma cacao (plants/m²)	
# garden_egg										density of Solanum macrocarpon (plants/m²)	
# crin-crin												density of Corchorus spp. (plants/m²)
# cassava												density of Manihot esculenta (plants/m²)
# banana												density of Musa AAA (plants/m²)
# pineapple											density of Ananas comosus (plants/m²)
# cola														density of Cola acuminata (plants/m²)
# safou													density of Dacryodes edulis (plants/m²)
# avocado												density of Persea americana (plants/m²)
# guava													density of Psidium guajava (plants/m²)
# vernonia												density of Vernonia spp. (plants/m²)
# amaranth											density of Amaranthus spp. (plants/m²)
# eru															density of Gnetum africanum (plants/m²)
# tomato												density of Solanum lycopersicum (plants/m²)
# gombo													density of Abelmoschus esculentus (plants/m²)
# niebe													density of Vigna unguiculata (plants/m²)
# sweet_potato									density of Ipomea batatas (plants/m²)
# triumphetta_pentadra				density of Triumphetta pentadra (plants/m²)
# mango													density of Mangifera indica (plants/m²)
# divp														diversity of plants in subplot
# Paratrechina_longicornis			Abundance of ants of the species Paratrechina longicornis at bait (number of individuals)
# Tetramorium_sp							Abundance of ants of the species Tetramorium sp. at bait (number of individuals)
# Axinidris_murielae						Abundance of ants of the specie Axinidris sp. at bait (number of individuals)
# Monomorium_spp							Abundance of ants of the genus Monomorium at bait (number of individuals)
# Monomorium_bicolor					Abundance of ants of the species Monomorium bicolor at bait (number of individuals)
# Monomorium_sp1							Abundance of ants of the species Monomorium sp. 1 at bait (number of individuals)
# Monomorium_sp2				    		Abundance of ants of the sepcies Monomorium sp. 2 at bait (number of individuals)
# Pheidole_spp									Abundance of ants of the genus Pheidole at bait (number of individuals)
# Pheidole_sp1							    	Abundance of ants of the genus Pheidole sp. 1 at bait (number of individuals)
# Pheidole_megacephala				Abundance of ants of the species Pheidole megacephala at bait (number of individuals)
# Camponotus_spp							Abundance of ants of the genus Camponotus at bait (number of individuals)
# Camponotus_acvapimensis		Abundance of ants of the species Camponotus acvapimensis at bait (number of individuals)
# Camponotus_brutus						Abundance of ants of the species Camponotus brutus at bait (number of individuals)
# Camponotus_sp1							Abundance of ants of the species Camponotus sp. 1 at bait (number of individuals)
# Camponotus_sp2							Abundance of ants of the species Camponotus sp. 2 at bait (number of individuals)
#Odontomachus_troglodytes		Abundance of ants of the species Odontomachus troglodytes at bait (number of individuals)
# Lpara													Log of abundance of ants of the species Paratrechina longicornis at bait
# Ltetra													Log of abundance of ants of the species Tetramorium sp. at bait 
# Laxin													Log of abundance of ants of the specie Axinidris sp. at bait 
# Lmono													Log of abundance of ants of the genus Monomorium at bait 
# Lphei													Log of abundance of ants of the genus Pheidole at bait
# Lcamp													Log of abundance of ants of the genus Camponotus at bait
# rankPara											rank values for the species Paratrechina longicornis according to its abundance 
# rankTetra											rank values for the species Tetramorium sp. according to its abundance 
# rankAxin											rank values for the species Axinidris sp. according to its abundance 
# rankMono											rank values for the genus Monomorium according to its abundance 
# rankPhei												rank values for the genus Pheidole according to its abundance 
# rankCamp											rank values for the genus Camponotus longicornis according to its abundance 
# DomSpecies										Identification number of the Dominant Species (see below line 105)
# High														density of the high stratum (plants/m²) 
# Intermediate									density of intermediate stratum (plants/m²)
# Low														density of low stratum (plants/m²)
# Musa													density of Musa (plants/m²)
}

#------------------------------------------------------------------------------#
{# Dominance analysis
#------------------------------------------------------------------------------#

	#------------------------------------------------------------------------------#
	{# Occurence of ant species
	#------------------------------------------------------------------------------#
		Bait<-rawdata
		BaitBin<-Bait

		nbrow<-length(Bait[,1])

		for(i in 1:nbrow)
		{
			for(j in 38:51)
			{
				if (Bait[i,j] > 0)
				{
					BaitBin[i,j] = 1
				}
			}
		}

		#Overall frequency
		dataAll=BaitBin[,c(38:51)]
		recorded<- colMeans(dataAll, na.rm = T)

		{#Frequency by season of sampling
		#Rainy season (sampling = 1)
		subdata1NA<-subset(BaitBin, BaitBin$sampling == 1)
		#Dry season (sampling = 2)
		subdata2NA<-subset(BaitBin, BaitBin$sampling == 2)

		#Selection of the column of species
		dataNA1=subdata1NA[,c(38:51)]
		dataNA2=subdata2NA[,c(38:51)]

		tempsNA1<- colMeans(dataNA1, na.rm = T)
		tempsNA2<- colMeans(dataNA2, na.rm = T)

		BinDate<-rbind(tempsNA1, tempsNA2)
		row.names(BinDate)<-c("Rainy season","Dry season")
		}
		#To display the results
		BinDate
	}

	#------------------------------------------------------------------------------#
	{# % of controled baits 
	# Baits were considered controlled by a species (i) if the number 
	# of individuals is >20 and no other ant was present (monopolization) 
	# or (ii) if one species was at least twice as numerous as the second 
	# numerous taxon when several species were present and the total 
	# number of individuals was >20)
	#------------------------------------------------------------------------------#

		Bait<-rawdata
		# We remove the abundance at the genus level except for Pheidole
		# (i.e. Pheidole_spp, Camponotus_spp & Monomorium_spp)
		Bait<-Bait[,-c(41,46)]
		BaitCont1<-Bait

		nbrow<-length(Bait[,1])

		for(i in 1:nbrow)
		{
			for(j in 38:49)
			{
				# if the species "j" is present
				if (Bait[i,j] > 0) 
				{
					#if there are more thant 20 individuals of the species "j" & if it is the only species present
					if (Bait[i,j] >= 20 & sum(Bait[i,c(38:49)]) - Bait[i,j] == 0) 
					{
						BaitCont1[i,j] = 1 # = monopolization
					}
					else 
					{
						# if there are more thant 20 individuals (all species considered)
						if  (sum(Bait[i,c(38:49)] ) > 20) 
						{
							second = 0
							for (k in c(38:49))
							{ 
							#We consider abundance of other species than "j" ...
							if (k != j) 
								# ... to see which species is the most abundant
								second = max(second, Bait[i,k]) 
							}
							# if the abundance of "j" is at least twice as numerous as the second numerous species
							if (Bait[i,j] > 2 * second)
							{
								BaitCont1[i,j] = 1 # = control
							}
							else
								BaitCont1[i,j] = 0 # = no control
						}
						else
							BaitCont1[i,j] = 0 
					} 
				}
				else BaitCont1[i,j] = NA
			}
		}

		#Overall control
		mean(BaitCont1$Pheidole_spp, na.rm = T)
		mean(BaitCont1$Paratrechina_longicornis, na.rm = T)
		mean(BaitCont1$Tetramorium_sp, na.rm = T)
		mean(BaitCont1$Axinidris_murielae, na.rm = T)
		mean(BaitCont1$Monomorium_bicolor, na.rm = T)
		mean(BaitCont1$Monomorium_sp1, na.rm = T)
		mean(BaitCont1$Monomorium_sp2, na.rm = T)
		mean(BaitCont1$Camponotus_acvapimensis, na.rm = T)
		mean(BaitCont1$Camponotus_brutus, na.rm = T)
		mean(BaitCont1$Camponotus_sp1, na.rm = T)
		mean(BaitCont1$Camponotus_sp2, na.rm = T)
		mean(BaitCont1$Odontomachus_troglodytes, na.rm = T)


		{#Frequency by season of sampling
		#Rainy season (sampling = 1)
		subdata1NA<-subset(BaitCont1, BaitCont1$sampling == 1)
		#Dry season (sampling = 2)
		subdata2NA<-subset(BaitCont1, BaitCont1$sampling == 2)

		#Selection of the column of species
		dataNA1=subdata1NA[,c(38:49)]
		dataNA2=subdata2NA[,c(38:49)]

		tempsNA1<- colMeans(dataNA1, na.rm = T)
		tempsNA2<- colMeans(dataNA2, na.rm = T)

		BaitContDate<-rbind(tempsNA1, tempsNA2)
		row.names(BaitContDate)<-c("Rainy season","Dry season")
		}
		#To display the results
		BaitContDate
	}

	#------------------------------------------------------------------------------#
	{# Mean abundance scores
	#------------------------------------------------------------------------------#
		Bait<-rawdata
		BaitScore<-Bait

		nbrow<-length(Bait[,1])

		#Creation of scores according to the 6 point abundnace scale
		for(i in 1:nbrow)
		{
			  for(j in 38:51)
			  {
					if (Bait[i,j] == 1)
					{
						BaitScore[i,j] = 1
					}
					if (Bait[i,j] > 1 & Bait[i,j] <= 5)
					{
						BaitScore[i,j] = 2
					}
					if (Bait[i,j] > 6  & Bait[i,j] <= 10)
					{
						BaitScore[i,j] = 3
					}
					if (Bait[i,j] > 11 & Bait[i,j] <= 20)
					{
						BaitScore[i,j] = 4
					}
					if (Bait[i,j] > 21 & Bait[i,j] <= 49)
					{
						BaitScore[i,j] = 5
					}
					if (Bait[i,j] >= 50)
					{
						BaitScore[i,j] = 6
					}
					if (Bait[i,j] == 0)
					{
						BaitScore[i,j] = NA
					}
				}
		}

			#Overall mean abundance scores
			mean(BaitScore$Pheidole_spp, na.rm = T)
			mean(BaitScore$Paratrechina_longicornis, na.rm = T)
			mean(BaitScore$Tetramorium_sp, na.rm = T)
			mean(BaitScore$Axinidris_murielae, na.rm = T)
			mean(BaitScore$Monomorium_spp, na.rm = T)
			mean(BaitScore$Monomorium_bicolor, na.rm = T)
			mean(BaitScore$Monomorium_sp1, na.rm = T)
			mean(BaitScore$Monomorium_sp2, na.rm = T)
			mean(BaitScore$Camponotus_spp, na.rm = T)
			mean(BaitScore$Camponotus_acvapimensis, na.rm = T)
			mean(BaitScore$Camponotus_brutus, na.rm = T)
			mean(BaitScore$Camponotus_sp1, na.rm = T)
			mean(BaitScore$Camponotus_sp2, na.rm = T)
			mean(BaitScore$Odontomachus_troglodytes, na.rm = T)

		{#Frequency by season of sampling
		#Rainy season (sampling = 1)
		subdata1NA<-subset(BaitScore, BaitScore$sampling == 1)
		#Dry season (sampling = 2)
		subdata2NA<-subset(BaitScore, BaitScore$sampling == 2)

		#Selection of the column of species
		dataNA1=subdata1NA[,c(38:51)]
		dataNA2=subdata2NA[,c(38:51)]

		tempsNA1<- colMeans(dataNA1, na.rm = T)
		tempsNA2<- colMeans(dataNA2, na.rm = T)

		BaitScoreDate<-rbind(tempsNA1, tempsNA2)
		row.names(BaitScoreDate)<-c("Rainy season","Dry season")
		}
		#To display the results
		BaitScoreDate
	}

#------------------------------------------------------------------------------#
# End of dominance analysis
}
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
{# Influence of the season on the occurence of genus at baits
#------------------------------------------------------------------------------#
	library(lme4)
	#Rainy season (sampling = 1)
	#Dry season (sampling = 2)
	
	#Pheidole
		Ph0<-glm(Pheidole_spp~sampling, family="binomial", data=BaitBin)
		drop1(Ph0, test = "Chisq")
		
		Ph1<-glm(Pheidole_megacephala~sampling, family="binomial", data=BaitBin)
		drop1(Ph1, test = "Chisq")

		Ph2<-glm(Pheidole_sp1~sampling, family="binomial", data=BaitBin)
		drop1(Ph2, test = "Chisq")
	
	#Paratrechina
		Pl0<-glm(Paratrechina_longicornis~sampling, family="binomial", data=BaitBin)
		drop1(Pl0, test = "Chisq")

	#Axinidris
		Axin0<-glm(Axinidris_murielae~sampling, family="binomial", data=BaitBin)
		drop1(Axin0, test = "Chisq")

	#Tetramorium
		Tetra0<-glm(Tetramorium_sp~sampling, family="binomial", data=BaitBin)
		drop1(Tetra0, test = "Chisq")

	#Monomorium
		Mon0<-glm(Monomorium_spp~sampling, family="binomial", data=BaitBin)
		drop1(Mon0, test = "Chisq")
		
		Mon1<-glm(Monomorium_bicolor~sampling, family="binomial", data=BaitBin)
		drop1(Mon1, test = "Chisq")
		
		Mon2<-glm(Monomorium_sp1~sampling, family="binomial", data=BaitBin)
		drop1(Mon2, test = "Chisq")
		
		Mon3<-glm(Monomorium_sp2~sampling, family="binomial", data=BaitBin)
		drop1(Mon3, test = "Chisq")

	#Camponotus
		Camp0<-glm(Camponotus_spp~sampling, family="binomial", data=BaitBin)
		drop1(Camp0, test = "Chisq")
		
		Camp1<-glm(Camponotus_acvapimensis~sampling, family="binomial", data=BaitBin)
		drop1(Camp1, test = "Chisq")
		
		Camp2<-glm(Camponotus_brutus~sampling, family="binomial", data=BaitBin)
		drop1(Camp2, test = "Chisq")
		
		Camp3<-glm(Camponotus_sp1~sampling, family="binomial", data=BaitBin)
		drop1(Camp3, test = "Chisq")
		
		Camp4<-glm(Camponotus_sp2~sampling, family="binomial", data=BaitBin)
		drop1(Camp4, test = "Chisq")

		#Odontomachus troglodytes
		Odont0<-glm(Odontomachus_troglodytes~sampling, family="binomial", data=BaitBin)
		drop1(Odont0, test = "Chisq")

#------------------------------------------------------------------------------#
# End of analysis on Influence of the season on the occurence 
# of genus at baits
}
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
{# Figure 1 - Frequency of numerical dominance                                    
#------------------------------------------------------------------------------#

	#------------------------------------------------------------------------------#
	# Load the rawdata file
	#------------------------------------------------------------------------------#
	# !!! Set your own work directory that contains the rawdata file
	setwd("D:/Mes donnees/Dropbox/Articles/Revised/Dassou-et-al")

	rawdata<-read.table("rawdata.txt",h=T,sep="")
	rawdata$plot=as.factor(rawdata$plot)
	rawdata$sampling=as.factor(rawdata$sampling)
	rawdata$subplot<-as.factor(as.character(rawdata$subplot))
	rawdata$siteid<-as.factor(rawdata$siteid)

	#------------------------------------------------------------------------------#
	# All codominace cases (i.e. DomSpecies > 6) are turned 
	# into one category
	for (i in 1:length(rawdata[,1]))
	{
		if (rawdata$DomSpecies[i] >= 7)
		rawdata$DomSpecies[i]<-7
	}

	#We remove the codominance cases
	subrawdata<-subset(rawdata, rawdata$DomSpecies < 7)

	#We revalue the factor
	subrawdata$DomSpecies<-as.factor(subrawdata$DomSpecies)
	library(plyr)
	subrawdata$DomSpecies<-revalue(subrawdata$DomSpecies, c("5"="5", "1"="2", "2"="1", "3"="0", "4"="3", "6"="4", "0"="6"))
	subrawdata$DomSpecies
	subrawdata$DomSpecies<-as.numeric(as.character(subrawdata$DomSpecies))

	#Graphical parameters
	breakmap<-c(-0.5,0.5,1.5,2.5,3.5,4.5,5.5,6.5)
	coloration=c("gray0","gray17","gray34","gray51","gray68","gray85","gray100")
	cnames=c("gray0","gray17","gray34","gray51","gray68","gray85","gray100")

	#Legend names
	vector.names<-c(expression(italic(Axinidris~murielae)), expression(italic(Tetramorium)~sp.), expression(italic(P.~longicornis)), expression(italic(Monomorium)~spp), expression(italic(Camponotus)~spp.), expression(italic(Pheidole)~spp.),expression(Empty))

	x11()
	#Uncomment the code line below if you want to save image in .tiff
	tiff("Figure1Unscaled.tiff", width=6.5, height=5, units="in", res=900)
	#To draw 2 graphics, one for each season of sampling:
	opar=par(ps=8,mfrow=c(1,2), mar = c(0, 0, 0, 1))
	subrawdata1<-subset(subrawdata, subrawdata$sampling=="1")  
	hist(subrawdata1$DomSpecies, breaks=breakmap, col=coloration,freq=F,main=NULL,  axes=F, xaxt = 'n', labels=T,xlab="", ylab = NULL)
	par(ps=10)
	mtext("Rainy season", side=1, line=-1)
	legend("topleft",vector.names, col=coloration, ,fill=coloration, bty="n", border="black", box.lty=0)
	par(ps=8)
	subrawdata2<-subset(subrawdata, subrawdata$sampling=="2")
	hist(subrawdata2$DomSpecies, breaks=breakmap, col=coloration,freq=F,main=NULL,  axes=F, xaxt = 'n', labels=T,xlab="", ylab = NULL)
	par(ps=10)
	mtext("Dry season", side=1, line=-1)
	#Uncomment the code line below if you want to save image in .tiff
	dev.off()

#------------------------------------------------------------------------------#
# End of Figure 1.
}
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
{# Multinomial modelling & Figure 2.                                
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
# Load the rawdata file
#------------------------------------------------------------------------------#
# !!! Set your own work directory that contains the rawdata file
setwd("D:/Mes donnees/Dropbox/Articles/Revised/Dassou-et-al")

library(VGAM)
rawdata<-read.table("rawdata.txt",h=T,sep="")

rawdata$plot=as.factor(rawdata$plot)
rawdata$sampling=as.factor(rawdata$sampling)
rawdata$subplot<-as.factor(as.character(rawdata$subplot))
rawdata$siteid<-as.factor(rawdata$siteid)

#------------------------------------------------------------------------------#
# Remove the codominance or empty data
#------------------------------------------------------------------------------#
Dominance<-subset(rawdata, rawdata$DomSpecies < 7 & rawdata$DomSpecies != 0)
Dominance$DomSpecies<-as.factor(Dominance$DomSpecies)

length(Dominance[,1])

#--> There must remain 844 lines on the initial 1000

################################################################################
# Best Multinomial Model Selection
################################################################################

#------------------------------------------------------------------------------#
# Pheidole is taken as a category of reference
#------------------------------------------------------------------------------#
levels(Dominance$DomSpecies) <- list(Pheidole = "5", Paratrechina = "1", Tetramorium = "2", Axinidris = "3", Monomorium = "4", Camponotus = "6")
myrefLevel <- with(Dominance, DomSpecies[4])

#------------------------------------------------------------------------------#
# Complete model : Intercepts + Low stratum + Intermediate stratum + High stratum + Musa stratum + plant diversity
#------------------------------------------------------------------------------#
library(VGAM)

Selection1<-data.frame(Covariate = c("Plant diversity", "Musa", "High stratum", "Intermediate stratum", "Low stratum", "Intercept"), Deltadf = 0, Chi2 = 0, Pvalue = 0, Significance = 0)
fit1 <- vglm(DomSpecies ~ Low + Intermediate + High + Musa + divp, family = multinomial(refLevel = myrefLevel), Dominance) #complete
fit2 <- vglm(DomSpecies ~ Low + Intermediate + High + Musa, family = multinomial(refLevel = myrefLevel), Dominance) # test divp
fit3 <- vglm(DomSpecies ~ Low + Intermediate + High + divp, family = multinomial(refLevel = myrefLevel), Dominance) #test Musa
fit4 <- vglm(DomSpecies ~ Low + Intermediate + Musa + divp, family = multinomial(refLevel = myrefLevel), Dominance) #test High
fit5 <- vglm(DomSpecies ~ Low + High + Musa + divp, family = multinomial(refLevel = myrefLevel), Dominance) #test Medium
fit6 <- vglm(DomSpecies ~ Intermediate + High + Musa + divp, family = multinomial(refLevel = myrefLevel), Dominance) # test Low
fit7 <- vglm(DomSpecies ~ -1 + Low + Intermediate + High + Musa + divp, family = multinomial(refLevel = myrefLevel), Dominance) #test intercept

LRT1<-lrtest(fit1,fit2)
LRT2<-lrtest(fit1,fit3)
LRT3<-lrtest(fit1,fit4)
LRT4<-lrtest(fit1,fit5)
LRT5<-lrtest(fit1,fit6)
LRT6<-lrtest(fit1,fit7)

Selection1$Deltadf[1]<-LRT1@Body$Df[2]
Selection1$Chi2[1]<-LRT1@Body$Chisq[2]
Selection1$Pvalue[1]<-LRT1@Body$"Pr(>Chisq)"[2]
Selection1$Deltadf[2]<-LRT2@Body$Df[2]
Selection1$Chi2[2]<-LRT2@Body$Chisq[2]
Selection1$Pvalue[2]<-LRT2@Body$"Pr(>Chisq)"[2]
Selection1$Deltadf[3]<-LRT3@Body$Df[2]
Selection1$Chi2[3]<-LRT3@Body$Chisq[2]
Selection1$Pvalue[3]<-LRT3@Body$"Pr(>Chisq)"[2]
Selection1$Deltadf[4]<-LRT4@Body$Df[2]
Selection1$Chi2[4]<-LRT4@Body$Chisq[2]
Selection1$Pvalue[4]<-LRT4@Body$"Pr(>Chisq)"[2]
Selection1$Deltadf[5]<-LRT5@Body$Df[2]
Selection1$Chi2[5]<-LRT5@Body$Chisq[2]
Selection1$Pvalue[5]<-LRT5@Body$"Pr(>Chisq)"[2]
Selection1$Deltadf[6]<-LRT6@Body$Df[2]
Selection1$Chi2[6]<-LRT6@Body$Chisq[2]
Selection1$Pvalue[6]<-LRT6@Body$"Pr(>Chisq)"[2]


for (i in 1:length(Selection1[,1]))
{
  if (Selection1$Pvalue[i] < 0.05) Selection1$Significance[i] = "**"
  else  Selection1$Significance[i] = "NS"
}

Selection1
#write.table(Selection1, "Selection1.txt", sep=";")

#------------------------------------------------------------------------------#
# Remove Low stratum and plant diversity from the model
# New sub-model : Intercepts + Intermediate stratum + High stratum + Musa stratum
#------------------------------------------------------------------------------#

Selection2<-data.frame(Covariate = c("Musa", "High stratum", "Intermediate" , "Intercept"), Deltadf = 0, Chi2 = 0, Pvalue = 0, Significance = 0)

fit1 <- vglm(DomSpecies ~ Intermediate + High + Musa, family = multinomial(refLevel = myrefLevel), Dominance) #complete
fit2 <- vglm(DomSpecies ~ Intermediate + High, family = multinomial(refLevel = myrefLevel), Dominance) # test Musa
fit3 <- vglm(DomSpecies ~ Intermediate + Musa, family = multinomial(refLevel = myrefLevel), Dominance) #test High
fit4 <- vglm(DomSpecies ~ High + Musa, family = multinomial(refLevel = myrefLevel), Dominance) #test Intermediate
fit5 <- vglm(DomSpecies ~- 1 + Intermediate + High + Musa, family = multinomial(refLevel = myrefLevel), Dominance) #test Intercept


LRT1<-lrtest(fit1,fit2)
LRT2<-lrtest(fit1,fit3)
LRT3<-lrtest(fit1,fit4)
LRT4<-lrtest(fit1,fit5)

Selection2$Deltadf[1]<-LRT1@Body$Df[2]
Selection2$Chi2[1]<-LRT1@Body$Chisq[2]
Selection2$Pvalue[1]<-LRT1@Body$"Pr(>Chisq)"[2]
Selection2$Deltadf[2]<-LRT2@Body$Df[2]
Selection2$Chi2[2]<-LRT2@Body$Chisq[2]
Selection2$Pvalue[2]<-LRT2@Body$"Pr(>Chisq)"[2]
Selection2$Deltadf[3]<-LRT3@Body$Df[2]
Selection2$Chi2[3]<-LRT3@Body$Chisq[2]
Selection2$Pvalue[3]<-LRT3@Body$"Pr(>Chisq)"[2]
Selection2$Deltadf[4]<-LRT3@Body$Df[2]
Selection2$Chi2[4]<-LRT3@Body$Chisq[2]
Selection2$Pvalue[4]<-LRT3@Body$"Pr(>Chisq)"[2]

for (i in 1:length(Selection2[,1]))
{
  if (Selection2$Pvalue[i] < 0.05) Selection2$Significance[i] = "**"
  else  Selection2$Significance[i] = "NS"
}

Selection2
#write.table(Selection2, "Selection2.txt", sep=";")

#------------------------------------------------------------------------------#
# Remove Musa stratum from the model
# New sub-model : Intercepts + Intermediate stratum + High stratum 
#------------------------------------------------------------------------------#

Selection3<-data.frame(Covariate = c("High stratum", "Intermediate stratum" , "Intercept"), Deltadf = 0, Chi2 = 0, Pvalue = 0, Significance = 0)

fit1 <- vglm(DomSpecies ~ Intermediate + High, family = multinomial(refLevel = myrefLevel), Dominance)
fit2 <- vglm(DomSpecies ~ Intermediate, family = multinomial(refLevel = myrefLevel), Dominance)
fit3 <- vglm(DomSpecies ~ High, family = multinomial(refLevel = myrefLevel), Dominance)
fit4 <- vglm(DomSpecies ~ -1 + Intermediate + High, family = multinomial(refLevel = myrefLevel), Dominance)

LRT1<-lrtest(fit1,fit2)
LRT2<-lrtest(fit1,fit3)
LRT3<-lrtest(fit1,fit4)

Selection3$Deltadf[1]<-LRT1@Body$Df[2]
Selection3$Chi2[1]<-LRT1@Body$Chisq[2]
Selection3$Pvalue[1]<-LRT1@Body$"Pr(>Chisq)"[2]
Selection3$Deltadf[2]<-LRT2@Body$Df[2]
Selection3$Chi2[2]<-LRT2@Body$Chisq[2]
Selection3$Pvalue[2]<-LRT2@Body$"Pr(>Chisq)"[2]
Selection3$Deltadf[3]<-LRT3@Body$Df[2]
Selection3$Chi2[3]<-LRT3@Body$Chisq[2]
Selection3$Pvalue[3]<-LRT3@Body$"Pr(>Chisq)"[2]

for (i in 1:length(Selection3[,1]))
{
  if (Selection3$Pvalue[i] < 0.05) Selection3$Significance[i] = "**"
  else  Selection3$Significance[i] = "NS"
}

Selection3
#write.table(Selection3, "Selection3.txt", sep=";")

#------------------------------------------------------------------------------#
# PREDICTIONS
#------------------------------------------------------------------------------#
#levels(Dominance$DomSpecies) <- list(Pheidole = "5", Paratrechina = "1", Tetramorium = "2", Axinidris = "3", Monomorium = "4", Camponotus = "6")
#myrefLevel <- with(Dominance, DomSpecies[4])

#High stratum density effect
GoodFit <- vglm(DomSpecies ~ High + Intermediate, family = multinomial(refLevel = myrefLevel), Dominance)
MyGrid<-expand.grid(DomSpecies = "Pheidole", High = seq(0,3.5,by=0.001), Intermediate = 0, Musa = 0, stratum = 1)
predictions1<- predictvglm(GoodFit, MyGrid, type="response")
predictions1<-as.data.frame(predictions1)
test1<-as.data.frame(cbind(MyGrid$High,predictions1[,1]))
test1$species<-"Pheidole"
test1$group<-"Trees"

test2<-as.data.frame(cbind(MyGrid$High,predictions1[,2]))
test2$species<-"Paratrechina"
test2$group<-"Trees"

test3<-as.data.frame(cbind(MyGrid$High,predictions1[,3]))
test3$species<-"Tetramorium"
test3$group<-"Trees"

test4<-as.data.frame(cbind(MyGrid$High,predictions1[,4]))
test4$species<-"Axinidris"
test4$group<-"Trees"

test5<-as.data.frame(cbind(MyGrid$High,predictions1[,5]))
test5$species<-"Monomorium"
test5$group<-"Trees"

test6<-as.data.frame(cbind(MyGrid$High,predictions1[,6]))
test6$species<-"Camponotus"
test6$group<-"Trees"

MyPred<-rbind(test1, test2, test3, test4, test5, test6)
MyPred$species<-as.factor(MyPred$species)
MyPred$group<-as.factor(MyPred$group)
names(MyPred)[1]<-"PlantDensity"
names(MyPred)[2]<-"Values"

#Intermediate stratum density effect
MyGrid<-expand.grid(DomSpecies = "Pheidole", High = 0, Intermediate = seq(0,3.5,by=0.01), Musa = 0, stratum = 2)
predictions2<- predictvglm(GoodFit, MyGrid, type="response")
predictions2<-as.data.frame(predictions2)
test1<-as.data.frame(cbind(MyGrid$Intermediate,predictions2[,1]))
test1$species<-"Pheidole"
test1$group<-"Shrubs"

test2<-as.data.frame(cbind(MyGrid$Intermediate,predictions2[,2]))
test2$species<-"Paratrechina"
test2$group<-"Shrubs"

test3<-as.data.frame(cbind(MyGrid$Intermediate,predictions2[,3]))
test3$species<-"Tetramorium"
test3$group<-"Shrubs"

test4<-as.data.frame(cbind(MyGrid$Intermediate,predictions2[,4]))
test4$species<-"Axinidris"
test4$group<-"Shrubs"

test5<-as.data.frame(cbind(MyGrid$Intermediate,predictions2[,5]))
test5$species<-"Monomorium"
test5$group<-"Shrubs"

test6<-as.data.frame(cbind(MyGrid$Intermediate,predictions2[,6]))
test6$species<-"Camponotus"
test6$group<-"Shrubs"

MyPred2<-rbind(test1, test2, test3, test4, test5, test6)
MyPred2$species<-as.factor(MyPred2$species)
MyPred2$group<-as.factor(MyPred2$group)
names(MyPred2)[1]<-"PlantDensity"
names(MyPred2)[2]<-"Values"

MyPred<-rbind(MyPred, MyPred2)
write.table(MyPred, "MyPred.txt", sep=";")

#------------------------------------------------------------------------------#
# PLOT OF PREDICTIONS
#------------------------------------------------------------------------------#
library(lattice)
vector.names<-c(expression(italic(A.~murielae)),expression(italic(Camponotus)~spp.), expression(italic(Monomorium)~spp.), expression(italic(P.~longicornis)),expression(italic(Pheidole)~spp.) ,expression(italic(Tetramorium)~sp.))

stripcol<-rgb(0.95,0.95,0.95)
opar=par(ps=10, mar = c(5, 4, 0, 0))
#Uncomment the code line below if you want to save image in .tiff
#tiff("Figure2Unscaled.tiff", width=6.5, height=5, units="in", res=900)
xyplot(Values~PlantDensity|species, groups=group, type = "l", lwd=2, xlab=list("Density of plants (individuals/m²)", cex=0.9), ylab=list("Probability of numerical dominance", cex = 0.9), lty = c(1, 1), col = c("black","darkgray"),  strip = strip.custom(strip.names = F, factor.levels = vector.names), par.strip.text = list(cex = 0.75),  par.settings = list(strip.background=list(col=stripcol)), scales = list(x = list(alternating=c(1,1), tck=c(-0.5,0)), y = list(alternating=c(1,1)), tck= c(-0.5,0)), main = NULL,data=MyPred)
#Uncomment the code line below if you want to save image in .tiff
#dev.off()

#------------------------------------------------------------------------------#
}# End of Multinomial modelling & Figure 2.                                
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
{# Figure S1.
#------------------------------------------------------------------------------#

	#------------------------------------------------------------------------------#
	# Load the rawdata file
	#------------------------------------------------------------------------------#
	# !!! Set your own work directory that contains the rawdata file
	setwd("D:/Mes donnees/Dropbox/Articles/Revised/Dassou-et-al")
	rawdata<-read.table("rawdata.txt",h=T,sep="")
	
	#------------------------------------------------------------------------------#
	# Graphical parameters
	#------------------------------------------------------------------------------#

	ParaCol="coral"
	PheiCol ="azure2"
	MonoCol = "darkseagreen1"
	CampCol="mistyrose1" 
	TetraCol="plum3"
	AxinCol= "indianred"

	#------------------------------------------------------------------------------#
	# Barplots
	#------------------------------------------------------------------------------#

	library(gplots)
	rawdata$sampling<-as.numeric(as.character(rawdata$sampling))

	x11()
	#Uncomment the code line below if you want to save image in .tiff
	#tiff("FigS1Unscaled.tiff", width=6.5, height=5, units="in", res=600)
	par(mfrow=c(2,3), mar = c(0, 3, 1, 1))           
	#Axinidris
	means<-(t(tapply(rawdata$Laxin,list(rawdata$sampling),mean)))
	stdev<-(t(tapply(rawdata$Laxin,list(rawdata$sampling),sd)))
	nsampling<-t(table(rawdata$sampling))
	ciw   <- qt(0.95, nsampling) * stdev / sqrt(nsampling)
	xscale<-(t(tapply(rawdata$sampling,list(rawdata$sampling),mean)))
	barplot2(means, plot.ci=TRUE,col=AxinCol, ci.l= means - ciw, ci.u= means + ciw,ci.color="black", axisnames = F, axes=F,xlab=NULL, ylab="", ci.lty = "solid", ci.lwd = 1)
	par(ps=9)
	axis(side = 2, tcl = 0.25, las = 1, labels = F)
	axis(side = 2, lwd = 0, las = 1, labels = T, line = -0.5)

	ynames<-c(expression(Abundance~of~italic(Axinidris~murielae)))
	par(ps=8)
	mtext(ynames, side=2, line=1.6)

	#Camponotus
	means<-(t(tapply(rawdata$Lcamp,list(rawdata$sampling),mean)))
	stdev<-(t(tapply(rawdata$Lcamp,list(rawdata$sampling),sd)))
	nsampling<-t(table(rawdata$sampling))
	ciw   <- qt(0.95, nsampling) * stdev / sqrt(nsampling)
	xscale<-(t(tapply(rawdata$sampling,list(rawdata$sampling),mean)))
	barplot2(means, plot.ci=TRUE,col=CampCol, ci.l= means - ciw, ci.u= means + ciw,ci.color="black", axisnames = F, axes=F, xlab="", ylab="", ci.lty = "solid", ci.lwd = 1)

	par(ps=9)
	axis(side = 2, tcl = 0.25, las = 1, labels = F)
	axis(side = 2, lwd = 0, las = 1, labels = T, line = -0.5)

	ynames<-c(expression(Abundance~of~italic(Camponotus)~spp.))
	par(ps=8)
	mtext(ynames, side=2, line=1.6)

	#Tetramorium
	means<-(t(tapply(rawdata$Ltetra,list(rawdata$sampling),mean)))
	stdev<-(t(tapply(rawdata$Ltetra,list(rawdata$sampling),sd)))
	nsampling<-t(table(rawdata$sampling))
	ciw   <- qt(0.95, nsampling) * stdev / sqrt(nsampling)
	xscale<-(t(tapply(rawdata$sampling,list(rawdata$sampling),mean)))
	barplot2(means, plot.ci=TRUE,col=TetraCol, ci.l= means - ciw, ci.u= means + ciw,ci.color="black",axisnames = F, axes=F, xlab="", ylab="", ci.lty = "solid", ci.lwd = 1)

	par(ps=9)
	axis(side = 2, tcl = 0.25, las = 1, labels = F)
	axis(side = 2, lwd = 0, las = 1, labels = T, line = -0.5)

	ynames<-c(expression(Abundance~of~italic(Tetramorium)~sp.))
	par(ps=8)
	mtext(ynames, side=2, line=1.6)

	par(mar = c(2, 3, 1, 1))

	#Pheidole
	means<-(t(tapply(rawdata$Lphei,list(rawdata$sampling),mean)))
	stdev<-(t(tapply(rawdata$Lphei,list(rawdata$sampling),sd)))
	nsampling<-t(table(rawdata$sampling))
	ciw   <- qt(0.95, nsampling) * stdev / sqrt(nsampling)
	xscale<-(t(tapply(rawdata$sampling,list(rawdata$sampling),mean)))
	barplot2(means, plot.ci=TRUE,col=PheiCol, ci.l= means - ciw, ci.u= means + ciw,ci.color="black",axisnames = F, axes=F, xlab="", ylab="", ci.lty = "solid", ci.lwd = 1)
	par(ps=8)
	mtext("Dry Season            Rainy season", side=1, line=0.25)

	par(ps=9)
	axis(side = 2, tcl = 0.25, las = 1, labels = F)
	axis(side = 2, lwd = 0, las = 1, labels = T, line = -0.5)

	ynames<-c(expression(Abundance~of~italic(Pheidole)~spp.))
	par(ps=8)
	mtext(ynames, side=2, line=1.6)

	#Paratrechina
	means<-(t(tapply(rawdata$Lpara,list(rawdata$sampling),mean)))
	stdev<-(t(tapply(rawdata$Lpara,list(rawdata$sampling),sd)))
	nsampling<-t(table(rawdata$sampling))
	ciw   <- qt(0.95, nsampling) * stdev / sqrt(nsampling)
	xscale<-(t(tapply(rawdata$sampling,list(rawdata$sampling),mean)))
	barplot2(means, plot.ci=TRUE,col=ParaCol, ci.l= means - ciw, ci.u= means + ciw,ci.color="black",axisnames = F, axes=F, xlab="", ylab="", ci.lty = "solid", ci.lwd = 1)
	par(ps=8)
	mtext("Dry Season            Rainy season", side=1, line=0.25)

	par(ps=9)
	axis(side = 2, tcl = 0.25, las = 1, labels = F)
	axis(side = 2, lwd = 0, las = 1, labels = T, line = -0.5)

	ynames<-c(expression(Abundance~of~italic(P.~longicornis)))
	par(ps=8)
	mtext(ynames, side=2, line=1.6)

	#Monomorium
	means<-(t(tapply(rawdata$Lmono,list(rawdata$sampling),mean)))
	stdev<-(t(tapply(rawdata$Lmono,list(rawdata$sampling),sd)))
	nsampling<-t(table(rawdata$sampling))
	ciw   <- qt(0.95, nsampling) * stdev / sqrt(nsampling)
	xscale<-(t(tapply(rawdata$sampling,list(rawdata$sampling),mean)))
	barplot2(means, plot.ci=TRUE,col=MonoCol, ci.l= means - ciw, ci.u= means + ciw,ci.color="black",axisnames = F, axes=F, xlab="", ylab="", ci.lty = "solid", ci.lwd = 1)
	par(ps=8)
	mtext("Dry Season            Rainy season", side=1, line=0.25)


	par(ps=9)
	axis(side = 2, tcl = 0.25, las = 1, labels = F)
	axis(side = 2, lwd = 0, las = 1, labels = T, line = -0.5)

	ynames<-c(expression(Abundance~of~italic(Monomorium)~spp.))
	par(ps=8)
	mtext(ynames, side=2, line=1.6)
	#Uncomment the code line below if you want to save image in .tiff
	#dev.off()

#------------------------------------------------------------------------------#
# End of Figure S1.
}
#------------------------------------------------------------------------------#

#------------------------------------------------------------------------------#
{# Figure S2
#------------------------------------------------------------------------------#

	#------------------------------------------------------------------------------#
	# Load the rawdata file
	#------------------------------------------------------------------------------#
	# !!! Set your own work directory that contains the rawdata file
	setwd("D:/Mes donnees/Dropbox/Articles/Revised/Dassou-et-al")
	rawdata<-read.table("rawdata.txt",h=T,sep="")

	library(lattice)
	
	#------------------------------------------------------------------------------#
	# All codominace cases (i.e. DomSpecies > 6) are turned 
	# into one category
	for (i in 1:length(rawdata[,1]))
	{
	  if (rawdata$DomSpecies[i] > 7)
	   rawdata$DomSpecies[i]<-7
	}

	Dominance<-rawdata
	Dominance$DomSpecies<-as.factor(Dominance$DomSpecies)

	#Modify DomSpecies factor values
	Dominance$DomSpecies<-as.factor(Dominance$DomSpecies)
	library(plyr)
	Dominance$DomSpecies<-revalue(Dominance$DomSpecies, c("5"="5", "1"="2", "2"="1", "3"="0", "4"="3", "6"="4", "0"="6", "7"="7"))
	Dominance$DomSpecies
	Dominance$DomSpecies<-as.numeric(as.character(Dominance$DomSpecies))

	#We plot only for the dry season (sampling = 2) 
	subrawdata<-subset(Dominance, Dominance$sampling == 2)

	coldom=c("indianred","mistyrose1","plum3","azure2","coral","white","darkseagreen1","gray")
	coldomL=c("indianred","mistyrose1","plum3","azure2","coral","darkseagreen1","white","gray")
	breakkeydom<-seq(0,7, by=1)
	Valuesdom<-c("Axi","Cam","Tet","Phe","Par","Mon","Emp","Cod")
	breakmapdom<-seq(-0.5,7.5, by=1)
	vectornames<-c("1","2","3","4","5","6","7","8","9","10","11","12","13","14","15","16","17","18","19", "20")
	stripcol<-"white"
	
	#Uncomment the code line below if you want to save image in .tiff
	#tiff("FigS2Unscaled.tiff", width=6.5, height=8.5, units="in", res=600)
	levelplot(DomSpecies~x*y|plot, xlab="", ylab= "", col.regions=coldom, colorkey=list(space="top",at=breakmapdom,col=coldomL,labels=list(at=breakkeydom,labels=Valuesdom), width=1, height=0.5),strip = strip.custom(strip.names =  FALSE, factor.levels = vectornames), par.strip.text = list(cex = 0.75),  par.settings = list(strip.background=list(col=stripcol)),scales = list(x = list(alternating=c(1,0), cex = 0.6, at=c(0, 300, 600, 900, 1200), labels=c("1.2","", "6.0", "", "10.8"), tck=c(0.25,0)), y = list(alternating=c(0,1), cex = 0.6, at=c(0, 300, 600, 900, 1200), labels=c("1.2","", "6.0", "", "10.8"), tck= c(0.25,0))), data=subrawdata, aspect.ratio=1)
	#Uncomment the code line below if you want to save image in .tiff
	#dev.off()
#------------------------------------------------------------------------------#
# End of Figure S2.
}
#------------------------------------------------------------------------------#