# R functions  used in the manuscript “A Method for Estimating Abundance of Mobile Populations Using 
# Telemetry and Counts of Unmarked Individuals” by Clement et al. 

# Load the functions contained in the file “coincidence_functions.R” into R. 
# To do this, save “coincidence_functions.R” to your computer. 
# Then open R and select File->Source R code, and select the file “coincidence_functions.R”.

# Contents
# 2          	R functions
# 2.8                      	field.N.foo
# 2.9                      	field.boot.foo
# 2.1                      	sim.data.foo
# 2.2                      	summ.data.foo
# 2.3                      	est.N.foo
# 2.4                      	boot.data.foo
# 2.5                      	fager.foo
# 2.6                      	move.foo
# 2.7                      	pow.analysis




#####	2.8	field.N.foo	#######################################################################

#----------	ESTIMATE ABUNDANCE FROM FIELD DATA  --------------------------------------------

field.N.foo	<- function(bat.loc, grp.size) {


   days		<- dim(bat.loc)[1]		# number of days in data
   num.bats	<- dim(bat.loc)[2]		# number of tagged bats


      # now estimate abundance

         # first, roost size
      mu		<- 1:num.bats
      for (nb in 1:num.bats) {
	 denom			<- if(sum(bat.loc[,nb,], na.rm=T)==0) NA else sum(bat.loc[,nb,], na.rm=T)
         mu[nb]			<- sum(bat.loc[,nb,] * grp.size[,], na.rm=T)/denom
      }

         # second, coincidence rate
      num.pairs		<- choose(num.bats,2)
      coin		<- 1:num.pairs
      c3		<- 0
      for (c in 1:(num.bats-1)) {
         for (c2 in (c+1):num.bats) {
            c3			<- c3+1
		# I had to add the extra tests to deal with bats in an unknown loc.
            coin[c3]	<- mean('&'('&'(max.col(bat.loc[,c,], ties.method="f") == max.col(bat.loc[,c2,], ties.method="f"), max.col(bat.loc[,c,], ties.method="f") == max.col(bat.loc[,c,], ties.method="l")), max.col(bat.loc[,c2,], ties.method="f") == max.col(bat.loc[,c2,], ties.method="l")), na.rm=T)
         }
      }

   field.N	<- (mean(mu)-1)/mean(coin)+1	


   field.N
}	# ends the function
#------------------------------------------------------------------------------------------------------------------





#####	 2.9		field.boot.foo   ###############################################################

#----------	BOOTSTRAP FUNCTION TO ESTIMATE PRECISION OF FIELD DATA  --------------------------------------------

field.boot.foo	<- function(bat.loc, grp.size, strap=1000) {


# bootstrap the tagged bats to get an estimate of precision

   N.strap		<- 1:strap

   days		<- dim(bat.loc)[1]		# number of days in data
   num.bats	<- dim(bat.loc)[2]		# number of tagged bats

   for (s in 1:strap) {

      # bootstrap data
      strap.days	<- sample(1:days, days, replace = T)
      strap.loc		<- array(bat.loc[strap.days,,], dim=dim(bat.loc))
   
      # now estimate abundance with new data
      N.strap[s]	<- field.N.foo(strap.loc,grp.size[strap.days,]) 

   }	# ends the strap loop

   N.strap
}	# ends the function
#------------------------------------------------------------------------------------------------------------------




##### 2.1 	sim.data.foo	######################################################################

####################################################
####################################################
#==================================================
# function to simulate bats and estimate abundance 
#==================================================
#
#	days: number of days of surveys;				colony: size of bat colony
#	trees: number of potential roosting locations		num.bats: number of bats that are radio-tagged
#	cliques: number of different social units		attract: amount of attraction between bats in the same social unit
#	repel: attraction bt bats in diff social units (negative means attraction)
#	base: level of attraction if no social attraction		stay: attraction to the roost used yesterday
####################################################
####################################################


