##### read in main file
X = read.csv("RNR_Master_May2018.csv", header=TRUE)
colnames(X)[13]= "stage"
X = X[!(X$stage == ""),]

Z = read.csv("stage_numbered.csv", header=TRUE)
XZ = merge(X, Z, by="stage")

	# subset to only neccessary fields
Y = XZ[ , c("class", "order", "family", "genus", "period_updated", "Num", "Reef_Nonreef")]

	# Split the indeces of the dataframe (drop=TRUE means non-real combos of period and genus name are left out)
Z = split(seq(nrow(Y)) , list(Y$Num, Y$genus), drop=TRUE)


# loop over each genus and extract stats
Q = lapply(Z, function(i){
		# subset indeces to actual data
	z = Y[i, ]
		# remove the occs with NAs for reef_nonreef
	q = z[!is.na(z$Reef_Nonreef), ]

		# ratio of reef to non reef
	RnRF = nrow(q[ q$Reef_Nonreef=="reef", ] ) / (nrow(q[ q$Reef_Nonreef=="nonreef", ] ) + nrow(q[ q$Reef_Nonreef=="reef", ] ))

	return(list(period=as.character(unique(z$period_updated)), stage=as.character(unique(z$Num)), class=as.character(unique(z$class)), order=as.character(unique(z$order)), family=as.character(unique(z$family)), genus=as.character(unique(z$genus)), Reef_Nonreef=RnRF))
})

	# transform list to data.frame
W = do.call("rbind", Q)

	# write the file to a csv
write.csv(W, "RNR_stage_outfile_2018_Exp.csv", row.names=FALSE)
#Probably need to clean manually for duplicate periods and missing stages
#cleaned for Quaternary/Holocene/Neogene, primarily deletion of stageless rows with multiple periods-- adding unnecessary columns
# Known as "RNR_stage_outfile_2018_Exp_clean.csv", used later

#This element takes the original full data set and produces genus level Reef_Nonreef  percentage (or NaN's)


####Now to get randomized data to determine what would be the expectation of reefal occurrences given sampling
####Start with original data
R=read.csv(file="RNR_Master_May2018.csv", header=TRUE)
Z=read.csv("stage_numbered.csv", header=TRUE)
colnames(R)[10]="period"
colnames(R)[13]="stage" #12 or 13? Rogue version
RZ=merge(R,Z,c("period", "stage"))

#clean for empty stages
RZ = RZ[!(RZ$stage == ""),]
#Input for resamples
RZ_FamStub=unique(RZ[c("Num", "family", "genus")])


Sp=split(RZ, RZ$Num, drop=TRUE)
ii=1

#loop here

for (j in 1:200){
	RZ_Shuffle=lapply(seq_along(Sp), function(ii){
	oo=Sp[[ii]]
	tmpENV=sample(oo$Reef_Nonreef, replace=FALSE)
	data.frame(unique(oo$Num), Reef_Nonreef=tmpENV, oo$collection_no)	
		})

	y=do.call("rbind", RZ_Shuffle)
	colnames(y)=c("Num", "Shuf_Reef_Nonreef", "collection_no")
	RZy=merge(RZ,y, by=c("Num", "collection_no"))
	#Now that we have shuffled environments for all occurrences, let's hget ratios for individual genera/stage
	
	# subset to only neccessary fields
	Y = RZy[ , c("class", "order", "family", "genus", "period", "collection_no", "Num", "Shuf_Reef_Nonreef")]

	# Split the indeces of the dataframe (drop=TRUE means non-real combos of period and genus name are left out)
	P = split(seq(nrow(Y)) , list(Y$Num, Y$genus), drop=TRUE)


# loop over each genus and extract percent reefal per genus/stage
		Q = lapply(P, function(i){
			# subset indeces to actual data
		z = Y[i, ]
			# remove the occs with NAs for reef_nonreef
		q = z[!is.na(z$Shuf_Reef_Nonreef), ]

			# ratio of reef to non reef
		RnRF = nrow(q[ q$Shuf_Reef_Nonreef=="reef", ] ) / (nrow(q[ q$Shuf_Reef_Nonreef=="nonreef", ] ) + nrow(q[ q$Shuf_Reef_Nonreef=="reef", ] ))

		return(list(period=as.character(unique(z$period_updated)), Num=as.numeric(unique(z$Num)), class=as.character(unique(z$class)), order=as.character(unique(z$order)), family=as.character(unique(z$family)), genus=as.character(unique(z$genus)), Reef_Nonreef=RnRF))
		})

	# transform list to data.frame
	W = do.call("rbind", Q)
	Wsub=W[,c("Num", "family", "genus", "Reef_Nonreef")]
	#Should Reef_Nonreef be renamed to include loop number?
	RZ_FamStub=merge(RZ_FamStub, Wsub, by=c("Num", "family", "genus"))
	

}

write.csv(RZ_FamStub, file="RNR_simulation_200.csv")

#### Let's get the 95th percentile
RNR_simulation_200= read.csv(file="RNR_simulation_200.csv")
RNRquantileest=apply(RNR_simulation_200[5:204], 1, quantile, probs = c(0.05, 0.95),  na.rm = TRUE)
#Need to rotate transpose ot for stitching onto taxonomy
RNRquantileTranspose=t(RNRquantileest)

write.csv(RNRquantileTranspose, file="RNR200_Quantile_Est.csv")
#stitch onto list of taxonomy
RNR_sim_stub=RNR_simulation_200[1:4]
RNR_sim_stub=cbind(RNR_sim_stub, RNRquantileTranspose)

###"RNR_stage_outfile..." is the outfile containing Period, stage_no, genus-level taxonomy, their total occurences in that stage, the number of occurrences without environmental data, and a % covering reefal occurrences of the remainder

R=read.csv(file="RNR_stage_outfile_2018_Exp_clean.csv", header=TRUE)
#rename stage column to Num, for consistency
names(R)[2] = paste("Num")

RNR_sim_percentile=merge(RNR_sim_stub, R, c("Num", "family", "genus") )
write.csv(RNR_sim_percentile, file="RNR_sim_percentile_new.csv")

# subset of only the genera/stage combinations with Reef_Nonreef above 95th percentile
select=which(RNR_sim_percentile$Reef_Nonreef >= RNR_sim_percentile[6])
RNR_ReefalGeneraStage= RNR_sim_percentile[select,]
write.csv(RNR_ReefalGeneraStage, file="RNR_ReefalGenera_byStage_new.csv")

#and once more for minimum of 4 occurrences in a stage
select=(which((RNR_ReefalGeneraStage$total.occ-RNR_ReefalGeneraStage$na.occ) > 3))
RNR_ReefalGeneraStage_4occ= RNR_ReefalGeneraStage[select,]
write.csv(RNR_ReefalGeneraStage_4occ, file= "RNR_ReefalGenera_byStage_new_min4occ.csv")


#and once more for minimum of 8 occurrences in a stage
select=(which((RNR_ReefalGeneraStage$total.occ-RNR_ReefalGeneraStage$na.occ) > 7))
RNR_ReefalGeneraStage_8occ= RNR_ReefalGeneraStage[select,]
write.csv(RNR_ReefalGeneraStage_8occ, file= "RNR_ReefalGenera_byStage_new_min8occ.csv")


#### same as above, but for those definitely not reef

select=which(RNR_sim_percentile$Reef_Nonreef <= RNR_sim_percentile[5])
RNR_NonReefalGeneraStage= RNR_sim_percentile[select,]
write.csv(RNR_NonReefalGeneraStage, file="RNR_NonReefalGenera_byStage_new.csv")

#and once more for minimum of 4 occurrences in a stage
select=(which((RNR_NonReefalGeneraStage$total.occ-RNR_NonReefalGeneraStage$na.occ) > 3))
RNR_NonReefalGeneraStage_4occ= RNR_NonReefalGeneraStage[select,]
write.csv(RNR_NonReefalGeneraStage_4occ, file= "RNR_NonReefalGenera_byStage_min4occ.csv")

#and once more for minimum of 8 occurrences in a stage
select=(which((RNR_NonReefalGeneraStage$total.occ-RNR_NonReefalGeneraStage$na.occ) > 7))
RNR_NonReefalGeneraStage_8occ= RNR_NonReefalGeneraStage[select,]
write.csv(RNR_NonReefalGeneraStage_8occ, file= "RNR_NonReefalGenera_byStage_min8occ.csv")


#--------------
#### Lets identify extinction actoss stage boundaries

###First, let's identify all genus/stage combinations to make the search a bit easier

X = read.csv("RNR_stage_outfile_2018_Exp_clean.csv", header=TRUE)


	# subset to only neccessary fields
Extn.Test.Merge = X[ , c("period", "stage", "order", "family", "genus")]

	#Define reef builders
Reef.Builder.Orders=c("Actinocerida","Actinostromatida", "Ajacicyathida","Agelasida", "Amphidiscosa", "Amphiporida", "Archaeocyathida","Auloporida", "Axinellida", "Capsulocyathida","Chaetetida", "Clathrodictyida", "Clavulina", "Cystiphyllida", "Favositida", "Hadromerida", "Hexactinosa", "Hippuritoida", "Labechiida", "Lyssacinosa", "Monocyathida","Orthocladina", "Orthotetida", "Sarcinulida", "Scleractinia", "Sphaerocoeliida", "Spirosclerophorida", "Stauriida", "Stellispongiida", "Stromatoporida", "Syringostromatida", "Tabulaconida", "unranked(reef)", "Vaceletida")

	#Define reef fauna
Reef.Fauna=read.csv("RNR_ReefalGenera_byStage_new_min8occ.csv")
Reef.Builders=Reef.Fauna[which(Reef.Fauna$order %in% Reef.Builder.Orders),]
'%!in%' <- function(x,y)!('%in%'(x,y))
Reef.Dwellers=Reef.Fauna[which(Reef.Fauna$order %!in% Reef.Builder.Orders),]

	#Define crises
