##Supplement 1##
##sampling model of the species-time-area relationship (STAR) 

##Written for R v2.6.2, 2008

#Daniel J. McGlinn
#117 Life Sciences East
#Oklahoma State University
#Stillwater, OK, 74074
#daniel.mcglinn@okstate.edu

##################################################
##create the nine relative abundance distributions (RADs)##
S<-800 #specify size of the species pool
##even RAD##
p1<-rep(1/S,S)
#############
##now generate the lognormals as averages of random lognormal distributions##
reps<-1000 ##number of replictates to average over
vars<-c(1,2,4) ##variances to use for the 3 LOGN RADs
X<-array(NA,dim=c(S,reps,length(vars))) #create holding array
for(k in 1:length(vars)){ #loop through # of variances
 for(i in 1:1000){ #loop through replicates
  X[,i,k]<-sort(rlnorm(S,0,vars[k]),decreasing=TRUE) #populate array with sorted values from largest to smallest
 }
}
X.m<-array(NA,dim=c(S,length(vars))) #create holding array
for(k in 1:length(vars)) 
 X.m[,k]<-sort(apply(X[,,k],1,mean),decreasing=TRUE)##calcualte mean and then order each distribution from largest to smallest
P.m<-array(NA,dim=c(S,length(vars))) #holding arrary for RAD
for(k in 1:length(vars))
 P.m[,k]<-X.m[,k]/sum(X.m[,k]) #now convert the abundance distribution into relative abundance distributions
##########################
##create the uneven RAD##
prop<-.99  #proportion of species pool dominated by one species
p5<-rep(0,S) #create holding vector
p5[1]<-prop #domiant species set
p5[-1]<-rep((1-prop)/(S-1),S-1) #the remaining S-1 species set
########################
##geometric series##
geom <- function(a1, k, S){
 ai <- a1*k^((1:S)-1)
 pi <- ai/sum(ai)
 pi
}
p6<-geom(10,.9,S)
##################
##broken stick##
brok <- function(abar, S){
 ai <- rep(NA,S)
 for (i in 1:S)
  ai[i] <- abar*sum(1/i:S)
 pi <- ai/sum(ai)
 pi
}
p7<-brok(10,S)
##################
##Zipf##
zipf <- function(a1,gamma,S){
 ai <- a1*(1:S)^-gamma
 pi <- ai / sum(ai)
 pi
}
p8<-zipf(10,1.3,S)
###################
##Zipf-Mandelbrot##
zipf.man <- function(a1,gamma,beta,S){
 ai <- a1*((1:S)+beta)^-gamma
 pi <- ai / sum(ai)
 pi
}
p9<-zipf.man(10,1.3,100,S)
#########################
P.mat<-cbind(p1,P.m,p5,p6,p7,p8,p9) #bind all the RAD vectors together into one matrix

##alternative binding for just the RADs that were primarily discussed in manuscript##
P.mat<-cbind(p1,P.m[,1],p6,P.m[,3],p5) #bind all the RAD vectors together into one matrix

###################################################
##Calculate Pielou's evenness (E) for each RAD##

E<-rep(0,dim(P.mat)[2]) #create a holder variable
for(i in 1:dim(P.mat)[2])
 for(ii in 1:S)
  E[i]<-E[i] -(P.mat[ii,i]*log(P.mat[ii,i])) #this effectively calculates the Shannon information (H')
E<-E/log(S) #scale the Shannon information by the maximum information
E ##values should be very similar to those in Table 1. unless changes were made above

####################################################
##RADs plot## 
##Figure C1##

xlabel<-"species rank"
ylabel<-expression("log probability")
plot(1:S,log10(P.mat[,1]),xlab=xlabel,ylab=ylabel,ylim=c(-12,log10(max(P.mat))),pch=19,type='n',lwd=2,frame.plot=F)
for(i in 1:dim(P.mat)[2])
 lines(1:S,log10(P.mat[,i]),col=i,lty=i,lwd=3)
###the next lines generate a legend### 
legend("bottomleft",c(expression("even: "*italic(p[1])*"..."*italic(p[S[p]])==frac(1,italic(S[p]))),
expression("log normal: "*italic(p)*" ~ "*"LOGN(0,1)"),
expression("geometric("*italic(k)*" = 0.9): "*italic(p[i])*" = "*frac(italic(k)*(1-italic(k))^(italic(i)-1),1-(1-italic(k))^italic(S[p]))),
expression("log normal: "*italic(p)*" ~ "*"LOGN(0,4)"),
expression("uneven: "*italic(p[1])*" = 0.99, "*italic(p[2])*"..."*italic(p[S[p]])*" = "*frac(0.01,italic(S[p])-1))),
lty=1:5,lwd=rep(2,5),col=1:5,bty='n',cex=1)