sim.data.foo	<- function(days=20, colony=100, trees=20, attract=800, base=20, cliques=1, repel=0, stay=0, BI=0, method="log") {


#----------	SIMULATE DATA  --------------------------------------------

# create social groups of (nearly) the same size
# 	So if 101 animals and 2 groups, then 51 and 50 animals in the groups 
extra 			<- cliques*ceiling(colony/cliques) - colony	
clique.ind		<- matrix(c(1:colony,rep(0,extra)), ncol=cliques, nrow=ceiling(colony/cliques), byrow=T)


###--------------------- assign all bats to groups ------------------------------

# this array will track the location of each animal at each time
bat.loc			<- array(0, dim=c(BI+days,colony,trees))

# this array contains a score for the attraction to each location at each time
bat.probs		<- base + bat.loc				# prob proportional to available trees



for (d in 1:(BI+days)) {

   # mix up the colony every day so it is not the same bat assigned first every day. 
   # mixing needed bc otherwise bats tend to end up with bats near them in the list.
   mix.colony		<- sample(1:colony, colony)

   # the two lines below make a bat more likely to select the same group as yesterday
   # i.e. which group was the first bat in yesterday
   if (d>1)	AR.i <- which(bat.loc[d-1,mix.colony[1],]==1)					
   # then adjust bat.probs accordingly
   if (d>1)	bat.probs[d,mix.colony[1],AR.i] <- stay+bat.probs[d,mix.colony[1],AR.i]		

   # assign the first animal to a group based on bat.probs score
   bat.loc[d,mix.colony[1],]		<- rmultinom(1, 1, bat.probs[d,mix.colony[1],])		

   for (c in 2:colony) {

   # calculate new probs to assign subsequent animals
   # prob depends on number of bats in same social group, or another social group
   # default attraction is logarithmic, so declining marginal value of more bats
   # cliq is which clique is the about-to-be-assigned bat is in
   cliq		<- ceiling(which(clique.ind==mix.colony[c])/dim(clique.ind)[1])
   # same.cliq is all the bats in that same cliq
   same.cliq		<- clique.ind[which(clique.ind[,cliq]>0),cliq]

# Section below adjusts bat probabilities by the locations of all the previous bats
#----------------------------------------------------------------
if (method == "log") {
   # LOGARITHMIC
      bat.probs[d,mix.colony[c],]	<- pmax(0, base+attract*log(colSums(bat.loc[d,same.cliq,])+1) - repel*log(colSums(bat.loc[d,-same.cliq,])+1))
}

if (method == "linear") {
   # LINEAR
        bat.probs[d,mix.colony[c],]	<- pmax(0, base+attract*colSums(bat.loc[d,same.cliq,]) - repel*colSums(bat.loc[d,-same.cliq,]))
}
#----------------------------------------------------------------


      # the two lines below make a bat more likely to select the same group as yesterday
      if (d>1) AR.i <- which(bat.loc[d-1,mix.colony[c],]==1)
      if (d>1) bat.probs[d,mix.colony[c],AR.i] <- stay+bat.probs[d,mix.colony[c],AR.i]

      # then assign the next bat based on bat.probs score
      bat.loc[d,mix.colony[c],]		<- rmultinom(1, 1, bat.probs[d,mix.colony[c],])		
  
      }	# end of colony loop

}	# end of Days loop


groups	<- apply(bat.loc[(BI+1):(BI+days),,], c(1,3),sum)				# this gives group sizes, excluding burn-in

  return(list(colony=colony, days=days, trees=trees, cliques=cliques, attract=attract, repel=repel, stay=stay, BI=BI, base=base, groups=groups, bat.loc=bat.loc))

}	# end sim.data function





##### 2.2		summ.data.foo	##########################################################

summ.data.foo	<- function(data) {

# here we get some descriptive statistics to summarize the groups

# how many groups are there? how large?
grp.sort	<- t(apply(data$groups,1,sort,decreasing=TRUE))
num.grps	<- max.col(-1*grp.sort, ties.method="first") - 1
	
mean.num.grps		<- mean(num.grps)
num.grps.range	<- quantile(num.grps, c(0.025,0.975))
mean.grp.size		<- colMeans(grp.sort)
grp.size.range	<- apply(grp.sort, 2, quantile, c(0.025,0.975))

  return(list(mean.num.grps = mean.num.grps, num.grps.range = num.grps.range, mean.grp.size =  mean.grp.size, grp.size.range = grp.size.range))

}	# end summ.data function





##### 2.3	est.N.foo		################################################################