Crises=c(11, 21, 25, 31, 35, 37, 44, 47, 55, 56, 62)
Big5=c(11, 21, 37, 44, 67)
OtherCrises=c(25, 31, 35, 47, 55, 56, 62)

	#define printout stub
Reef.Transition = data.frame(matrix(ncol=16, nrow=0))

	
	#get looping

for(a in 1:73){
#Identify all occurences after stage
Extn.Tester=subset(Extn.Test.Merge, Extn.Test.Merge$stage > a)

Reef.Fauna.St=subset(Reef.Fauna, Num==a)
	Reef.Gen=unique(Reef.Fauna.St$genus)
Reef.Gen.No=length(Reef.Gen)
Reef.Fauna.Tester=subset(Reef.Fauna, Num > a)
	Ex.Count=0
	Surv.Count=0
	Reef.Count=0

	for(b in 1:Reef.Gen.No){
		Extinction = length(which(Extn.Tester$genus == as.character(Reef.Gen[b])))
		Persistence= length(which(Reef.Fauna.Tester$genus == as.character(Reef.Gen[b])))
				if (Extinction==0){Ex.Count=Ex.Count+1}
			else if (Extinction>0){Surv.Count=Surv.Count+1}
				if(Persistence>0){Reef.Count=Reef.Count+1}			
			}
			if (Reef.Gen.No==0){Ex.Count=0}
			if (Reef.Gen.No==0){Surv.Count=0}
			if (Reef.Gen.No==0){Reef.Count=0}
	Ex.pc=Ex.Count/(Ex.Count+Surv.Count)
	Surv.pc=Surv.Count/(Ex.Count+Surv.Count)
	Reef.pc=Reef.Count/Surv.Count
	
	#Calculate error
	z<-1
Ex.plo<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.pup<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.se.p<-0.5*(Ex.pup-Ex.plo)

Surv.plo<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.pup<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.se.p<-0.5*(Surv.pup-Surv.plo)

Reef.plo<-(2*(Surv.Count)*Reef.pc+z^2-sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.pup<-(2*(Surv.Count)*Reef.pc+z^2+sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.se.p<-0.5*(Reef.pup-Reef.plo)	
	
	Reef.Transition.List= c(a,Ex.Count, Surv.Count, Reef.Count, Ex.pc, Surv.pc, Reef.pc, Ex.se.p, Surv.se.p, Reef.se.p, Ex.plo, Surv.plo, Reef.plo, Ex.pup, Surv.pup, Reef.pup)
	Reef.Transition=rbind(Reef.Transition, Reef.Transition.List)
	
}
x=c("Stage", "Extinction", "Survival", "Reef_Persist", "Extinction_%", "Survival_%", "Reef_Persist_%", "Ex_Standard_Error", "Surv_Standard_Error", "Reef_Pers_Standard Error", "Ex_error_lo", "Surv_error_lo", "Reef_error_lo", "Ex_error_up", "Surv_error_up", "Reef_error_up")
colnames(Reef.Transition)=x

write.csv(Reef.Transition, file="Reef.Transition.2018.8occ.csv")



#Version for nonreef

NonReef.Fauna=read.csv("RNR_NonReefalGenera_byStage_min8occ.csv")


	#Define crises
Crises=c(11, 21, 25, 31, 35, 37, 44, 47, 55, 56, 62)
Big5=c(11, 21, 37, 44, 67)
OtherCrises=c(25, 31, 35, 47, 55, 56, 62)

	#define printout stub
NonReef.Transition = data.frame(matrix(ncol=10, nrow=0))

	
	#get looping

for(a in 1:73){
#Identify all occurences after stage
Extn.Tester=subset(Extn.Test.Merge, Extn.Test.Merge$stage > a)

NonReef.Fauna.St=subset(NonReef.Fauna, Num==a)
	NonReef.Gen=unique(NonReef.Fauna.St$genus)
NonReef.Gen.No=length(NonReef.Gen)
NonReef.Fauna.Tester=subset(NonReef.Fauna, Num > a)
	Ex.Count=0
	Surv.Count=0
	Reef.Count=0

	for(b in 1:NonReef.Gen.No){
		Extinction = length(which(Extn.Tester$genus == as.character(NonReef.Gen[b])))
		Persistence= length(which(NonReef.Fauna.Tester$genus == as.character(Reef.Gen[b])))
				if (Extinction==0){Ex.Count=Ex.Count+1}
			else if (Extinction>0){Surv.Count=Surv.Count+1}
				if(Persistence>0){Reef.Count=Reef.Count+1}			
			}
			if (NonReef.Gen.No==0){Ex.Count=0}
			if (NonReef.Gen.No==0){Surv.Count=0}
			if (NonReef.Gen.No==0){Reef.Count=0}
	Ex.pc=Ex.Count/(Ex.Count+Surv.Count)
	Surv.pc=Surv.Count/(Ex.Count+Surv.Count)
	Reef.pc=Reef.Count/Surv.Count
	
	#Calculate error
	z<-1
Ex.plo<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.pup<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.se.p<-0.5*(Ex.pup-Ex.plo)

Surv.plo<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.pup<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.se.p<-0.5*(Surv.pup-Surv.plo)

Reef.plo<-(2*(Surv.Count)*Reef.pc+z^2-sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.pup<-(2*(Surv.Count)*Reef.pc+z^2+sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.se.p<-0.5*(Reef.pup-Reef.plo)	
	
	NonReef.Transition.List= c(a,Ex.Count, Surv.Count, Reef.Count, Ex.pc, Surv.pc, Reef.pc, Ex.se.p, Surv.se.p, Reef.se.p, Ex.plo, Surv.plo, Reef.plo, Ex.pup, Surv.pup, Reef.pup)
	NonReef.Transition=rbind(NonReef.Transition, NonReef.Transition.List)
	
}
x=c("Stage", "Extinction", "Survival", "Reef_Persist", "Extinction_%", "Survival_%", "Reef_Persist_%", "Ex_Standard_Error", "Surv_Standard_Error", "Reef_Pers_Standard Error", "Ex_error_lo", "Surv_error_lo", "Reef_error_lo", "Ex_error_up", "Surv_error_up", "Reef_error_up")
colnames(NonReef.Transition)=x

write.csv(NonReef.Transition, file="NonReef.Transition.2018.8occ.csv")




	#Version for Reef Builders
Reef.Builder.Transition = data.frame(matrix(ncol=10, nrow=0))

	for(a in 1:73){
#Identify all occurences after stage
Extn.Tester=subset(Extn.Test.Merge, Extn.Test.Merge$stage > a)

Reef.Builder.St=subset(Reef.Builders, Num==a)
	Reef.Gen=unique(Reef.Builder.St$genus)
Reef.Builder.Gen.No=length(Reef.Gen)
Reef.Builder.Tester=subset(Reef.Builders, Num > a)
	Ex.Count=0
	Surv.Count=0
	Reef.Count=0

	for(b in 1:Reef.Builder.Gen.No){
		Extinction = length(which(Extn.Tester$genus == as.character(Reef.Gen[b])))
		Persistence= length(which(Reef.Builder.Tester$genus == as.character(Reef.Gen[b])))
				if (Extinction==0){Ex.Count=Ex.Count+1}
			else if (Extinction>0){Surv.Count=Surv.Count+1}
				if(Persistence>0){Reef.Count=Reef.Count+1}			
			}
			if (Reef.Builder.Gen.No==0){Ex.Count=0}
			if (Reef.Builder.Gen.No==0){Surv.Count=0}
			if (Reef.Builder.Gen.No==0){Reef.Count=0}			
	Ex.pc=Ex.Count/(Ex.Count+Surv.Count)
	Surv.pc=Surv.Count/(Ex.Count+Surv.Count)
	Reef.pc=Reef.Count/Surv.Count
	
		#Calculate error
	z<-1
Ex.plo<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.pup<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.se.p<-0.5*(Ex.pup-Ex.plo)

Surv.plo<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.pup<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.se.p<-0.5*(Surv.pup-Surv.plo)

Reef.plo<-(2*(Surv.Count)*Reef.pc+z^2-sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.pup<-(2*(Surv.Count)*Reef.pc+z^2+sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.se.p<-0.5*(Reef.pup-Reef.plo)
	
	Reef.Builder.Transition.List= c(a,Ex.Count, Surv.Count, Reef.Count, Ex.pc, Surv.pc, Reef.pc, Ex.se.p, Surv.se.p, Reef.se.p, Ex.plo, Surv.plo, Reef.plo, Ex.pup, Surv.pup, Reef.pup)
	Reef.Builder.Transition=rbind(Reef.Builder.Transition, Reef.Builder.Transition.List)
	
}
x=c("Stage", "Extinction", "Survival", "Reef_Persist", "Extinction_%", "Survival_%", "Reef_Persist_%", "Ex_Standard_Error", "Surv_Standard_Error", "Reef_Pers_Standard Error", "Ex_error_lo", "Surv_error_lo", "Reef_error_lo", "Ex_error_up", "Surv_error_up", "Reef_error_up")
colnames(Reef.Builder.Transition)=x

write.csv(Reef.Builder.Transition, file="Reef.Builder.Transition.2018.8occ.csv")
	
	#Version for Reef Dwellers
	
	Reef.Dweller.Transition = data.frame(matrix(ncol=10, nrow=0))

	for(a in 1:73){
#Identify all occurences after stage
Extn.Tester=subset(Extn.Test.Merge, Extn.Test.Merge$stage > a)

Reef.Dweller.St=subset(Reef.Dwellers, Num==a)
	Reef.Gen=unique(Reef.Dweller.St$genus)
Reef.Dweller.Gen.No=length(Reef.Gen)
Reef.Dweller.Tester=subset(Reef.Dwellers, Num > a)
	Ex.Count=0
	Surv.Count=0
	Reef.Count=0

	for(b in 1:Reef.Dweller.Gen.No){
		Extinction = length(which(Extn.Tester$genus == as.character(Reef.Gen[b])))
		Persistence= length(which(Reef.Dweller.Tester$genus == as.character(Reef.Gen[b])))
				if (Extinction==0){Ex.Count=Ex.Count+1}
			else if (Extinction>0){Surv.Count=Surv.Count+1}
				if(Persistence>0){Reef.Count=Reef.Count+1}			
			}
			if (Reef.Dweller.Gen.No==0){Ex.Count=0}
			if (Reef.Dweller.Gen.No==0){Surv.Count=0}
			if (Reef.Dweller.Gen.No==0){Reef.Count=0}			
	Ex.pc=Ex.Count/(Ex.Count+Surv.Count)
	Surv.pc=Surv.Count/(Ex.Count+Surv.Count)
	Reef.pc=Reef.Count/Surv.Count
	
		#Calculate error
	z<-1
Ex.plo<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.pup<-(2*(Ex.Count+Surv.Count)*Ex.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Ex.pc*(Ex.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Ex.se.p<-0.5*(Ex.pup-Ex.plo)

Surv.plo<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2-sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.pup<-(2*(Ex.Count+Surv.Count)*Surv.pc+z^2+sqrt(z^2*(z^2-4*(Ex.Count+Surv.Count)*Surv.pc*(Surv.pc-1))))/(2*((Ex.Count+Surv.Count)+z^2))

Surv.se.p<-0.5*(Surv.pup-Surv.plo)

Reef.plo<-(2*(Surv.Count)*Reef.pc+z^2-sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.pup<-(2*(Surv.Count)*Reef.pc+z^2+sqrt(z^2*(z^2-4*(Surv.Count)*Reef.pc*(Reef.pc-1))))/(2*((Surv.Count)+z^2))

Reef.se.p<-0.5*(Reef.pup-Reef.plo)

	
	Reef.Dweller.Transition.List= c(a,Ex.Count, Surv.Count, Reef.Count, Ex.pc, Surv.pc, Reef.pc, Ex.se.p, Surv.se.p, Reef.se.p, Ex.plo, Surv.plo, Reef.plo, Ex.pup, Surv.pup, Reef.pup)
	Reef.Dweller.Transition=rbind(Reef.Dweller.Transition, Reef.Dweller.Transition.List)
	
}
x=c("Stage", "Extinction", "Survival", "Reef_Persist", "Extinction_%", "Survival_%", "Reef_Persist_%", "Ex_Standard_Error", "Surv_Standard_Error", "Reef_Pers_Standard Error", "Ex_error_lo", "Surv_error_lo", "Reef_error_lo", "Ex_error_up", "Surv_error_up", "Reef_error_up")
colnames(Reef.Dweller.Transition)=x

write.csv(Reef.Dweller.Transition, file="Reef.Dweller.Transition.2018.8occ.csv")



#Stenotopy vs Eurytopy

#stage by stage-- get the average %reefal for nonreef and the avg %reefal for reef taxa. Plot against each other

ReefalGenera=read.csv(file="RNR_ReefalGenera_byStage_new_min4occ.csv")
NonReefalGenera=read.csv(file="RNR_NonReefalGenera_byStage_min4occ.csv")
StenTest=data.frame(matrix(ncol=5, nrow=0))


for (i in 1:73){
	
	StenTestReef=subset(ReefalGenera, Num==i)
	StenTestNonReef=subset(NonReefalGenera, Num==i)
	StenTestReef.Mean=mean(as.numeric(StenTestReef$Reef_Nonreef))
	StenTestNonReef.Mean=mean(as.numeric(StenTestNonReef$Reef_Nonreef))
	
	p=StenTestReef.Mean
	q=StenTestNonReef.Mean
	n=length(StenTestReef$Reef_NonReef)
	m=length(StenTestNonReef$Reef_NonReef)
	z=1
	
	plo=(2*n*p+z^2 -sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
	pup=(2*n*p+z^2 +sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
	plo.nr=(2*m*q+z^2 -sqrt(z^2*(z^2-4*m*q*(q-1))))/(2*(m+z^2))
	pup.nr=(2*m*q+z^2 +sqrt(z^2*(z^2-4*m*q*(q-1))))/(2*(m+z^2))
	se.p=0.5*(pup-plo)	
	se.p.nr=0.5*(pup.nr-plo.nr)
	
	
		
StenTester=c(i, StenTestReef.Mean, StenTestNonReef.Mean, se.p, se.p.nr)
StenTest=rbind(StenTester, StenTest)
print(i)	
}

plot(StenTest[,2], StenTest[,3], xlim=c(0,1), ylim=c(0,1), ylab="Nonreef occupancy")

#############################################
#############################################
#### PLOTS PLOTS PLOTS PLOTS PLOTS PLOTS#####


Reef.Transition=read.csv(file="Reef.Transition.2018.csv", header=TRUE)
Reef.Builder.Transition=read.csv(file="Reef.Builder.Transition.2018.csv", header=TRUE)
Reef.Dweller.Transition=read.csv(file="Reef.Dweller.Transition.2018.csv", header=TRUE)

#OR

Reef.Transition=read.csv(file="Reef.Transition.2018.8occ.csv", header=TRUE)
Reef.Builder.Transition=read.csv(file="Reef.Builder.Transition.2018.8occ.csv", header=TRUE)
Reef.Dweller.Transition=read.csv(file="Reef.Dweller.Transition.2018.8occ.csv", header=TRUE)

#OR
Reef.Transition=read.csv(file="Reef.Transition.2018.4occ.csv", header=TRUE)
Reef.Builder.Transition=read.csv(file="Reef.Builder.Transition.2018.4occ.csv", header=TRUE)
Reef.Dweller.Transition=read.csv(file="Reef.Dweller.Transition.2018.4occ.csv", header=TRUE)

#OR

Reef.Transition=read.csv(file="Reef.Transition.2018.4occ_2stages.csv", header=TRUE)
Reef.Builder.Transition=read.csv(file="Reef.Builder.Transition.2018.4occ_2stage.csv", header=TRUE)
Reef.Dweller.Transition=read.csv(file="Reef.Dweller.Transition.2018.4occ_2stage.csv", header=TRUE)

#OR

Reef.Transition=read.csv(file="Reef.Transition.2018.8occ_2stages.csv", header=TRUE)
Reef.Builder.Transition=read.csv(file="Reef.Builder.Transition.2018.8occ_2stage.csv", header=TRUE)
Reef.Dweller.Transition=read.csv(file="Reef.Dweller.Transition.2018.8occ_2stage.csv", header=TRUE)


# ExtinctionFor all reef fauna
	Reef.Transition.Crises=subset(Reef.Transition, Reef.Transition[,2] %in% Crises)
	Reef.Transition.Big5=subset(Reef.Transition, Reef.Transition[,2] %in% Big5)

	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Transition[,2], y=Reef.Transition[,6], type='p', xlab="Stage Transition", ylab="Proportional Extinction", main="Extinction of Reefal Genera", ylim=c(0,1))
	points(x=Crises, y=Reef.Transition.Crises[,6], pch=16, col="red")
	points(x=Big5, y=Reef.Transition.Big5[,6], pch=5, col="blue")
	arrows(Reef.Transition[,2], Reef.Transition[,12], Reef.Transition[,2], Reef.Transition[,15], length=0.05, angle=90, code=3)
	legend(62,1, c("Stage Boundary", "Reef Crisis", "Big 5 Extn"),pch=c(1,16, 5), col= c("black", "red", "blue"))	
		

#Extinction For Builders
	Reef.Builder.Transition.Crises=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Crises)
	Reef.Builder.Transition.Big5=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Big5)
	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Builder.Transition[,2], y=Reef.Builder.Transition[,6], type='p', xlab="Stage Boundary", ylab="Proportional Extinction", main="Extinction of Reef-builder Genera")
	points(x=Crises, y=Reef.Builder.Transition.Crises[,6], pch=16, col="red")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,6], pch=5, col="blue")

# Extinction For Dwellers
	Reef.Dweller.Transition.Crises=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Crises)
	Reef.Dweller.Transition.Big5=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Big5)
	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Dweller.Transition[,2], y=Reef.Dweller.Transition[,6], type='p', xlab="Stage Boundary", ylab="Proportional Extinction", main="Extinction of Reef-dweller Genera")
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,6], pch=16, col="red")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,6], pch=5, col="blue")
	
	#Double plot
		Reef.Builder.Transition.Crises=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Crises)
	Reef.Builder.Transition.Big5=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Big5)
		Reef.Dweller.Transition.Crises=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Crises)
	Reef.Dweller.Transition.Big5=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Big5)	
	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Builder.Transition[,2], y=Reef.Builder.Transition[,6], type='p', xlab="Stage Boundary", ylab="Proportional Extinction", main="Extinction of Reef-Builders and Reef-Dwellers", col="blue")
	points(x=Crises, y=Reef.Builder.Transition.Crises[,6], pch=16, col="blue")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,6], pch=5, col="blue")	
	
	points(x=Reef.Dweller.Transition[,2], y=Reef.Dweller.Transition[,6], col="darkorange")
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,6], pch=16, col="darkorange")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,6], pch=5, col="darkorange")	
	
	arrows(Reef.Builder.Transition[,2], Reef.Builder.Transition[,12], Reef.Builder.Transition[,2], Reef.Builder.Transition[,15], length=0.05, angle=90, code=3, col="blue")	
	arrows(Reef.Dweller.Transition[,2], Reef.Dweller.Transition[,12], Reef.Dweller.Transition[,2], Reef.Dweller.Transition[,15], length=0.05, angle=90, code=3, col="darkorange")
	legend(8,1, c("Stage Boundary","Reef Crisis", "Big 5 Extn", "Reef-Builders", "Reef-Dwellers"), pch=c(1, 16, 5, 15, 15), col= c("black", "black", "black", "blue", "darkorange"), cex=.9)	
	
