# Toward ecologically explicit null models of nestedness
# J.E. Moore and R. K. Swihart
# Oecologia 152:763-777

# R Code for calculating a null presence-absence matrix of expected occurrence probabilities, based on a given null model
# Null model is based on a prior ecological hypothesis (e.g., species-area relationship, random placement model, etc)
# Inputs are expected value data (expected marginal [row/column] sums)
# Input data are expected marginals for songbirds in the Upper Wabash River Basin in northern Indiana, as presented in Moore and Swihart (2007)

# Created 7 November 2006, by Jeffrey E. Moore
# Modified and added to website 21 March 2008

# Contact jemoore@duke.edu





#### INPUT EXPECTED-VALUE DATA FOR CONSTRUCTING RANDOM MATRICES ####

# Input data are the expected richness (row) and incidence (col) totals, based on a priori model
# These expectations can be input here, generated here (by fitting a model here), or read in as vectors

# expected richness in each site
rich <-  c(6.997,8.071,6.573,4.825,7.977,5.553,8.319,8.787,8.598,10.270,11.082,5.716,7.243,7.873,6.744,6.863,7.650,2.364,8.554,6.202, 9.855,8.224,6.728,9.299,6.711,8.302,9.244,6.571,4.260,9.133)

# expected incidence for each species (# sites where it will occur)
inc <- c(28.135,22.544,1.953,4.562,17.3,2.863,16.317,12.671,17.299,1,29.406,25.538,28.618,7.537,8.845)





#### CALCULATE INITIAL PARAMETERS ####


Nsites <- length(rich)		# number of sites sampled
Nspp <- length(inc)		# number of species observed in data

JF <- matrix(NA,Nsites,Nspp)	# empty site x spp matrix to hold joint-frequencies

for(i in 1:Nsites){		# calculate joint frequencies
	for(j in 1:Nspp){
		JF[i,j] <- rich[i]*inc[j]
	}
}

EF <- JF/sum(JF)*sum(rich)	# expected frequencies (=element probs x expected fill)

probDens <- EF/sum(EF)		# matrix sums to 1





#### ROUND 1 OF RE-ALLOCATING "EXTRA PROBABILITIES" ####


extra <- matrix(NA,Nsites,Nspp)

for(i in 1:Nsites){		# calculate "extra probability" in each element (amount in EF that is > 1)
	for(j in 1:Nspp){
		extra[i,j] <- max(0,(EF[i,j]-1))
	}
}

R.extra1 <- rowSums(extra)	# marginal sums of "extra probability"
C.extra1 <- colSums(extra)


idx <- matrix(NA,Nsites,Nspp)

for(i in 1:Nsites){		# index elements where original expectation < 1
	for(j in 1:Nspp){
		if(extra[i,j] > 0)
		idx[i,j] <- 0
		else idx[i,j] <- 1
	}
}


pDens1 <- matrix(NA,Nsites,Nspp)	# empty matrix to hold prob densities where EF < 1
pDens1.Rrescale <- pDens1
pDens1.Crescale <- pDens1		# empty matrices to hold rescaled row and col densities where EF < 1

for(i in 1:Nsites){		# probability densities (from probDens) where EF < 1 (put zeros elsewhere in this matrix)
	for(j in 1:Nspp){
		if(idx[i,j]==1)
		pDens1[i,j] <- probDens[i,j]
		else
		pDens1[i,j] <- 0
	}
}

pD1.Rmarg <- rowSums(pDens1)	# marginal sums of densities from elements where EF < 1
pD1.Cmarg <- colSums(pDens1)

for(i in 1:Nsites){			# calculate proportion of excess in row i or col j that will go in each element 
	for(j in 1:Nspp){
		pDens1.Rrescale[i,j] <- pDens1[i,j]/pD1.Rmarg[i]
		pDens1.Crescale[i,j] <- pDens1[i,j]/pD1.Cmarg[j]
	}
}


R.alloc1 <- matrix(NA,Nsites,Nspp)	# empty matrices to hold added probs
C.alloc1 <- R.alloc1

for (i in 1:Nsites){		# allocate excess probability across elements
	for (j in 1:Nspp){
		R.alloc1[i,j] <- pDens1.Rrescale[i,j] * R.extra1[i]
		C.alloc1[i,j] <- pDens1.Crescale[i,j] * C.extra1[j]
	}
}


EF2.temp <- matrix(NA,Nsites,Nspp)	# empty 'temporary matrix' to hold new frequencies, 
					# following double allocation of excess across rows AND columns

for(i in 1:Nsites){			# calcuate the temporary matrix of new frequencies
	for(j in 1:Nspp){
		if(EF[i,j]>1)
		EF2.temp[i,j] <- 1
		else
		EF2.temp[i,j] <- EF[i,j] + R.alloc1[i,j] + C.alloc1[i,j]
	}
}