###########################################
##to create a cumulative probability plot##
p.c<-function(probs){
 p.cum<-probs[1]
 for (i in 1:(length(probs)-1)){
 p.cum[i+1]<-p.cum[i]+probs[i+1]
 }
p.cum 
}
##############################
#cumulative probability plots#
plot(1:S,p.c(P.mat[,1]),ylab="summed probability",ylim=c(0,1),main='even',type='l')
for(i in 2:dim(P.mat)[2])
 points(1:S,p.c(P.mat[,i]),col=i,lty=i,type='l',lwd=2)
legend('bottomright',legend=c(1:dim(P.mat)[2]),col=1:dim(P.mat)[2],lty=1:dim(P.mat)[2],lwd=rep(2,9))
############################################
##Set STAR parameters##

##to make scale increase multiplicatively
ncurves<- 15 #number of curves to calculate
areas<-1 #first area is at 1
times<-1 #first time is at 1
for(i in 2:ncurves){
 areas[i]<-areas[i-1]*2 #each sequential area is 2x a large as the first
 times[i]<-times[i-1]*2
}

R<-c(0,.01,.1,.5,1) #the values for the replacement rate

ExpS<-array(0,dim=c(ncurves,ncurves,length(R),dim(P.mat)[2])) #holder for expected richness

Z<-array(0,dim=c(ncurves,ncurves,length(R),dim(P.mat)[2])) #holder for the slope of the log-log SAR
W<-array(0,dim=c(ncurves,ncurves,length(R),dim(P.mat)[2])) #holder for the slope of the log-log STR
#################################################

##Calculate STAR##

for (k in 1:dim(P.mat)[2]){ #to loop through the RADs
 for (m in 1:length(R)){  #to loop through the values for the replacement rate
  for (i in 1:length(areas)){ #to loop through the areas
   for(ii in 1:length(times)){ #to loop through the times
    a<-areas[i]
    t<-times[ii]
    r<-R[m]
    ssums<-0 #set these three indices to zero before summing (see below)
    zsums<-0
    wsums<-0
    for (j in 1:S){ #now we will sum across the S species in the pool 
     ssums<-ssums+(1-P.mat[j,k])^(a+a*r*(t-1))
     zsums<-zsums+(1-P.mat[j,k])^(a+a*r*(t-1))*log(1-P.mat[j,k])*(1+r*t-r)
     wsums<-wsums+(1-P.mat[j,k])^(a+a*r*(t-1))*log(1-P.mat[j,k])*(a*r)
    }
    ExpS[i,ii,m,k]<-S-ssums
    Z[i,ii,m,k]<-(-a*zsums)/ExpS[i,ii,m,k]
    W[i,ii,m,k]<-(-t*wsums)/ExpS[i,ii,m,k]
}}}}

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

##Calculate time-by-area interaction##
#this does not have to be run seperately of the previous code but it increases
#computation time so I have seperated here so that you can chose to calcualate U or not

U<-array(0,dim=c(length(areas),length(times),length(R),dim(P.mat)[2]))

for (k in 1:dim(P.mat)[2]){
 for(m in 1:length(R)){
  for (i in 1:length(areas)){
   for(ii in 1:length(times)){
    a<-areas[i]
    t<-times[ii]
    r<-R[m]
    u1<-0 
    u2<-0 
    for (j in 1:S){
     u1<- u1 + ((1-P.mat[j,k])^(a+a*r*(t-1))*log(1-P.mat[j,k]))
     u2<- u2 + ((1-P.mat[j,k])^(a+a*r*(t-1))*(log(1-P.mat[j,k])^2)*(a*r))
    }
    U[i,ii,m,k]<- t*(((-a*r*u1 - u2*(a*(1+r*(t-1))))* ExpS[i,ii,m,k]^-1) + ((ExpS[i,ii,m,k]^-2*a*r*u1)*((-a*(1+r*(t-1)))*u1)))  ##It's a pretty long equation!##
}}}}
#########################################

###visuallizing the STAR##########

op <- par(no.readonly = TRUE) #save the default graphical parameters

#First chose the RAD for which you would like to see the STAR for

k<-3 #in this case we have chosen the third RAD
cls<- colorRampPalette(c("blue","skyblue","lightblue"),space="rgb")(ncurves)