#Double plot with just crises
	
		Reef.Builder.Transition.Crises=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Crises)
	Reef.Builder.Transition.Big5=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Big5)
		Reef.Dweller.Transition.Crises=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Crises)
	Reef.Dweller.Transition.Big5=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Big5)	
	dev.new(width=11, height=5, unit="in")
	plot(x=Crises, y=Reef.Builder.Transition.Crises[,6], pch=16, col="red",xlab="Reef Crises", ylab="Proportional Extinction",xlim=c(0,73), ylim=c(0,1), main="Extinction of Reef Fauna over Reef Crises")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,6], pch=5, col="blue")	
	
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,6], pch=16, col="green")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,6], pch=5, col="darkorange")	
	
		arrows(Crises, Reef.Builder.Transition.Crises[,12], Crises, Reef.Builder.Transition.Crises[,15], length=0.05, angle=90, code=3, col="grey")	
	arrows(Big5, Reef.Dweller.Transition.Big5[,12], Big5, Reef.Dweller.Transition.Big5[,15], length=0.05, angle=90, code=3)
	arrows(Crises, Reef.Dweller.Transition.Crises[,12], Crises, Reef.Dweller.Transition.Crises[,15], length=0.05, angle=90, code=3)	
	arrows(Big5, Reef.Builder.Transition.Big5[,12], Big5, Reef.Builder.Transition.Big5[,15], length=0.05, angle=90, code=3, col="grey")	
		
		
	legend(0,1, c("Reef Crisis", "Big 5 Extn", "Reef-builder", "Reef-dweller"),pch=c(16, 5, 5, 5), col= c("black", "black", "blue", "darkorange"))
	
	
	
	
	
	
	