rFreqDiff1 <- rowSums(EF2.temp)-rowSums(EF)	# calculate marginal sums of excess probability in the matrix, following 2x allocation
cFreqDiff1 <- colSums(EF2.temp)-colSums(EF)


removeProb1 <- matrix(NA,Nsites,Nspp)		# calculate amount of the new excess to remove from each element
						# do this proportionally to the product of marginal sums of the excess
for(i in 1:Nsites){
	for(j in 1:Nspp){
		removeProb1[i,j] <- rFreqDiff1[i]*cFreqDiff1[j]/sum(rFreqDiff1)
	}
}

EF2 <- EF2.temp - removeProb1		# remove the excess from each element, to get the updated expected frequency matrix

max(EF2)	# evaluate the new matrix
		# if max element value = 1, we're done.  We have the new expected frequency matrix
		# if max element value > 1, we need to do another iteration
		# in this example, max = 1.18, so we need to continue


probDens2 <- EF2/sum(EF2)		# matrix sums to 1





#### ROUND 2 OF RE-ALLOCATING "EXTRA PROBABILITIES" ####

# This round is MOSTLY identical to above code.
# Exception 1: object names have been changed, so that most objects in this step refer to objects output from previous step
# Exception 2 occurs in loop involving creation of "idx" object



extra2 <- matrix(NA,Nsites,Nspp)

for(i in 1:Nsites){		# calculate "extra probability" in each element (amount in EF > 1)
	for(j in 1:Nspp){
		extra2[i,j] <- max(0,(EF2[i,j]-1))
	}
}


R.extra2 <- rowSums(extra2)	# marginal sums of extra frequency
C.extra2 <- colSums(extra2)


idx2 <- matrix(NA,Nsites,Nspp)

for(i in 1:Nsites){		# index values in which original expectation < 1
	for(j in 1:Nspp){
		if(extra2[i,j] > 0 | idx[i,j]==0)	# Diff't from first iteration, with added argument "or idx[i,j]==0"
		idx2[i,j] <- 0
		else idx2[i,j] <- 1
	}
}


pDens2 <- matrix(NA,Nsites,Nspp)	# empty matrix to hold prob densities where EF < 1
pDens2.Rrescale <- pDens2
pDens2.Crescale <- pDens2		# matrices to hold rescaled row and col densities where EF < 1

for(i in 1:Nsites){
	for(j in 1:Nspp){
		if(idx2[i,j]==1)
		pDens2[i,j] <- probDens2[i,j]
		else
		pDens2[i,j] <- 0
	}
}


pD2.Rmarg <- rowSums(pDens2)
pD2.Cmarg <- colSums(pDens2)

for(i in 1:Nsites){			# calculate proportion of excess in row i or col j that will go in each element 
	for(j in 1:Nspp){
		pDens2.Rrescale[i,j] <- pDens2[i,j]/pD2.Rmarg[i]
		pDens2.Crescale[i,j] <- pDens2[i,j]/pD2.Cmarg[j]
	}
}


R.alloc2 <- matrix(NA,Nsites,Nspp)	# empty matrices to hold added probs
C.alloc2 <- R.alloc2

for (i in 1:Nsites){
	for (j in 1:Nspp){
		R.alloc2[i,j] <- pDens2.Rrescale[i,j] * R.extra2[i]
		C.alloc2[i,j] <- pDens2.Crescale[i,j] * C.extra2[j]
	}
}

EF3.temp <- matrix(NA,Nsites,Nspp)	# matrix of new frequencies, before taking some probability "back out"

for(i in 1:Nsites){
	for(j in 1:Nspp){
		if(EF2[i,j]>1)
		EF3.temp[i,j] <- 1
		else
		EF3.temp[i,j] <- EF2[i,j] + R.alloc2[i,j] + C.alloc2[i,j]
	}
}


rFreqDiff2 <- rowSums(EF3.temp)-rowSums(EF2)
cFreqDiff2 <- colSums(EF3.temp)-colSums(EF2)

removeProb2 <- matrix(NA,Nsites,Nspp)

for(i in 1:Nsites){
	for(j in 1:Nspp){
		removeProb2[i,j] <- rFreqDiff2[i]*cFreqDiff2[j]/sum(rFreqDiff2)
	}
}

EF3 <- EF3.temp - removeProb2

max(EF3)


probDens3 <- EF3/sum(EF3)		# matrix sums to 1





#### ROUND 3 OF RE-ALLOCATING "EXTRA PROBABILITIES" ####

# repeat as in ROUND 2
# For this dataset, four iterations were required to get all expected freq elements <= 1