##SAR and STR curves
par(mfrow=c(2,length(R))) 
for (m in 1:length(R)){ 
 plot(log(areas),log(ExpS[,1,m,k]),ylim=c(min(log(ExpS[,,,k])),max(log(ExpS[,,,k]))),type='n',xlab='ln area',ylab='ln richness')
 for(i in 1:length(times)){
  points(log(areas),log(ExpS[,i,m,k]),type='l',col=cls[i])
 }
}
for (m in 1:length(R)){
 plot(log(times),log(ExpS[1,,m,k]),type='n',ylim=c(min(log(ExpS[,,,k])),max(log(ExpS[,,,k]))),xlab='ln time',ylab='ln richness')
 for(i in 1:length(areas)){
  points(log(times),log(ExpS[i,,m,k]),type='l',col=cls[i])
 }
}

###W as a function of ln area AND Z as a function of ln time
par(mfrow=c(2,length(R)))
for (m in 1:length(R)){
 plot(log(areas),W[,1,m,k],type='n',ylim=c(min(W[,,,k]),max(W[,,,k])),xlab='ln area',ylab='w')
 for(i in 1:length(times)){
  points(log(areas),W[,i,m,k],type='l',col=cls[i])
 }
}
for (m in 1:length(R)){
 plot(log(times),Z[1,,m,k],ylim=c(min(Z[,,,k]),max(Z[,,,k])),type='n',xlab='ln time',ylab='z')
 for(i in 1:length(areas)){
  points(log(times),Z[i,,m,k],type='l',col=cls[i])
 }
}

###U as a function of area and then time##########
par(mfrow=c(2,length(R)))
for (m in 1:length(R)){
 plot(log(areas),U[,1,m,k],ylim=c(min(U[,,,k]),max(U[,,,k])),type='n',xlab='',ylab='')
 for(i in 1:length(times)){
  points(log(areas),U[,i,m,k],type='l',col=cls[i])
 }
abline(h=0,col='grey',lty=2)
}
for (m in 1:length(R)){
 plot(log(times),U[1,,m,k],ylim=c(min(U[,,,k]),max(U[,,,k])),type='n',xlab='',ylab='')
 for(i in 1:length(areas)){
  points(log(times),U[i,,m,k],type='l',col=cls[i])
 }
abline(h=0,col='grey',lty=2)
}

###########################################################
#####color versions for online appendices########################
##to elminate spaces between plots use the "mar" command within "par()"
library(dichromat)
cls<- colorRampPalette(colorschemes$BrowntoBlue.10,space="Lab")(ncurves) #set up colors template
##SARs
##Figure C2
linewd<-1.5 #width of lines
par(mfrow=c(dim(P.mat)[2],length(R))) #set graphical space 
#the replacement rate will increase to the right across the columns
#evenness will decrease downward across the rows
for(k in dim(P.mat)[2]:1){
 for(m in 1:length(R)){
  par(mar=rep(2,4))
  plot(log10(areas),log10(ExpS[,1,m,k]),xlim=c(0,log10(max(areas))),ylim=c(min(log10(ExpS)),ceiling(max(log10(ExpS)))),type='n',xlab='',ylab='',frame.plot=F)
  for(i in 1:ncurves)
   points(log10(areas),log10(ExpS[,i,m,k]),type='l',lwd=linewd,col=cls[i])
  abline(h=log10(S),col='grey',lty=2,lwd=linewd)
}}
##STRs
##Figure C3
linewd<-1.5
par(mfrow=c(dim(P.mat)[2],length(R)))
for(k in dim(P.mat)[2]:1){
 for(m in 1:length(R)){
  par(mar=rep(2,4))
  plot(log10(times),log10(ExpS[1,,m,k]),xlim=c(0,log10(max(areas))),ylim=c(min(log10(ExpS)),ceiling(max(log10(ExpS)))),type='n',xlab='',ylab='',frame.plot=F)
  for(i in 1:ncurves)
   points(log10(times),log10(ExpS[i,,m,k]),type='l',lwd=linewd,col=cls[i])
  abline(h=log10(S),col='grey',lty=2,lwd=linewd)
}}
###############################################################
##U as a function of accumulated richess (ExpS)
##hold area constant and vary time
##Figure C4##
linewd<-1.5
par(mfrow=c(dim(P.mat)[2],length(R)))
for(k in dim(P.mat)[2]:1){
 for(m in 1:length(R)){
  par(mar=rep(2,4))
  plot(log10(ExpS[1,,m,k]),U[1,,m,k],xlim=c(0,log10(S)),ylim=c(min(U),max(U)),type='n',xlab='',ylab='',frame.plot=F)
  abline(h=0,col='grey',lty=2,lwd=linewd)
  abline(v=S,col='grey',lty=2,lwd=linewd)
  for(i in 1:ncurves)
   points(log10(ExpS[i,,m,k]),U[i,,m,k],type='l',lwd=linewd,col=cls[i])
}}
#########################################################

par(op) ##reset the graphical parameters to default or simply close the active graphic

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