#Reef Persistence of All Reef Fauna
	Reef.Transition.Crises=subset(Reef.Transition, Reef.Transition[,2] %in% Crises)
	Reef.Transition.Big5=subset(Reef.Transition, Reef.Transition[,2] %in% Big5)
	Reef.Transition.OtherCrises=subset(Reef.Transition, Reef.Transition[,2] %in% OtherCrises)

	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Transition[,2], y=Reef.Transition[,8], type='p', xlab="Stage Boundary", ylab="Proportion Sustained Reef Preference", main="Sustained Reef Preference of Surviving Reefal Genera", ylim=c(0,1))
	points(x=Crises, y=Reef.Transition.Crises[,8], pch=16, col="blue")
	points(x=Big5, y=Reef.Transition.Big5[,8], pch=5, col="blue")
	
	arrows(Reef.Transition[,2], Reef.Transition[,14], Reef.Transition[,2], Reef.Transition[,17], length=0.05, angle=90, code=3)
	legend(0,1, c("Stage boundary","Reef Crisis", "Big 5 Extn"),pch=c(1,16, 5), col= c("black", "blue", "blue"))	

#Reef Persistence of Builders
	Reef.Builder.Transition.Crises=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Crises)
	Reef.Builder.Transition.Big5=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Big5)
	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Builder.Transition[,2], y=Reef.Builder.Transition[,8], type='p', xlab="Stage Boundary", ylab="Proportion Sustained Reef Preference", main="Sustained Reef Preference of Surviving Reef-builder Genera")
	points(x=Crises, y=Reef.Builder.Transition.Crises[,8], pch=16, col="red")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,8], pch=5, col="blue")

# Reef Persistence of Dwellers
	Reef.Dweller.Transition.Crises=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Crises)
	Reef.Dweller.Transition.Big5=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Big5)
	dev.new(width=11, height=5, unit="in")
	plot(x=Reef.Dweller.Transition[,2], y=Reef.Dweller.Transition[,8], type='p', xlab="Stage Boundary", ylab="Proportion Sustained Reef Preference", main="Sustained Reef Preference of Surviving Reef-dweller Genera")
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,8], pch=16, col="red")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,8], pch=5, col="blue")	
	
	#Reef Persistence of Both
	Reef.Builder.Transition.Crises=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Crises)
	Reef.Builder.Transition.Big5=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Big5)
	Reef.Dweller.Transition.Crises=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Crises)
	Reef.Dweller.Transition.Big5=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Big5)	
	dev.new(width=11, height=5, unit="in")
	plot(x=Crises, y=Reef.Builder.Transition.Crises[,8], type='p', pch=16, col="blue", xlim=c(0,73), ylim=c(0,1), xlab="Stage Boundary", ylab="Proportion Sustained Reef Preference", main="Sustained Reef Preference for Reef-builders and Reef-Dwellers")
	
	points(x=Reef.Builder.Transition[,2], y=Reef.Builder.Transition[,8], col="blue")
	points(x=Reef.Dweller.Transition[,2], y=Reef.Dweller.Transition[,8], col="darkorange")
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,8], pch=16, col="darkorange")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,8], pch=5, col="blue")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,8], pch=5, col="darkorange")
	
	arrows(Reef.Builder.Transition[,2], Reef.Builder.Transition[,14], Reef.Builder.Transition[,2], Reef.Builder.Transition[,17], length=0.05, angle=90, code=3, col="blue")	
	arrows(Reef.Dweller.Transition[,2], Reef.Dweller.Transition[,14], Reef.Dweller.Transition[,2], Reef.Dweller.Transition[,17], length=0.05, angle=90, code=3, col="darkorange")	
	legend(-1,1, c("Stage boundary","Reef Crisis", "Big 5 Extn", "Reef-Builders", "Reef-Dwellers"),pch=c(1,16, 5, 15, 15), col= c("black","black", "black", "blue", "darkorange"), cex=.75)		
	

	
	
	
#Reef Persistence of Both during crises
	Reef.Builder.Transition.Crises=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Crises)
	Reef.Builder.Transition.Big5=subset(Reef.Builder.Transition, Reef.Builder.Transition[,2] %in% Big5)
	Reef.Dweller.Transition.Crises=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Crises)
	Reef.Dweller.Transition.Big5=subset(Reef.Dweller.Transition, Reef.Dweller.Transition[,2] %in% Big5)	
	dev.new(width=11, height=5, unit="in")
	plot(x=Crises, y=Reef.Builder.Transition.Crises[,8], type='p', pch=16, col="blue", xlim=c(0,73), ylim=c(0,1), xlab="Stage Boundary", ylab="Proportion Sustained Reef Preference", main="Sustained Reef Preference over Reef Crises")
	
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,8], pch=16, col="darkorange")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,8], pch=5, col="blue")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,8], pch=5, col="darkorange")
	arrows(Crises, Reef.Builder.Transition.Crises[,14], Crises, Reef.Builder.Transition.Crises[,17], length=0.05, angle=90, code=3, col="blue")	
	arrows(Big5, Reef.Dweller.Transition.Big5[,14], Big5, Reef.Dweller.Transition.Big5[,17], length=0.05, angle=90, code=3, col="darkorange")
	arrows(Crises, Reef.Dweller.Transition.Crises[,14], Crises, Reef.Dweller.Transition.Crises[,17], length=0.05, angle=90, code=3, col="darkorange")	
	arrows(Big5, Reef.Builder.Transition.Big5[,14], Big5, Reef.Builder.Transition.Big5[,17], length=0.05, angle=90, code=3, col="blue")	
	legend(-1,1, c("Reef Crisis", "Big 5 Extn", "Reef-Builders", "Reef-Dwellers"),pch=c(16, 5, 15, 15), col= c("black", "black", "blue", "darkorange"), cex=.9)		
	
	
RNR_ReefalGeneraStage=read.csv(file="RNR_ReefalGenera_byStage_new.csv")

#and once more for minimum of 8 occurrences in a stage, the apply to the above
select=(which(RNR_ReefalGeneraStage$total.occ > 7))
RNR_ReefalGeneraStage_8occ= RNR_ReefalGeneraStage[select,]
write.csv(RNR_ReefalGeneraStage_8occ, file= "RNR_ReefalGenera_byStage_new_min8occ.csv")

	
	
####Comparison reefal and SIGNIFICANTLY nonreefal, all stages

#NonReef.Transition=read.csv(file="NonReef.Transition.2018.4occ.csv", header=TRUE)
NonReef.Transition=read.csv(file="NonReef.Transition.2018.8occ.csv", header=TRUE)

# Extinction For all reef and nonreef fauna
	NonReef.Transition.Crises=subset(NonReef.Transition, NonReef.Transition[,2] %in% Crises)
	NonReef.Transition.Big5=subset(NonReef.Transition, NonReef.Transition[,2] %in% Big5)
	
	dev.new(width=11, height=5, unit="in")
	plot(x=NonReef.Transition[,2], y=NonReef.Transition[,6], type='p',col="darkorange", xlab="Stage Transition", ylab="Proportional Extinction", main="Extinction of Reefal and Non-Reefal Genera", ylim=c(0,1))
	points(x=Reef.Transition[,2], y=Reef.Transition[,6], type='p', col="blue")
	
	points(x=Crises, y=NonReef.Transition.Crises[,6], pch=16, col="darkorange")
	points(x=Big5, y=NonReef.Transition.Big5[,6], pch=5, col="darkorange")
	
	points(x=Crises, y=Reef.Transition.Crises[,6], pch=16, col="blue")
	points(x=Big5, y=Reef.Transition.Big5[,6], pch=5, col="blue")
	
	arrows(Reef.Transition[,2], Reef.Transition[,12], Reef.Transition[,2], Reef.Transition[,15], length=0.05, angle=90, code=3, col="blue")
	arrows(NonReef.Transition[,2], NonReef.Transition[,12], NonReef.Transition[,2], NonReef.Transition[,15], length=0.05, angle=90, code=3, col="darkorange")
	legend(62,1, c("Stage Boundary","Reef Crisis", "Big 5 Extn", "Reefal", "Non-Reefal"), pch=c(1, 16, 5, 15, 15), col= c("black", "black", "black", "blue", "darkorange"), cex=.9)
	