est.N.foo		<- function(data, num.bats=4) {

days	<- data$days
BI	<- data$BI
bat.loc<- data$bat.loc	
###--------------------- assign radio tagged bats ------------------------------

# random values for assigning tagged bats to groups
tags		<- sample(1:data$colony, num.bats)

#----------	ESTIMATE ABUNDANCE FROM DATA  --------------------------------------------

   N.est	<- field.N.foo(bat.loc[(BI+1):(BI+days),tags,], data$groups)


#return(list(N = N.est, mu = mean(mu), c = mean(coin), num.bats=num.bats, tags=tags))
return(list(N = N.est, num.bats=num.bats, tags=tags))



}	# end est.N function



##### 2.4	boot.data.foo	#########################################################

#----------	ESTIMATE PRECISION FROM DATA  --------------------------------------------
# sim.data is output of sim.data.foo function and N.data is output of est.N.foo function
boot.data.foo		<- function(sim.data, N.data, strap=100) {

BI		<- sim.data$BI
days		<- sim.data$days
bat.loc	<- sim.data$bat.loc
groups		<- sim.data$groups
tags		<- N.data$tags
num.bats	<- N.data$num.bats

# N.strap will hold bootstrapped estimates of abundance
N.strap		<- rep(0, strap)

tagged.bats		<- bat.loc[(BI+1):(BI+days),tags,]
 

for (s in 1:strap) {
   # bootstrap data
   strap.days		<- sample(1:days, days, replace = T)
   strap.loc		<- tagged.bats[strap.days,,]
   
   # now estimate abundance with new data

   N.strap[s]	<- field.N.foo(strap.loc, groups[strap.days,])

}	# ends the strap loop

return(list(N.strap = N.strap))

}	# end boot.data function





##### 2.5	fager.foo	#################################################################

#----------	ESTIMATE FAGER INDEX  --------------------------------------------
#----------	ONLY WORKS WITH TWO SOCIAL SUB-UNITS

#	calculated days i&j in same roost divided by total days. 

fager.foo		<- function(sim.data) {

colony		<- sim.data$colony
bat.loc	<- sim.data$bat.loc
BI		<- sim.data$BI
days		<- sim.data$days
cliques	<- sim.data$cliques

fager	<- matrix(0,nrow=colony, ncol=colony)

for (c in 1:colony) {
   for (c2 in 1:colony) {
      #	max.col returns the group that bat c is in
      fager[c,c2]		<- mean(max.col(bat.loc[(BI+1):(BI+days),c,]) == max.col(bat.loc[(BI+1):(BI+days),c2,]))
   }
}

extra 			<- cliques*ceiling(colony/cliques) - colony
clique.ind		<- matrix(c(1:colony,rep(0,extra)), ncol=cliques, nrow=ceiling(colony/cliques), byrow=T)

# compare the cliques **ASSUMES TWO SOCIAL GROUPS**
c11	<-	c12	<- c21	<- c22	<- 1:nrow(clique.ind)
for (i in 1:nrow(clique.ind)) {
	c11[i]	<- mean(fager[clique.ind[-i,1],clique.ind[i,1]])	
	c12[i]	<- mean(fager[clique.ind[,1],clique.ind[i,2]])
	c21[i]	<- mean(fager[clique.ind[,2],clique.ind[i,1]])
	c22[i]	<- mean(fager[clique.ind[-i,2],clique.ind[i,2]])
}
mn.c11 <- mean(c11, na.rm=T);	mn.c22 <- mean(c22, na.rm=T)	
mn.c12 <- mean(c12, na.rm=T);	mn.c21 <- mean(c21, na.rm=T)
social.ratio	<- mean(c(mn.c11,mn.c22))/mn.c12


return(list(fager = fager, ratio=social.ratio, c11=mn.c11, c22=mn.c22, c12=mn.c12, c21=mn.c21))

}	# end fager function

#---------------------------------------------------------------------------------





##### 2.6	move.foo	###############################################################

#------------- 	DID THE ANIMAL MOVE LOCATIONS?	-----------------------------------------------

