#A simple model to show how random walks of clades with different inherent turnover rates can lead to decreasing turnover over time
#This doesn't accurately model speciation/extinction and the fossil record certainly contains non-random events due to extinctions, evolutionary innovations, etc.

null_model <- function() {
  
  #Creates an empty matrix with 12 rows (for the 12 taxa) and 100 columns (for the 100 time-steps)
  taxa<-matrix(nrow=12,ncol=100)
  
  #Starting position - each taxon has a "diversity" of 10
  taxa[,1]<-rep(10,12)
  
  #Runs 99 iterations of the diversification model
  #There are three groups of taxa, four with a standard deviation (step size) of 2.5, four with a step size of 1, and four with a step size of 0.25
  #All taxa have a mean change of 0 (i.e. no directionality, an unbiased random walk)
  #In each iteration, a turnover step size is randomly generated for each taxon from the normal distribution
  #That step size is added to the previous diversity (time i) to generate the new diversity (time j)
  #the "which(taxa[,i]..." part just tells R to stop adding the step sizes once a taxon has reached 0 ("extinction")
  
  for (i in 1:99) {
    j=i+1
    turnover<-c(rnorm(4,0,2.5),rnorm(4,0,1),rnorm(4,0,0.25))
    taxa[which(taxa[,i]>0),j]<-taxa[which(taxa[,i]>0),i]+turnover[which(taxa[,i]>0)]
  }
  
  
  #Creates a plot window with two panels, and modifies the margin sizes
  par(mfcol=c(2,1),mar=c(4,4,1,2))
  
  #plots the diversity time-series for taxon 1
  plot(seq(1:length(taxa[1,])),taxa[1,],type="l",col="red",lwd=2,ylim=c(0,max(taxa,na.rm=T)),xlab="Time step",ylab="Diversity")
  
  #the next set of commands adds the diversity lines for the remaining 11 taxa
  lines(seq(1:length(taxa[2,])),taxa[2,],col="red",lwd=2)
  lines(seq(1:length(taxa[3,])),taxa[3,],col="red",lwd=2)
  lines(seq(1:length(taxa[4,])),taxa[4,],col="red",lwd=2)
  
  lines(seq(1:length(taxa[5,])),taxa[5,],col="orange",lwd=2)
  lines(seq(1:length(taxa[6,])),taxa[6,],col="orange",lwd=2)
  lines(seq(1:length(taxa[7,])),taxa[7,],col="orange",lwd=2)
  lines(seq(1:length(taxa[8,])),taxa[8,],col="orange",lwd=2)
  
  lines(seq(1:length(taxa[9,])),taxa[9,],col="blue",lwd=2)
  lines(seq(1:length(taxa[10,])),taxa[10,],col="blue",lwd=2)
  lines(seq(1:length(taxa[11,])),taxa[11,],col="blue",lwd=2)
  lines(seq(1:length(taxa[12,])),taxa[12,],col="blue",lwd=2)
  
  #calculates the average "turnover" (absolute value of the step size change) in each time step
  turnover.avg<-apply(abs(apply(taxa,1,diff)),1,function(x) mean(x,na.rm=T))
  
  #plots the turnover and adds a trendline
  plot(turnover.avg,xlab="Time step",ylab="Turnover Rate")
  abline(lm(turnover.avg~seq(1:99)),col="red")
  
  #End
  
}