#or just crises and big5
	
	dev.new(width=11, height=5, unit="in")
	plot(x=Crises, y=NonReef.Transition.Crises[,6], pch=16,col="darkorange", xlab="Stage Transition", ylab="Proportional Extinction",xlim=c(0,73), ylim=c(0,1), main="Extinction of Reefal and Significantly Non-Reefal Genera over Reef Crises")
	points(x=Crises, y=Reef.Transition.Crises[,6], pch=16, col="blue")
	points(x=Big5, y=NonReef.Transition.Big5[,6], pch=5, col="darkorange")
	points(x=Big5, y=Reef.Transition.Big5[,6], pch=5, col="blue")
	
	arrows(Crises, Reef.Transition.Crises[,12], Crises, Reef.Transition.Crises[,15], length=0.05, angle=90, code=3, col="blue")
	arrows(Crises, NonReef.Transition.Crises[,12], Crises, NonReef.Transition.Crises[,15], length=0.05, angle=90, code=3, col="darkorange")
	arrows(Big5, Reef.Transition.Big5[,12], Big5, Reef.Transition.Big5[,15], length=0.05, angle=90, code=3, col="blue")
	arrows(Big5, NonReef.Transition.Big5[,12], Big5, NonReef.Transition.Big5[,15], length=0.05, angle=90, code=3, col="darkorange")	
		
	legend(0,1, c("Reef Crisis", "Big 5 Extn", "Reefal", "Non-Reefal"), pch=c(16, 5, 5, 5), col= c("black", "black", "blue", "grey"))

		
		
		
#Comparison reef builder and dweller and nonreefal
	
	dev.new(width=11, height=5, unit="in")
	plot(x=Not.Reef.Transition[,2], y=Not.Reef.Transition[,6], type='p',col="grey", xlab="Stage Transition", ylab="Proportional Extinction", main="Extinction of Reefal and Non-Reefal Genera")
	points(x=Reef.Builder.Transition[,2], y=Reef.Builder.Transition[,6], type='p', col="blue")
	points(x=Reef.Dweller.Transition[,2], y=Reef.Dweller.Transition[,6], type='p', col="darkorange")
		
	points(x=Crises, y=Not.Reef.Transition.Crises[,6], pch=16, col="green")
	points(x=Big5, y=Not.Reef.Transition.Big5[,6], pch=5, col="green")
	
	points(x=Crises, y=Reef.Transition.Crises[,6], pch=16, col="red")
	points(x=Big5, y=Reef.Transition.Big5[,6], pch=5, col="red")
	
	
	
#or just crises and big5
	
	dev.new(width=11, height=5, unit="in")
	plot(x=Crises, y=Not.Reef.Transition.Crises[,6], pch=16,col="grey", xlab="Stage Transition", ylab="Extinction, per capita",xlim=c(0,73), ylim=c(0,1), main="Extinction of Reefal and Non-Reefal Genera over Reef Crises")
	
	points(x=Big5, y=Not.Reef.Transition.Big5[,6], pch=5, col="black")
	points(x=Crises, y=Reef.Builder.Transition.Crises[,6], pch=16, col="red")
	points(x=Big5, y=Reef.Builder.Transition.Big5[,6], pch=5, col="blue")	
	
	points(x=Crises, y=Reef.Dweller.Transition.Crises[,6], pch=16, col="green")
	points(x=Big5, y=Reef.Dweller.Transition.Big5[,6], pch=5, col="darkorange")	
	
	legend(0,1, c("Reef Crisis", "Big 5 Extn", "Reef-builder", "Reef-dweller", "All under reefal cutoff"),pch=c(16, 5, 5, 5, 5), col= c("black", "black", "blue", "darkorange", "black"))
	





##############################
##############################
###Determining geographic restriction or narrow habitat range("stenotopy")
###Work  continued off example code supplied by Professor Michael Foote. Thanks, once again, Michael.

##First, we need to determine reefal affinity, as the former method cannot be applied
##We do this purely statistically using the known distributions, rather than manual shuffling

cutoff<-4 #minimum number of occurrences to determine affinity

infile<-"RNR_Master_May2018_StageNo_Builders_LatLong.csv"
X<-read.csv(infile)
R<-as.character(X$Reef_Nonreef)
plat<-X$paleolat
plon<-X$paleolng

#ii<-which(!is.na(R)&!is.na(plat)&!is.na(plon)) #remove occurrences not known as reef or non-reef and those without coordinates
#X<-X[ii,]

genus<-as.character(X$genus)
R<-as.character(X$Reef_Nonreef)
plat<-X$paleolat
plon<-X$paleolng

bin<-X$Num
nbin<-max(bin)

ugen<-sort(unique(genus))
ngen<-length(ugen)

Occ<-table(genus,bin) #implicitly sorts genera alphabetically and bins numerically
C<-col(Occ)
SBN<-ifelse(Occ>0,C,NA) #sampled in bin number
FO<-apply(SBN,1,min,na.rm=TRUE) #first appearance bin
LO<-apply(SBN,1,max,na.rm=TRUE) #last appearance bin

#genus-by-bin affinities: 0=nonreef, 1=reef, NA=no affinity or insufficient data

AFF<-matrix(NA,ngen,nbin)
NR<-NNR<-matrix(NA,ngen,nbin)

for (i in 1:ngen)
{
 cat(i,ngen,"\n")
 for (j in FO[i]:LO[i])
 {
  ii<-which(genus==ugen[i]&bin==j)
  mm<-which(bin==j)
  #mreef<-sum(R[mm]=="reef")
  #mnonreef<-sum(R[mm]=="nonreef")
  mreef<-length(which((R[mm]=="reef")))
  mnonreef<-length(which((R[mm]=="nonreef")))
  n=length(ii)
  p<-(mreef/length(mm))
  q<-(mnonreef/length(mm))
  if (length(ii)>0)
  {
   #nr<-sum(R[ii]=="reef")
   #nnr<-sum(R[ii]=="nonreef")
   nr<-length(which((R[ii]=="reef")))
   nnr<-length(which((R[ii]=="nonreef")))
   NR[i,j]<-nr
   NNR[i,j]<-nnr
   
   Reef.affinity=1-pbinom(nr-1,n,p)
   Nonreef.affinity=1-pbinom(nnr-1,n,q)
   if (is.na(nr) & is.na(nnr)){AFF[i,j]<-NA} else if (nr+nnr>=cutoff)
   	{
    if (is.na(Reef.affinity)) AFF[i,j]<-NA else if (Reef.affinity<.05) AFF[i,j]<-1 else if (Nonreef.affinity<.05) AFF[i,j]<-0 else AFF[i,j]<-NA
    }
   }
  }
 }


write.csv(NR, file="Reef_Affinity_FooteHybrid_4occ_fulldataset.csv")
write.csv(NNR, file="Nonreef_Affinity_FooteHybrid_4occ_fulldataset.csv")
write.csv(AFF, file="Affinity_FooteHybrid_4occ_fulldataset.csv")



####Variation using Foote's suggested 50:50 reef determination

#for (i in 1:ngen)
#{
 #cat(i,ngen,"\n")
 #for (j in FO[i]:LO[i])
 #{
  #ii<-which(genus==ugen[i]&bin==j)
  #if (length(ii)>0)
  #{
   #nr<-sum(R[ii]=="reef")
   #nnr<-sum(R[ii]=="nonreef")
   #NR[i,j]<-nr
   #NNR[i,j]<-nnr
   #if (nr+nnr>=cutoff)
  #{
    #if (nr>nnr) AFF[i,j]<-1 else AFF[i,j]<-0
   #}
  #}
 #}
#}



#proportion of occurrences in reef vs. nonreef environments
PR<-NR/(NR+NNR)
PNR<-NNR/(NR+NNR)

#number and proportion of occurrences, by stage
NR.stage<-apply(NR,2,sum,na.rm=TRUE)
NNR.stage<-apply(NNR,2,sum,na.rm=TRUE)

NR.AFF<-NNR.AFF<-rep(0,nbin)
for (j in 1:nbin)
{
 NR.AFF[j]<-sum(AFF[,j]==1,na.rm=TRUE)
 NNR.AFF[j]<-sum(AFF[,j]==0,na.rm=TRUE)
}

#Geographic range (restricted vs. widespread)
GR<-matrix(NA,ngen,nbin)
min.deg<-5 #minimum range in degrees to consider genus "widespread" (both lat and lng)
for (i in 1:ngen)
{
 cat("GR",i,"\n")
 for (j in FO[i]:LO[i])
 {
  ii<-which(genus==ugen[i]&bin==j)
  if (length(ii)>0)
  {
   lat.rng<-range(plat[ii])
   lng.rng<-range(plon[ii])
   if (!is.na(lat.rng))
   {
   
   if (diff(lat.rng)<min.deg|diff(lng.rng)<min.deg) GR[i,j]<-0 else GR[i,j]<-1
   }
  }
 }
}

#write.csv(GR, file="GeogRange_FooteHybrid_4occ_fulldataset.csv")
#write.csv(NR.AFF, file="NR_AFF_4occ_fulldataset.csv")
#write.csv(NNR.AFF, file="NNR_AFF_4occ_fulldataset.csv")

#just in case
NR=read.csv(file="Reef_Affinity_FooteHybrid_4occ_fulldataset.csv")
NNR=read.csv(file="Nonreef_Affinity_FooteHybrid_4occ_fulldataset.csv")
AFF=read.csv(file="Affinity_FooteHybrid_4occ_fulldataset.csv")
GR=read.csv(file="GeogRange_FooteHybrid_4occ_fulldataset.csv")
NR.AFF=read.csv(file="NR_AFF_4occ_fulldataset.csv")
NNR.AFF=read.csv(file="NNR_AFF_4occ_fulldataset.csv")

##Addition--builder and dweller status

BD=matrix(NA, ngen,nbin)

S=X[,c("genus","Builder_Dweller")]
S.unique=unique(S)