move.foo		<- function(sim.data) {

colony		<- sim.data$colony
days		<- sim.data$days
BI		<- sim.data$BI
bat.loc	<- sim.data$bat.loc

no.move	<- rep(0,days-1)			

for (d in (BI+2):(BI+days)) {
  for (c in 1:colony) {

      if (which(bat.loc[d,c,]==1) == which(bat.loc[d-1,c,]==1)) no.move[d-1] <- no.move[d-1]+1

   }
}	# end loops

1-mean(no.move/colony)

}	# end move.foo function







##### 2.7	pow.analysis	#######################################################

##############################################################################

#	POWER ANALYSIS

##############################################################################


pow.analysis	<- function(sites=1, days=20, control=100, treat=40, trees=20, num.bats=4, attract=800, reps=300, rand=300, base=20, cliques=1, repel=0, stay=0, BI=0, method="log") {

#sites=2; days=5; control=10; treat=5; trees=4; num.bats=3; attract=800; reps=1; rand=3
#base=20;cliques=1;repel=0;stay=0;BI=0;method="log"

# set up objects to hold estimates of N
N.rand.sites		<- 1:(2*sites)
N.diff			<- matrix(0, nrow=rand, ncol=reps)

# repeat the analysis 'reps' times
for (p in 1:reps) {

   # set up objects to hold data from multiple sites
   loc.data			<- array(0, dim=c(BI+days, num.bats, trees, 2*sites))	# loc at multiple sites, but only tagged bats
   group.data			<- array(0, dim=c(BI+days, trees, 2*sites))

   # create several colonies of different sizes
   colony				<- 1:(2*sites)
   colony[1:sites]			<- rmultinom(1, control*sites, rep(1, sites))
   colony[(1+sites):(2*sites)]		<- rmultinom(1, (control+treat)*sites, rep(1, sites))


   #-------- SIMULATE ANIMAL LOCATIONS AND GROUP SIZES AT ALL SITES  ----------------
   for (s in 1:(2*sites)) {	
     sim.out	<- sim.data.foo(days=days, colony=colony[s], trees=trees, attract=attract, base=base, 
cliques=cliques, repel=repel, stay=stay, BI=BI, method=method)

     N.out		<- est.N.foo(sim.out, num.bats=num.bats)

        loc.data[,,,s]			<- sim.out$bat.loc[(BI+1):(BI+days),N.out$tags,]
        group.data[,,s]		<- sim.out$groups

   } 	# end sites loop


   #----------	RANDOMIZE DATA, ESTIMATE ABUNDANCE -------------------------------------------

   for (r in 1:rand) {
      # prepare objects to hold randomized data
      group.data.rand	<- group.data	
      loc.data.rand	<- loc.data

      # generate index to data and randomize it
      mc			<- matrix(1:(2*sites), nrow=days, ncol=2*sites, byrow=T)	# site names
      ind			<- t(apply(mc, 1, sample))					# randomize site names

      # randomize the data using index above 
      if ( r>1) {					# don't randomize first day bc this is the original data
         for (d in 1:days) {
            for (s in 1:(2*sites)) {
               group.data.rand[d,,s]	<- group.data[d,,ind[d,s]]
               loc.data.rand[d,,,s]		<- loc.data[d,,,ind[d,s]]
            }
         }
      }

      #---------------------- now estimate abundance

      for (s in 1:(2*sites)) {  		
      
	   N.rand.sites[s]	<- field.N.foo(loc.data.rand[,,,s], group.data.rand[,,s])

      } 	# end the sites loop



      # difference in abundance bt treatment and control is test statistic
      N.diff[r,p]		<- mean(N.rand.sites[(1+sites):(2*sites)]) - mean(N.rand.sites[1:sites])

   } 	# end the rand loop
 
} 	# end the reps loop


# Extract the one-tail percentile from N.diff and only save that to save disc space.
ord		<- apply(N.diff, 2, order)		# order estimates from low to high
perc		<- apply(ord==1, 2, which)/rand	# where is the first estimate in that list? (1st is not randomized)
power		<- mean(perc>=0.95)			# how often is data among the largest (one tailed)
power.2tail	<- mean(perc<=0.025 | perc>=0.975)	# two tailed version

  return(list(power=power, power.2tail=power.2tail, sites=sites, control=control, treat=treat, days=days, trees=trees, num.bats=num.bats, attract=attract, base=base, reps=reps, rand=rand))


} 	# end PA function