for (i in 1:ngen)
{
	cat("BD", i, "\n")
	S_choice=which(S.unique[,1]==ugen[i])
	BD_choice=as.character(S.unique[S_choice,2])
	for (j in FO[i]:LO[i]) 
	{ 
		if(is.na(AFF[i,j])==TRUE){BD[i,j]=NA} else if(AFF[i,j]==0){BD[i,j]=0} else if(AFF[i,j]==1 && BD_choice=="Builder") {BD[i,j]=1} else if(AFF[i,j]==1 && BD_choice=="Dweller"){BD[i,j]=2} else {BD[i,j]=NA}
		
	}
	
}

B.AFF<-D.AFF<-rep(0,nbin)
for (j in 1:nbin)
{
 B.AFF[j]<-sum(BD[,j]==1,na.rm=TRUE)
 D.AFF[j]<-sum(BD[,j]==2,na.rm=TRUE)
}

write.csv(B.AFF, file="B_AFF_FooteHybrid_4occ_fulldataset.csv")
write.csv(D.AFF, file="D_AFF_FooteHybris_4occ_fulldataset.csv")

#Extinction proportions of reef vs. nonreef.
#For each stage, find genera with given affinity.
#If they have an affinity, that implies they are sampled in the stage in question.
#To determine whether they survive, determine whether LO is after stage in question,
#regardless of subsequent affinity.
PE.R<-PE.NR<-rep(NA,nbin)
P.LO.R<-P.UP.R<-rep(NA,nbin) #lower and upper confidence limits approximating +/- 1 SE
P.LO.NR<-P.UP.NR<-rep(NA,nbin)
for (j in 1:nbin)
{
 n<-sum(AFF[,j]==0,na.rm=TRUE)
 nsurv<-sum(AFF[,j]==0&LO>j,na.rm=TRUE)
 p<-PE.NR[j]<-1-nsurv/n
 
 z<-1
 P.LO.NR[j]<-(2*n*p+z^2-sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
 P.UP.NR[j]<-(2*n*p+z^2+sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))

 n<-sum(AFF[,j]==1,na.rm=TRUE)
 nsurv<-sum(AFF[,j]==1&LO>j,na.rm=TRUE)
 p<-PE.R[j]<-1-nsurv/n
 
 z<-1
 P.LO.R[j]<-(2*n*p+z^2-sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
 P.UP.R[j]<-(2*n*p+z^2+sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
}

#Extinction proportions of builder v dweller
#B and D replace R and NR, 0 for NR, 1 for R-Builder, 2 for R-Dweller
#For each stage, find genera with given affinity.
#If they have an affinity, that implies they are sampled in the stage in question.
#To determine whether they survive, determine whether LO is after stage in question,
#regardless of subsequent affinity.
PE.B<-PE.D<-rep(NA,nbin)
P.LO.B<-P.UP.B<-rep(NA,nbin) #lower and upper confidence limits approximating +/- 1 SE
P.LO.D<-P.UP.D<-rep(NA,nbin)
for (j in 1:nbin)
{
 n<-sum(BD[,j]==2,na.rm=TRUE)
 nsurv<-sum(BD[,j]==2&LO>j,na.rm=TRUE)
 p<-PE.D[j]<-1-nsurv/n
 
 z<-1
 P.LO.D[j]<-(2*n*p+z^2-sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
 P.UP.D[j]<-(2*n*p+z^2+sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))

 n<-sum(BD[,j]==1,na.rm=TRUE)
 nsurv<-sum(BD[,j]==1&LO>j,na.rm=TRUE)
 p<-PE.B[j]<-1-nsurv/n
 
 z<-1
 P.LO.B[j]<-(2*n*p+z^2-sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
 P.UP.B[j]<-(2*n*p+z^2+sqrt(z^2*(z^2-4*n*p*(p-1))))/(2*(n+z^2))
}


Crises=c(11, 21, 25, 31, 35, 37, 44, 47, 56, 57, 62)
Big5=c(11, 21, 37, 44, 67)
OtherCrises=c(25, 31, 35, 47, 56, 57, 62)

#plot extinction rates by reef-nonreef affinity
nmin<-8 #minimum number of genera with a given affinity to include point
xlab<-"Bin number"
ylab<-"Proportion extinct"
pch<-c(22,1, 8, 15, 16)
leg<-c("Reef","Nonreef", "Big 5", "Crisis")
col<-c("blue","darkorange", "black", "blue","darkorange")
x<-c(1:nbin)
xlim<-c(0,nbin+1)
ylim<-c(0,1)
offset<-0.1

y<-PE.R
ylo<-P.LO.R
yup<-P.UP.R
ii<-which(NR.AFF>=nmin)

plot(x[ii],y[ii],pch=pch[1],col=col[1],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Proportional Extinction, Min",nmin,"genera"))
segments(x[ii],ylo[ii],x[ii],yup[ii],col=col[1])
points(Crises, y[Crises], pch=15, col="blue")
points(Big5, y[Big5], pch=8, col="blue", cex=1.5)

x<-c(1:nbin)+offset
y<-PE.NR
ylo<-P.LO.NR
yup<-P.UP.NR
ii<-which(NNR.AFF>=nmin)

points(x[ii],y[ii],pch=pch[2],col=col[2],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim)
segments(x[ii],ylo[ii],x[ii],yup[ii],col=col[2])
points(Crises+offset, y[Crises], pch=16, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)

legend(-2, 0.25, leg, pch=pch, col=col, bty="n")
#legend(xlim[1],ylim[2],leg,pch=pch,col=col,bty="n")




#Do the same, this time plotting points where *both* categories meet the minimum number of genera
x=c(1:nbin)
y<-PE.R
ylo<-P.LO.R
yup<-P.UP.R
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)


plot(x[ii],y[ii],pch=pch[1],col=col[1],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Min",nmin,"genera in *both* categories"))
segments(x[ii],ylo[ii],x[ii],yup[ii],col="blue")
points(Crises, y[Crises], pch=15, col="blue")
points(Big5, y[Big5], pch=8, col="blue", cex=1.5)

x<-c(1:nbin)+offset
y<-PE.NR
ylo<-P.LO.NR
yup<-P.UP.NR
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)

points(x[ii],y[ii],pch=pch[2],col="darkorange",xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim)
segments(x[ii],ylo[ii],x[ii],yup[ii],col="darkorange")
points(Crises+offset, y[Crises], pch=16, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)

legend(-2, 0.25, leg, pch=pch, col=col, bty="n")
#legend(xlim[1],ylim[2],leg,pch=pch,col=c("blue", "darkorange"),bty="n")

#par(pty="s")
nmin<-8
xlab<-"Proportion extinct: reefal genera"
ylab<-"Proportion extinct: non-reef genera"
xlim<-ylim<-c(0,1)
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)
x<-PE.R
y<-PE.NR
plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Min",nmin,"genera in *both* categories"))
segments(0,0,1,1,lty=2)

#plot extinction rates by builder-dweller affinity
nmin<-8 #minimum number of genera with a given affinity to include point
xlab<-"Bin number"
ylab<-"Proportion extinct"
pch<-c(22,1, 8, 15, 16)
leg<-c("Builder","Dweller", "Big 5", "Crisis" )
col<-c("blue","darkorange", "black", "blue")
x<-c(1:nbin)
xlim<-c(0,nbin+1)
ylim<-c(0,1)
offset<-0.1

y<-PE.B
ylo<-P.LO.B
yup<-P.UP.B
ii<-which(B.AFF>=nmin)

plot(x[ii],y[ii],pch=pch[1],col=col[1],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Proportional Extinction, Builders and Dwellers, Min",nmin,"genera"))
segments(x[ii],ylo[ii],x[ii],yup[ii],col=col[1])
points(Crises+offset, y[Crises], pch=15, col="blue")
points(Big5+offset, y[Big5], pch=8, col="blue", cex=1.5)

#points(Crises, y[Crises], pch=16, col="red")


x<-c(1:nbin)+offset
y<-PE.D
ylo<-P.LO.D
yup<-P.UP.D
ii<-which(D.AFF>=nmin)

points(x[ii],y[ii],pch=pch[2],col=col[2],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim)
segments(x[ii],ylo[ii],x[ii],yup[ii],col=col[2])
points(Crises+offset, y[Crises], pch=16, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)

legend(-2, 0.25, leg, pch=pch, col=col, bty="n")
#legend(xlim[1],ylim[2],leg,pch=pch,col=col,bty="n")


#Do the same, this time plotting points where *both* categories meet the minimum number of genera-- need to figure out how to make an equivaent to NR.AFF and NNR.AFF, a B.AFF and D.AFF
x<-c(1:nbin)
y<-PE.B
ylo<-P.LO.B
yup<-P.UP.B
ii<-which(B.AFF>=nmin&D.AFF>=nmin)


plot(x[ii],y[ii],pch=pch[1],col=col[1],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Min",nmin,"genera in *both* categories"))
segments(x[ii],ylo[ii],x[ii],yup[ii],col=col[1])
points(Crises+offset, y[Crises], pch=15, col="blue")
points(Big5+offset, y[Big5], pch=8, col="blue", cex=1.5)

x<-c(1:nbin)+offset
y<-PE.D
ylo<-P.LO.D
yup<-P.UP.D
ii<-which(B.AFF>=nmin&D.AFF>=nmin)

points(x[ii],y[ii],pch=pch[2],col=col[2],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim)
segments(x[ii],ylo[ii],x[ii],yup[ii],col=col[2])
points(Crises+offset, y[Crises], pch=16, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)

legend(-2, 0.25, leg, pch=pch, col=col, bty="n")
#legend(xlim[1],ylim[2],leg,pch=pch,col=col,bty="n")

par(pty="s")
nmin<-8
xlab<-"Proportion extinct: builder genera"
ylab<-"Proportion extinct: dweller genera"
xlim<-ylim<-c(0,1)
ii<-which(B.AFF>=nmin&D.AFF>=nmin)
x<-PE.B
y<-PE.D
plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Min",nmin,"genera in *both* categories"))
segments(0,0,1,1,lty=2)


#Extinction vs. eurytopy, reef vs nonreef. For non-reef genera, proportion of non-reef occurrences ranges 0.5 to 1.0
#Likewise for reef genera and reef occurrences.
#NB: regression statistics can be highly unreliable for small sample sizes.
nmin<-8 #specify minimum number of genera to include for regression.
B.R<-B.NR<-rep(NA,nbin) #regression coefficient; positive means higher risk of extinction with increasing stenotopy, i.e. occurrences closer to 100% in environment of affinity.
p.R<-p.NR<-rep(NA,nbin) #p-value
XX.NR<-YY.NR<-XX.R<-YY.R<-numeric()
for (j in 1:nbin)
{
 ii<-which(AFF[,j]==0)
 if (length(ii)>=nmin)
 {
  x<-PNR[ii,j]
  y<-ifelse(LO[ii]>j,0,1) #0=survive, 1=extinct
  XX.NR<-c(XX.NR,x) #concatenate for one big analysis
  YY.NR<-c(YY.NR,y)
  g<-summary(glm(y~x,family=binomial)) #logistic regression
  if (dim(g$coef)[1]==2&dim(g$coef)[2]==4)
  {
   B.NR[j]<-g$coef[2,1]
   p.NR[j]<-g$coef[2,4]
  }
 }
 ii<-which(AFF[,j]==1)
 cat(j,length(ii),"\n")
 if (length(ii)>=nmin)
 {
  x<-PR[ii,j]
  y<-ifelse(LO[ii]>j,0,1) #0=survive, 1=extinct
  XX.R<-c(XX.R,x) #concatenate for one big analysis
  YY.R<-c(YY.R,y)
  g<-summary(glm(y~x,family=binomial)) #logistic regression
  if (dim(g$coef)[1]==2&dim(g$coef)[2]==4)
  {
  B.R[j]<-g$coef[2,1]
  p.R[j]<-g$coef[2,4]
  }
 }
}
#overall regressions on pooled data
#Results indicate that both types of genera have higher risk of extinction as they become
#more stenotopic (i.e. a higher proportion of occurrences in preferred environment).
g.R<-summary(glm(YY.R~XX.R,family=binomial))
g.NR<-summary(glm(YY.NR~XX.NR,family=binomial))
print(g.R)
print(g.NR)


#And the following plots show that, on average, non-reef genera are more stenotopic.
#plot steno/eurytopy by stage
par(pty="s")
Xtmp<-Ytmp<-rep(0,nbin)
for (j in 1:nbin)
{
 i0<-which(AFF[,j]==0)
 Xtmp[j]<-mean(PNR[i0,j])
 i1<-which(AFF[,j]==1)
 Ytmp[j]<-mean(PR[i1,j])
}
Xtmp<-ifelse(is.finite(Xtmp),Xtmp,NA)
Ytmp<-ifelse(is.finite(Ytmp),Ytmp,NA)
xlim<-ylim<-range(Xtmp,Ytmp,na.rm=TRUE)
xlab<-"Mean stenotopy: nonreef genera"
ylab<-"Mean stenotopy: reef genera"
plot(Xtmp,Ytmp,xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim)
segments(0,0,1,1,lty=2)

#Ditto, only for stages where each type represented by minimum number of genera
nmin<-8
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)
plot(Xtmp[ii],Ytmp[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Min",nmin,"genera in *both* categories"))
segments(0,0,1,1,lty=2)


#for stage by stage
xlim<-c(0,nbin+1)
xlab<-"Bin number"
ylab<-"Mean stenotopy"

nmin<-8 #minimum number of genera with given affinity to include
x<-c(1:nbin)
y<-Ytmp
ii<-which(NR.AFF>=nmin)

#plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,pch=pch[1],col=col[1],main=paste("Min",nmin,"genera"))
plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,pch=pch[1],col="blue",main=paste("Min",nmin,"genera"))
points(Crises+offset, y[Crises], pch=15, col="blue")
points(Big5+offset, y[Big5], pch=8, col="blue", cex=1.5)

y<-Xtmp
ii<-which(NNR.AFF>=nmin)
#points(x[ii],y[ii],pch=pch[2],col=col[2])
points(x[ii],y[ii],pch=pch[2],col="darkorange")
points(Crises+offset, y[Crises], pch=15, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)

legend(xlim[1],ylim[1],c("Reefal", "Nonreefal", "Big5", "Crisis"),pch=pch,col=col,bty="n",yjust=0)

#Do the same for stages where *both* categories meet minimum number of genera
x<-c(1:nbin)
y<-Ytmp
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)

plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,pch=pch[1],col=col[1],main=paste("Min",nmin,"genera in *both* categories"))
points(Crises+offset, y[Crises], pch=16, col="blue")
points(Big5+offset, y[Big5], pch=8, col="blue", cex=1.5)

y<-Xtmp
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)
points(x[ii],y[ii],pch=pch[2],col=col[2])
points(Crises+offset, y[Crises], pch=15, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)

legend(xlim[1],ylim[1],c("Reefal", "Nonreefal", "Big5", "Crisis"),pch=pch,col=col,bty="n",yjust=0)



#Do the same for geographic range, just using pooled data.
#Here 0=restricted, 1=widespread, so negative coefficient means restricted genera have higher extinction risk.

GR=read.csv(file="GeogRange_FooteHybrid_4occ.csv")
AFF=read.csv(file="Affinity_FooteHybrid_4occ.csv")

YY.R<-XX.R<-YY.NR<-XX.NR<-numeric()
for (j in 1:nbin)
{
 ii<-which(AFF[,j]==0)
 x<-GR[ii,j]
 y<-ifelse(LO[ii]>j,0,1) #0=survive, 1=extinct
 XX.NR<-c(XX.NR,x) #concatenate for one big analysis
 YY.NR<-c(YY.NR,y)
 
 ii<-which(AFF[,j]==1)
 x<-GR[ii,j]
 y<-ifelse(LO[ii]>j,0,1) #0=survive, 1=extinct
 XX.R<-c(XX.R,x) #concatenate for one big analysis
 YY.R<-c(YY.R,y)
}

g2.R<-summary(glm(YY.R~XX.R,family=binomial))
g2.NR<-summary(glm(YY.NR~XX.NR,family=binomial))
print(g2.R)
print(g2.NR)

#proportion restricted vs. widespread
par(pty="s")
Xtmp<-Ytmp<-rep(0,nbin)
for (j in 1:nbin)
{
 i0<-which(AFF[,j]==0)
 Xtmp[j]<-sum(GR[i0,j]==0)/length(i0)
 i1<-which(AFF[,j]==1)
 Ytmp[j]<-sum(GR[i1,j]==0)/length(i1)
}
Xtmp<-ifelse(is.finite(Xtmp),Xtmp,NA)
Ytmp<-ifelse(is.finite(Ytmp),Ytmp,NA)
xlim<-ylim<-range(Xtmp,Ytmp,na.rm=TRUE)
xlab<-"Proportion geographically restricted genera: nonreef"
ylab<-"Proportion geographically restricted genera: reef"
plot(Xtmp,Ytmp,xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim)
segments(0,0,1,1,lty=2)
#Ditto, only for stages where each type represented by minimum number of genera
nmin<-8
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)
plot(Xtmp[ii],Ytmp[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,main=paste("Proportion geographic restriction (Min",nmin,"genera in *both* categories)"))
segments(0,0,1,1,lty=2)


par(pin=c(6,4))
dev.off()
xlim<-c(0,nbin+1)
xlab<-"Bin number"
ylab<-"Proportion geographically restricted genera"

nmin<-8 #minimum number of genera with given affinity to include
x<-c(1:nbin)
y<-Ytmp
ii<-which(NR.AFF>=nmin)

plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,pch=pch[1],col="blue",main=paste("Min",nmin,"genera"))
points(Crises+offset, y[Crises], pch=16, col="blue")
points(Big5+offset, y[Big5], pch=8, col="blue", cex=1.5)

y<-Xtmp
ii<-which(NNR.AFF>=nmin)
points(x[ii],y[ii],pch=pch[2],col="darkorange")
points(Crises+offset, y[Crises], pch=15, col="darkorange")
points(Big5+offset, y[Big5], pch=8, col="darkorange", cex=1.5)
legend(xlim[1],ylim[2],c("Reefal", "Nonreefal", "Big 5", "Crisis"),pch=pch,col=col,bty="n",yjust=1)

#Do the same for stages where *both* categories meet minimum number of genera
x<-c(1:nbin)
y<-Ytmp
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)

plot(x[ii],y[ii],xlab=xlab,ylab=ylab,xlim=xlim,ylim=ylim,pch=pch[1],col="blue",main=paste("Min",nmin,"genera in *both* categories"))

y<-Xtmp
ii<-which(NR.AFF>=nmin&NNR.AFF>=nmin)
points(x[ii],y[ii],pch=pch[2],col="darkorange")
legend(xlim[1],ylim[2],c("Reefal", "Nonreefal"),pch=pch,col=col,bty="n",yjust=1)

boxplot(Ytmp[ii], Xtmp[ii])
par(mar=c(2,2,2,5))
boxplot(Ytmp[ii]-Xtmp[ii], ylim=c(-1,1), main=" ", ylab="Residuals from 1:1, proportion \n geographically-restricted genera", axes=FALSE)
axis(4)
mtext("Residuals from 1:1,\nProportion geographically-restricted genera", side=4, line=3.5)
median(Ytmp[ii]-Xtmp[ii])

dev.off()
system("open TierneyPlots.pdf")


#output
Genus<-ugen
output<-cbind(Genus,AFF)
write.csv(output,"Affinities.csv",row.names=FALSE)
output<-cbind(Genus,NR)
write.csv(output,"ReefOccurrences.csv",row.names=FALSE)
output<-cbind(Genus,NNR)
write.csv(output,"NonreefOccurrences.csv",row.names=FALSE)
output<-cbind(PE.R,P.LO.R,P.UP.R)
write.csv(output,"ReefExtinctionProportions.csv")
output<-cbind(PE.NR,P.LO.NR,P.UP.NR)
write.csv(output,"NoneefExtinctionProportions.csv")

system("rm PooledLogisticRegressionResults.txt")
sink("PooledLogisticRegressionResults.txt")
cat("Stenotopy vs. Extinction Risk: Reefal")
print(g.R)
cat("Stenotopy vs. Extinction Risk: Non-Reefal")
print(g.NR)
cat("Geographic range vs. Extinction Risk: Reefal")
print(g2.R)
cat("Geographic range vs. Extinction Risk: Non-Reefal")
print(g2.NR)
sink()
system("open PooledLogisticRegressionResults.txt")
 
#################
#################End Foote-based code


###Crossplots
##Reef vs nonreef
#plot 1:1 line, drop prop extn of either on there, with errors.

par(mfrow=c(1,2))

Reef.Transition.NonCrises= subset(Reef.Transition, !(Reef.Transition[,2] %in% c(Crises,Big5)))
NonReef.Transition.NonCrises= subset(NonReef.Transition, !(NonReef.Transition[,2] %in% c(Crises,Big5)))
plot(x=Reef.Transition.NonCrises[,6], y=NonReef.Transition.NonCrises[,6],pch=20,col="blue", xlim=c(0,1), ylim=c(0,1), xlab="Proportional Extinction, Reefal Taxa", ylab="Proportional Extinction Nonreefal Taxa", main="Non-Crisis, Non-Big 5 Stage Boundaries")
segments(0,0,1,1,lty=2)
arrows(Reef.Transition.NonCrises[,6], NonReef.Transition.NonCrises[,12], Reef.Transition.NonCrises[,6], NonReef.Transition.NonCrises[,15], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Transition.NonCrises[,12], NonReef.Transition.NonCrises[,6], Reef.Transition.NonCrises[,15],NonReef.Transition.NonCrises[,6], length=0.05, angle=90, code=3, col="grey")
points(x=Reef.Transition.NonCrises[,6], y=NonReef.Transition.NonCrises[,6],pch=20,col="blue")
legend(0, 1, c("Non-Crisis Stage Boundary"), pch=c(20), col=c("blue"), cex=.9)




plot(x=Reef.Transition.Crises[,6], y=NonReef.Transition.Crises[,6],pch=20,col="blue", xlim=c(0,1), ylim=c(0,1), xlab="Proportional Extinction, Reefal Taxa",ylab=" ", main="Crisis and Big 5 Extinction Boundaries")
segments(0,0,1,1,lty=2)
arrows(Reef.Transition.Crises[,6], NonReef.Transition.Crises[,12], Reef.Transition.Crises[,6], NonReef.Transition.Crises[,15], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Transition.Crises[,12], NonReef.Transition.Crises[,6], Reef.Transition.Crises[,15],NonReef.Transition.Crises[,6], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Transition.Big5[,6], NonReef.Transition.Big5[,12], Reef.Transition.Big5[,6], NonReef.Transition.Big5[,15], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Transition.Big5[,12], NonReef.Transition.Big5[,6], Reef.Transition.Big5[,15],NonReef.Transition.Big5[,6], length=0.05, angle=90, code=3, col="grey")
points(x=Reef.Transition.Crises[,6], y=NonReef.Transition.Crises[,6],pch=20,col="blue")
points(x=Reef.Transition.Big5[,6], y=NonReef.Transition.Big5[,6], pch=5, col="blue")


legend(0, 1, c("Reef Crisis", "Big 5 Extinction"), pch=c(20,5), col=c("blue", "blue"), cex=.9)


##Extinction, Builders versus dwellers

par(mfrow=c(1,2))

Reef.Builder.Transition.NonCrises= subset(Reef.Builder.Transition, !(Reef.Builder.Transition[,2] %in% c(Crises,Big5)))
Reef.Dweller.Transition.NonCrises= subset(Reef.Dweller.Transition, !(Reef.Dweller.Transition[,2] %in% c(Crises,Big5)))
plot(x=Reef.Builder.Transition.NonCrises[,6], y=Reef.Dweller.Transition.NonCrises[,6],pch=20,col="blue", xlim=c(0,1), ylim=c(0,1), xlab="Proportional Extinction, Reef-builder Taxa", ylab="Proportional Extinction Reef-dweller Taxa", main="Non-Crisis, Non-Big 5 Stage Boundaries")
segments(0,0,1,1,lty=2)
arrows(Reef.Builder.Transition.NonCrises[,6], Reef.Dweller.Transition.NonCrises[,12], Reef.Builder.Transition.NonCrises[,6], Reef.Dweller.Transition.NonCrises[,15], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.NonCrises[,12], Reef.Dweller.Transition.NonCrises[,6], Reef.Builder.Transition.NonCrises[,15], Reef.Dweller.Transition.NonCrises[,6], length=0.05, angle=90, code=3, col="grey")
points(x=Reef.Builder.Transition.NonCrises[,6], y=Reef.Dweller.Transition.NonCrises[,6],pch=20,col="blue")
legend(0, 1, c("Non-Crisis Stage Boundary"), pch=c(20), col=c("blue"), cex=.9)


plot(x=Reef.Builder.Transition.Crises[,6], y=Reef.Dweller.Transition.Crises[,6],pch=20,col="blue", xlim=c(0,1), ylim=c(0,1), xlab="Proportional Extinction, Reef-builder Taxa",ylab=" ", main="Crisis and Big 5 Extinction Boundaries")
segments(0,0,1,1,lty=2)
arrows(Reef.Builder.Transition.Crises[,6], Reef.Dweller.Transition.Crises[,12], Reef.Builder.Transition.Crises[,6], Reef.Dweller.Transition.Crises[,15], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.Crises[,12], Reef.Dweller.Transition.Crises[,6], Reef.Builder.Transition.Crises[,15],Reef.Dweller.Transition.Crises[,6], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.Big5[,6], Reef.Dweller.Transition.Big5[,12], Reef.Builder.Transition.Big5[,6], Reef.Dweller.Transition.Big5[,15], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.Big5[,12], Reef.Dweller.Transition.Big5[,6], Reef.Builder.Transition.Big5[,15], Reef.Dweller.Transition.Big5[,6], length=0.05, angle=90, code=3, col="grey")
points(x=Reef.Builder.Transition.Crises[,6], y=Reef.Dweller.Transition.Crises[,6],pch=20,col="blue")
points(x=Reef.Builder.Transition.Big5[,6], y=Reef.Dweller.Transition.Big5[,6], pch=5, col="blue")


legend(0, 1, c("Reef Crisis", "Big 5 Extinction"), pch=c(20,5), col=c("blue", "blue"), cex=.9)




##Persistence,Builders versus Dwellers

par(mfrow=c(1,2))

Reef.Builder.Transition.NonCrises= subset(Reef.Builder.Transition, !(Reef.Builder.Transition[,2] %in% c(Crises,Big5)))
Reef.Dweller.Transition.NonCrises= subset(Reef.Dweller.Transition, !(Reef.Dweller.Transition[,2] %in% c(Crises,Big5)))
plot(x=Reef.Builder.Transition.NonCrises[,8], y=Reef.Dweller.Transition.NonCrises[,8],pch=20,col="blue", xlim=c(0,1), ylim=c(0,1), xlab="Proportion Reef Persistence, Reef-builder Taxa", ylab="Proportion Reef Persistence, Reef-dweller Taxa", main="Non-Crisis, Non-Big 5 Stage Boundaries")
segments(0,0,1,1,lty=2)
arrows(Reef.Builder.Transition.NonCrises[,8], Reef.Dweller.Transition.NonCrises[,14], Reef.Builder.Transition.NonCrises[,8], Reef.Dweller.Transition.NonCrises[,17], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.NonCrises[,14], Reef.Dweller.Transition.NonCrises[,8], Reef.Builder.Transition.NonCrises[,17], Reef.Dweller.Transition.NonCrises[,8], length=0.05, angle=90, code=3, col="grey")
points(x=Reef.Builder.Transition.NonCrises[,8], y=Reef.Dweller.Transition.NonCrises[,8],pch=20,col="blue")
legend(-.02, .98, c("Non-Crisis Stage Boundary"), pch=c(20), col=c("blue"), cex=.76)


plot(x=Reef.Builder.Transition.Crises[,8], y=Reef.Dweller.Transition.Crises[,8],pch=20,col="blue", xlim=c(0,1), ylim=c(0,1), xlab="Proportional Reef Persistence, Reef-builder Taxa",ylab=" ", main="Crisis and Big 5 Extinction Boundaries")
segments(0,0,1,1,lty=2)
arrows(Reef.Builder.Transition.Crises[,8], Reef.Dweller.Transition.Crises[,14], Reef.Builder.Transition.Crises[,8], Reef.Dweller.Transition.Crises[,17], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.Crises[,14], Reef.Dweller.Transition.Crises[,8], Reef.Builder.Transition.Crises[,17],Reef.Dweller.Transition.Crises[,8], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.Big5[,8], Reef.Dweller.Transition.Big5[,14], Reef.Builder.Transition.Big5[,8], Reef.Dweller.Transition.Big5[,17], length=0.05, angle=90, code=3, col="grey")
arrows(Reef.Builder.Transition.Big5[,14], Reef.Dweller.Transition.Big5[,8], Reef.Builder.Transition.Big5[,17], Reef.Dweller.Transition.Big5[,8], length=0.05, angle=90, code=3, col="grey")
points(x=Reef.Builder.Transition.Crises[,8], y=Reef.Dweller.Transition.Crises[,8],pch=20,col="blue")
points(x=Reef.Builder.Transition.Big5[,8], y=Reef.Dweller.Transition.Big5[,8], pch=5, col="blue")