load(file="Input objects-Western Atlantic.Rdata") #the input file has these objects #3 objects - species, genera, genusnames - vectors of 928 species names, 390 genus names, and 928 genus names for each species #2 objects - Slongextent, Slatextent - empirical latitudinal and longitudinal ranges in degrees for 928 Western Atlantic species #2 objects - lats, longs - latitudes and longitudes of 1-degree cells with species records in the Western Atlantic database #3 objects - shelflats, shelflongs, originationrates - latitudes and longitudes of all 1-degree cells in the Western Atlantic (shelf and island < 200 m) #and origination rates for each cell - based on temperature dependence of speciations, using mean sea-surface temperature for each cell require(permute) require(grid) require(gridBase) Ntaxa<-length(species) Ngenera<-length(genera) SG=tapply(species, genusnames,length) gridlat=seq(-89.5,89.5,by=1) gridlong=seq(-179.5,179.5,by=1) gridmap<-matrix(data=0, nrow=length(gridlat), ncol=length(gridlong)); rownames(gridmap)=gridlat colnames(gridmap)=gridlong #TRUE - model with spatial dependency between range centroids of ancestors and descendants #FALSE - model without spatial dependency between range centroids of ancestors and descendants SPATIALDEP=TRUE #SPECIES RANGE LIMITS Slimit<-rep(NA,Ntaxa);Nlimit<-rep(NA,Ntaxa);Elimit<-rep(NA,Ntaxa);Wlimit<-rep(NA,Ntaxa) #SPECIES LATITUDINAL RANGE Slatrange<-rep(NA,Ntaxa); #ANCESTORS-DESCENDANTS transition=numeric(); descendant=numeric(); ancmidpoint=numeric(); absdescendant=numeric(); absancmidpoint=numeric() #LATITUDINAL LOCATION OF ANCESORS FOR SG DYNAMIC anclatlocation<-rep(NA,Ntaxa) #IDENTIFY CELLS WITH LATITUDINAL LOCATION OF ANCESORS FOR SG DYNAMIC anclocationID<-rep(NA,Ntaxa) #SPECIES SAMPLED LATITUDINAL RANGE sampledSlatrange<-rep(NA,Ntaxa) #MODELED ARRAY WITH PRESENCE-ABSENCE IN EACH CELL FOR CELLS THAT WILL BE SAMPLED simcomp<-array(0, dim=c(length(lats),Ntaxa)) #SOME VECTORS FOR BOOKKEEPING secondlongint<-numeric() addlongitude<-numeric() ###################################################################################################################### #distance represens LENGTHS OF ONE LONGITUDINAL DEGREE IN KM FOR 0.5 to 89.5 latitude longitudelats=seq(0.5,90,by=1) distance=c(111319,111303,111252,111168,111050,110899,110714,110495,110243,109958, 109639,109288,108903,108485, 108034,107550,107034,106486,105905,105292,104647,103970,103262,102523,101752,100950,100118,99255,98362,97439, 96486,95504,94493,93453,92385,91288,90164,89012,87832,86626,85394,84135,82851,81541,80206,78847,77463,76056, 74625,73172,71696,70198,68678,67137,65576,63994,62393,60772,59133,57475,55800,54107, 52398,50673,48932,47176, 45405,43620,41822,40010,38187,36351,34504,32647,30779, 28902,27016,25121,23219,21310,19393,17471,15544,13611, 11675,9735,7791,5846,3898, 1949,0)/1000 ###################################################################################################################### longitudinalweight<-numeric() proporiginationrates=originationrates/max(originationrates,na.rm=T) proplongitudinalweight=distance/max(distance) originationweight<-numeric() for (i in 1:length(shelflats)) { position<-which(abs(shelflats[i])==longitudelats) longitudinalweight[i]<-proplongitudinalweight[position] originationweight[i]<-proporiginationrates[i] } #longitudinalweight accounts for a larger longitudinal component of cells in the tropics finalweight=originationweight*longitudinalweight SGDYNAMIC=TRUE #SHUFFLE SPECIES LATITUDINAL RANGES if (SGDYNAMIC==TRUE) { reshuffle=shuffle(length(Slatextent)) newSlatextent=Slatextent[reshuffle] newSlongextent=Slongextent[reshuffle] } if (SGDYNAMIC==FALSE) { reshuffle=shuffle(length(Slatextent)) newSlatextent=Slatextent[reshuffle] newSlongextent=Slongextent[reshuffle] } ######################################################################################################################### # SIMULATION RUN ###################################################################################################################### #find the locations of first species in each genus for j in 1:indicesofgenera #after these initial species of genera are placed, their progeny is randomly placed again (i.e., omit #indicesofgenera from Ntaxa), one can visualize initial gradient in richness just by looking at initial genera specfactor<-as.integer(factor(species)) genfactor<-as.integer(factor(genera)) #locations of first species in each genus if (SGDYNAMIC==TRUE) {indicesofgenera=tapply(specfactor, genusnames,min)} if (SGDYNAMIC==FALSE) {indicesofgenera=1:Ntaxa} genusassignment=rep(NA, Ntaxa) #this loop is for all initial founding species for each genus for (i in 1:length(indicesofgenera)) { genusassignment[i]=indicesofgenera[i] #RANDOMLY SELECT ONE LATITUDE AND ONE LONGITUDE FROM THE VECTOR OF 180 LONGITUDINAL AND LATITUDINAL VALUES AT 5-DEGREE RESOLUTION #THESE VALUES WILL REPRESENT RANDOMLY LOCATED MIDPOINTS OF SPECIES RANGES idFORxlocation<-sample(c(1:length(shelflongs)),1, replace=T,prob=finalweight) xlocation<-shelflongs[idFORxlocation] ylocation<-shelflats[idFORxlocation] anclatlocation[i]<-ylocation anclocationID[i]<-idFORxlocation #ASSIGN LATITUDINAL AND LONGITUDINAL RADIUS TO EACH SPECIES, USING EMPIRICAL LATITUDINAL AND LONGITUDINAL RANGES DIVIDED BY 2 in degrees yradius=newSlatextent[i]/2 xradius=newSlongextent[i]/2 #IDENTIFY LATITUDINAL AND LONGITUDINAL VECTORS OF 5-DEGREE GRID CELLS #THAT WILL CAPTURE MAXIMUM LATITUDINAL AND LONGITUDINAL RANGE SIZE FOR EACH SPECIES #IDENTIFY NORTHERN AND SOUTHERN RANGE LIMIT Nlimit[i]<-ylocation+yradius Slimit[i]<-ylocation-yradius #IDENTIFY EASTERN AND WESTERN RANGE LIMITS Elimit[i]<-xlocation+xradius Wlimit[i]<-xlocation-xradius #IF RANGE LIMITS ARE BEYOND +-90 DEGREES, TRUNCATE THEM. HOWEVER, KEEP THE UNTRUNCATED LIMITS FOR LATER Nlimit2<-Nlimit[i] Slimit2<-Slimit[i]; if (Slimit[i]<=(-89.5)) { Slimit2<-Slimit[i]; Slimit[i]<--89.5 } if (Nlimit[i]>=90) { Nlimit2<-Nlimit[i]; Nlimit[i]<-89.5 } Slatrange[i]<-Nlimit[i]-Slimit[i] ################################################################################### #CORRECTING FOR TRUNCATION EFFECTS ON THE WESTERN AND EASTERN EDGES #IF RANGE LIMITS ARE BEYOND +-180 DEGREES, TRUNCATE THEM. HOWEVER, KEEP THE #UNTRUNCATED LIMITS FOR LATER COMPUTATION ################################################################################### if (Wlimit[i]>(-180) & Elimit[i]<180) { #EWmake is a vector that identifies the longitudes in the gridmap #longint is a vector that identifies the column locations in the gridmap EWmake<-c(seq(Wlimit[i], Elimit[i],by=1)) longint<-findInterval(EWmake+0.5,gridlong) } #IF LONGITUDINAL RANGE LIMITS ARE TRUNCATED ON THE WESTERN SIDE if (Wlimit[i]<=(-180)) { Wlimit2<-Wlimit[i] Wlimit[i]<--180 remain<-180+(Wlimit2+180) id<-findInterval(remain+0.5, gridlong) if (Wlimit[i]+0.5Elimit[i]) { EWmake<-c(seq(gridlong[id],180,by=1),Wlimit[i]+0.5) } longint<-findInterval((EWmake),gridlong) } #IF LONGITUDINAL RANGE LIMITS ARE TRUNCATED ON THE EASTERN SIDE if (Elimit[i]>=180) { Elimit2<-Elimit[i] Elimit[i]<-180 remain<--180+(Elimit2-180) id<-findInterval(remain+0.5, gridlong); if (Wlimit[i]Elimit[i]) { EWmake<-c(Elimit[i], seq(-179.5,gridlong[id], by=1)) } longint<-findInterval(EWmake,gridlong) } ################################################################################### #CORRECTING FOR TRUNCATION EFFECTS ON THE NORTHERN AND SOUTHERN EDGES ################################################################################### first=findInterval(Slimit[i],gridlat) second=findInterval(Nlimit[i],gridlat) #NSmake is a vector that identifies the longitudes (columns) in the gridmap if (first == second) {NSmake<-gridlat[first]} if (first < second) {NSmake<-seq(gridlat[first],gridlat[second],by=1)} #latint is a vector that identifies the latitudes (rows) in the gridmap latint<-findInterval(NSmake,gridlat) NSint<-latint #NORTHERN EDGE EFFECT #compute how much does the range extend into the other hemisphere, #and generate the vector "secondNSint" that defines S and N limits in the #other #hemisphere if (Nlimit2 >= 90) { secondNlimit<-90-(Nlimit2-90) first=findInterval(secondNlimit+0.5,gridlat) second=findInterval(89.5,gridlat) make<-seq(gridlat[first],gridlat[second],by=1) secondNSint<-findInterval(make, gridlat) firstNSinterval<-min(latint):180 latint<-firstNSinterval #concatenate ranges from both hemispheres into "NSint" that #identifies the latitudes in gridmap NSint<-c(min(latint):180, rev(secondNSint)) } #SOUTHERN EDGE EFFECT if (Slimit2 <= -90) { thirdSlimit<--90-(Slimit2+90) first=findInterval(thirdSlimit+0.5,gridlat) second=findInterval(-89.5,gridlat) make<-seq(gridlat[second],gridlat[first],by=1) secondNSint<-findInterval(make, gridlat) NSint<-c(rev(secondNSint), 1:max(latint)) } ################################################################################### #IDENTIFY LONGITUDINAL MIDPOINT OF RANGES - FURTHER CORRECTING FOR TRUNCATION ################################################################################### if (length (EWmake) == 1) { midlongitude<-findInterval(EWmake,gridlong) } if (length (EWmake) > 1) { midlongitude<-EWmake[((length(EWmake))/2)]; midlongitude<-findInterval(midlongitude,gridlong) } #IDENTIFY LATITUDINAL MIDPOINT OF SPECIES RANGE - IF THE NUMBER OF THE OCCUPIED #CELLS ALONG LATITUDE IS EVEN OR THERE IS JUST ONE CELL OCCUPIED if(length(NSint) == 1) { latitudecenter=c(NSint,NSint) } if(length(NSint) > 1) { latitudecenter<-c(NSint[floor(length(NSint)/2)], NSint [ceiling(length(NSint)/2)]) } #IDENTIFY LATITUDINAL MIDPOINT OF SPECIES RANGE - IF THE NUMBER OF THE OCCUPIED #CELLS ALONG LATITUDE IS NOT EVEN if (!is.integer(length(NSint)/2)) { latitudecenter=c(NSint [ceiling(length(NSint)/2)], NSint [ceiling(length(NSint)/2)]) } ################################################################################### #FILL ALL CELLS IN A GRIDDED MAP WHERE A GIVEN SPECIES #SHOULD OCCUR #IDENTIFY OCCUPIED CELLS (WITHOUT TRUNCATION) ################################################################################### #GRIDDED MAP gridmap<-matrix(data=0, nrow=length(gridlat), ncol=length(gridlong)) rownames(gridmap)=gridlat colnames(gridmap)=gridlong gridmap[latitudecenter[1],longint]<-longint helpint<-longint if (length(helpint) > 1) { if(latitudecenter[1] < 180) { for (f in (latitudecenter[2]+1): max(latint)) { gridmap[f,helpint]<-helpint } } } helpint<-longint #THIS LOOP ADDS CELLS TO EACH LATITUDE if (length(helpint) > 1) { for (f in ((latitudecenter[1]-1)):min(latint)) { gridmap[f,helpint]<-helpint } } ################################################################################### #IDENTIFY OCCUPIED CELLS IN RANGES WITH NORTHERN TRUNCATION ################################################################################### if (Nlimit2 > 90) { for (z in 1:length(longint)) { if (longint[z] <= 180) {secondlongint[z]<-longint[z]+180} if (longint[z] > 180) {secondlongint[z]<-longint[z]-180} } helpint<-secondlongint addlongitude=helpint if (length(helpint) == 1) {gridmap[180,helpint]<-helpint} if (length(helpint) > 1) { for (f in 180:min(secondNSint)) { gridmap[f,helpint]<-helpint } } } ################################################################################### #IDENTIFY OCCUPIED CELLS IN RANGES WITH SOUTHERN TRUNCATION ################################################################################### if (Slimit2 < -90) { for (z in 1:length(longint)) { if (longint[z] <= 180) {secondlongint[z]<-longint[z]+180} if (longint[z] > 180) {secondlongint[z]<-longint[z]-180} } helpint<-secondlongint addlongitude=helpint if (length(helpint) > 1) { for (f in 1:max(secondNSint)) { gridmap[f,helpint]<-helpint } } } #################################################################################### #GENERATE SPECIES-BY-CELL TABLES WITH PRESENCE-ABSENCE OF SPECIES (simcomp) ################################################################################### x<-numeric(); y<-numeric(); for (v in 1:length(lats)) { idx<-which(gridlat==lats[v]) idy<-which(gridlong==longs[v]) if (gridmap[idx,idy] > 0) { simcomp[v,i]=1 y<-append(y,lats[v],length(y)) x<-append(x,longs[v],length(x)) } } ################################################################################### #COMPUTE LATITUDINAL RANGE OF THE SAMPLED SPECIES, CONDITIONED BY THE GRID CELLS #THAT ARE ACTUALLY SAMPLED IN THE EMPIRICAL DATASET ################################################################################### polys<-cbind(x,y) if (nrow(polys) > 1) { sampledSlatrange[i]<-max(y)-min(y) } if (nrow(polys)==1) { sampledSlatrange[i]=0.5 } print(i) } #finish the loop for species i #keep initial latitudinal locations of founding species at time of their origination (placement) ancestorplacements=anclatlocation ancestorplacements=ancestorplacements[!is.na(ancestorplacements)] Slatrange2=Slatrange indicesofSGspecies=c(1:Ntaxa) ################################################################################################# ################################################################################################# #ORIGINATION OF CONGENERIC SPECIES ################################################################################################# ################################################################################################# for (zzz in 2:max(SG)) { #locate ids of ancestors (here, simply the first species in a string of all species within each genus) oldindicesofgenera=tapply(specfactor[indicesofSGspecies], genusnames[indicesofSGspecies],min) anclatitudes=anclatlocation[!is.na(anclatlocation)] anclongitudinalweight=numeric() ancoriginationweight<-numeric() Slatrange2[Slatrange2==0]=0.1 ancSlatrange=Slatrange2/max(Slatrange2, na.rm=T) ancSlatrange=ancSlatrange[!is.na(ancSlatrange)] #anclatitudes is a vector with latitudes of origination of all initial ancestor species (e.g., 390 in Western Atlantic) #loop here over all initial species #HERITABILITY OF RANGE LOCATION for (i in 1:length(anclatitudes)) { #shelflats is a vector of all possible latitudes in latitude-longitude cells in a given transect #the position object finds the cells of the ancestor #THIS DETERMINES WHAT ANCESTORS WILL PRODUCE NEW SPECIES position<-which(anclatitudes[i]==shelflats) ancoriginationweight[i]<-proporiginationrates[position] } #origination probability is weighted by latitudinal range size of ancestor ancfinalweight=ancoriginationweight*ancSlatrange #IF ORIGINATIONS FOLLOW LATITUDINAL GRADIENT, THEN TROPICAL SPECIES WILL PRODUCE SPECIES MORE LIKELY THAN EXTRATROPICAL SPECIES if (length(which(SG>(zzz-1)))>1) { newindicesofSGspecies=sample(c(1:Ntaxa)[oldindicesofgenera],length(which(SG>(zzz-1))),replace=F, prob=ancfinalweight+0.00001) } if (length(which(SG>(zzz-1)))==1) { newindicesofSGspecies=oldindicesofgenera[1] } #locate ids of descendant species belonging to genera with at least zzz species #these species will originate along a latitudinal gradient, not all ancestors have descendants newindicesofgenera=c(1:Ntaxa)[newindicesofSGspecies] #identify 1-degree cells where ancestors originated ancestralLATSids=anclocationID[!is.na(anclocationID)] ancestralLATS=shelflats[ancestralLATSids] ancestralweights=finalweight[ancestralLATSids] #sample descendants according to latitude of ancestors if (length(ancfinalweight)==1) {descendants=oldindicesofgenera} if (length(ancfinalweight)>1) {descendants=newindicesofgenera} #descendants will be ancestors in the next step, in the loop for genera with 3 and more species indicesofSGspecies=newindicesofSGspecies starting=length(genusassignment[!is.na(genusassignment)]) ################################################################################################### anclatlocation=rep(NA,Ntaxa) anclocationID=rep(NA,Ntaxa) ancSlatrange=rep(NA,Ntaxa) Slatrange2=rep(NA,Ntaxa) for (i in (starting+1):(starting+length(descendants))) { ii=i-starting ID=which(descendants[ii]==genusassignment) ancestorSlimit=Slimit[ID[length(ID)]] ancestorNlimit=Nlimit[ID[length(ID)]] #identify cells within the ancestor range #ancestormidpoint is equal to anclatitudes[ID] ancestormidpoint=mean(c(Slimit[ID[length(ID)]],Nlimit[ID[length(ID)]])) #THIS IS THE RANGE OF LATITUDES FROM WHICH NEW LOCATION SHOULD BE DRAWN #latitude of origination is weighted by latitudinal differences in shelf extent identifylats=which(shelflats >= ancestorSlimit & shelflats <= ancestorNlimit) if (SPATIALDEP==TRUE) {identifylats=which(shelflats >= ancestorSlimit & shelflats <= ancestorNlimit)} if (SPATIALDEP==FALSE) {identifylats=1:length(shelflats)} #RANDOMLY SELECT ONE LATITUDE AND ONE LONGITUDE FROM THE VECTOR OF 180 LONGITUDINAL AND LATITUDINAL VALUES AT 5-DEGREE RESOLUTION #THESE VALUES WILL REPRESENT RANDOMLY LOCATED MIDPOINTS OF SPECIES RANGES if (length(identifylats)>1 & SPATIALDEP==TRUE) {idFORxlocation<-sample(identifylats,1, replace=T)} if (length(identifylats)>1 & SPATIALDEP==FALSE) {idFORxlocation<-sample(identifylats,1, replace=T, prob=originationweight)} if (length(identifylats)==1) {idFORxlocation<-identifylats} xlocation<-shelflongs[idFORxlocation] ylocation<-shelflats[idFORxlocation] transition[i]=abs(ylocation)-abs(ancestormidpoint) absdescendant[i]=abs(ylocation) absancmidpoint[i]=abs(ancestormidpoint) descendant[i]=(ylocation) ancmidpoint[i]=(ancestormidpoint) anclatlocation[i]<-ylocation anclocationID[i]<-idFORxlocation #ASSIGN LATITUDINAL AND LONGITUDINAL RADIUS TO EACH SPECIES, USING EMPIRICAL LATITUDINAL AND LONGITUDINAL RANGES DIVIDED BY 2 in degrees yradius=newSlatextent[i]/2 xradius=newSlongextent[i]/2 #IDENTIFY NORTHERN AND SOUTHERN RANGE LIMIT, EASTERN AND WESTERN RANGE LIMITS Nlimit[i]<-ylocation+yradius Slimit[i]<-ylocation-yradius Elimit[i]<-xlocation+xradius Wlimit[i]<-xlocation-xradius #IF RANGE LIMITS ARE BEYOND +-90 DEGREES, TRUNCATE THEM. HOWEVER, KEEP THE #UNTRUNCATED LIMITS FOR LATER Nlimit2<-Nlimit[i] Slimit2<-Slimit[i]; if (Slimit[i]<=(-89.5)) { Slimit2<-Slimit[i]; Slimit[i]<--89.5 } if (Nlimit[i]>=90) { Nlimit2<-Nlimit[i]; Nlimit[i]<-89.5 } Slatrange[i]<-Nlimit[i]-Slimit[i] Slatrange2[i]<-Nlimit[i]-Slimit[i] ################################################################################### #CORRECTING FOR TRUNCATION EFFECTS ON THE WESTERN AND EASTERN EDGES #IF RANGE LIMITS ARE BEYOND +-180 DEGREES, TRUNCATE THEM. HOWEVER, KEEP THE #UNTRUNCATED LIMITS FOR LATER COMPUTATION ################################################################################### if (Wlimit[i]>(-180) & Elimit[i]<180) { #EWmake is a vector that identifies the longitudes in the gridmap #longint is a vector that identifies the column locations in the gridmap EWmake<-c(seq(Wlimit[i], Elimit[i],by=1)) longint<-findInterval(EWmake+0.5,gridlong) } #IF LONGITUDINAL RANGE LIMITS ARE TRUNCATED ON THE WESTERN SIDE if (Wlimit[i]<=(-180)) { Wlimit2<-Wlimit[i] Wlimit[i]<--180 remain<-180+(Wlimit2+180) id<-findInterval(remain+0.5, gridlong) if (Wlimit[i]+0.5Elimit[i]) { EWmake<-c(seq(gridlong[id],180,by=1),Wlimit[i]+0.5) } longint<-findInterval((EWmake),gridlong) } #IF LONGITUDINAL RANGE LIMITS ARE TRUNCATED ON THE EASTERN SIDE if (Elimit[i]>=180) { Elimit2<-Elimit[i] Elimit[i]<-180 remain<--180+(Elimit2-180) id<-findInterval(remain+0.5, gridlong); if (Wlimit[i]Elimit[i]) { EWmake<-c(Elimit[i], seq(-179.5,gridlong[id], by=1)) } longint<-findInterval(EWmake,gridlong) } ################################################################################### #CORRECTING FOR TRUNCATION EFFECTS ON THE NORTHERN AND SOUTHERN EDGES ################################################################################### first=findInterval(Slimit[i],gridlat) second=findInterval(Nlimit[i],gridlat) #NSmake is a vector that identifies the longitudes (columns) in the gridmap if (first == second) {NSmake<-gridlat[first]} if (first < second) {NSmake<-seq(gridlat[first],gridlat[second],by=1)} #latint is a vector that identifies the latitudes (rows) in the gridmap latint<-findInterval(NSmake,gridlat) NSint<-latint #NORTHERN EDGE EFFECT #compute how much does the range extend into the other hemisphere, #and generate the vector "secondNSint" that defines S and N limits in the #other #hemisphere if (Nlimit2 >= 90) { secondNlimit<-90-(Nlimit2-90) first=findInterval(secondNlimit+0.5,gridlat) second=findInterval(89.5,gridlat) make<-seq(gridlat[first],gridlat[second],by=1) secondNSint<-findInterval(make, gridlat) firstNSinterval<-min(latint):180 latint<-firstNSinterval #concatenate ranges from both hemispheres into "NSint" that #identifies the latitudes in gridmap NSint<-c(min(latint):180, rev(secondNSint)) } #SOUTHERN EDGE EFFECT if (Slimit2 <= -90) { thirdSlimit<--90-(Slimit2+90) first=findInterval(thirdSlimit+0.5,gridlat) second=findInterval(-89.5,gridlat) make<-seq(gridlat[second],gridlat[first],by=1) secondNSint<-findInterval(make, gridlat) NSint<-c(rev(secondNSint), 1:max(latint)) } ################################################################################### #IDENTIFY LONGITUDINAL MIDPOINT OF RANGES - FURTHER CORRECTING FOR TRUNCATION ################################################################################### if (length (EWmake) == 1) { midlongitude<-findInterval(EWmake,gridlong) } if (length (EWmake) > 1) { midlongitude<-EWmake[((length(EWmake))/2)]; midlongitude<-findInterval(midlongitude,gridlong) } #IDENTIFY LATITUDINAL MIDPOINT OF SPECIES RANGE - IF THE NUMBER OF THE OCCUPIED #CELLS ALONG LATITUDE IS EVEN OR THERE IS JUST ONE CELL OCCUPIED if(length(NSint) == 1) { latitudecenter=c(NSint,NSint) } if(length(NSint) > 1) { latitudecenter<-c(NSint[floor(length(NSint)/2)], NSint [ceiling(length(NSint)/2)]) } #IDENTIFY LATITUDINAL MIDPOINT OF SPECIES RANGE - IF THE NUMBER OF THE OCCUPIED #CELLS ALONG LATITUDE IS NOT EVEN if (!is.integer(length(NSint)/2)) { latitudecenter=c(NSint [ceiling(length(NSint)/2)], NSint [ceiling(length(NSint)/2)]) } ################################################################################### #FILL ALL CELLS IN A GRIDDED MAP WHERE A GIVEN SPECIES #SHOULD OCCUR #IDENTIFY OCCUPIED CELLS (WITHOUT TRUNCATION) ################################################################################### #GRIDDED MAP gridmap<-matrix(data=0, nrow=length(gridlat), ncol=length(gridlong)) rownames(gridmap)=gridlat colnames(gridmap)=gridlong gridmap[latitudecenter[1],longint]<-longint helpint<-longint if (length(helpint) > 1) { if(latitudecenter[1] < 180) { for (f in (latitudecenter[2]+1): max(latint)) { gridmap[f,helpint]<-helpint } } } helpint<-longint #THIS LOOP ADDS CELLS TO EACH LATITUDE if (length(helpint) > 1) { for (f in ((latitudecenter[1]-1)):min(latint)) { gridmap[f,helpint]<-helpint } } ################################################################################### #IDENTIFY OCCUPIED CELLS IN RANGES WITH NORTHERN TRUNCATION ################################################################################### if (Nlimit2 > 90) { for (z in 1:length(longint)) { if (longint[z] <= 180) {secondlongint[z]<-longint[z]+180} if (longint[z] > 180) {secondlongint[z]<-longint[z]-180} } helpint<-secondlongint addlongitude=helpint if (length(helpint) == 1) {gridmap[180,helpint]<-helpint} if (length(helpint) > 1) { for (f in 180:min(secondNSint)) { gridmap[f,helpint]<-helpint } } } ################################################################################### #IDENTIFY OCCUPIED CELLS IN RANGES WITH SOUTHERN TRUNCATION ################################################################################### if (Slimit2 < -90) { for (z in 1:length(longint)) { if (longint[z] <= 180) {secondlongint[z]<-longint[z]+180} if (longint[z] > 180) {secondlongint[z]<-longint[z]-180} } helpint<-secondlongint addlongitude=helpint if (length(helpint) > 1) { for (f in 1:max(secondNSint)) { gridmap[f,helpint]<-helpint } } } ################################################################################### #CHECK RANGE REJECTION - HERE, ASSESSING JUST NORTHERN AND SOUTHERN COORDINATES AT #GEOGRAPHICAL RANGE MIDPOINT ################################################################################### #IDENTIFY LATITUDE AND LONGITUDE OF NORTHERN AND SOUTHERN LIMITS TO CHECK WHETHER #THEY ARE LOCATED ON SHELF if (Slimit2 < -90) { Sfocal<-paste(gridlat[min(NSint)],gridlong[round(median(addlongitude),0)]) } if (Nlimit2 > 90) { Nfocal<-paste(gridlat[max(NSint)],gridlong[round(median(addlongitude),0)]) } if (Slimit2 >= -90) { Sfocal<-paste(gridlat[min(NSint)],gridlong[midlongitude]) } if (Nlimit2 <= 90) { Nfocal<-paste(gridlat[max(NSint)],gridlong[midlongitude]) } ################################################################################### #GENERATE SPECIES-BY-CELL TABLES WITH PRESENCE-ABSENCE OF SPECIES (simcomp) ################################################################################### x<-numeric(); y<-numeric(); for (v in 1:length(lats)) { idx<-which(gridlat==lats[v]) idy<-which(gridlong==longs[v]) if (gridmap[idx,idy] > 0) { simcomp[v,i]=1 y<-append(y,lats[v],length(y)) x<-append(x,longs[v],length(x)) } } ################################################################################### #COMPUTE LATITUDINAL RANGE OF THE SAMPLED SPECIES, CONDITIONED BY THE GRID CELLS #THAT ARE ACTUALLY SAMPLED IN THE EMPIRICAL DATASET ################################################################################### polys<-cbind(x,y) if (nrow(polys) > 1) { sampledSlatrange[i]<-max(y)-min(y) } if (nrow(polys)==1) { sampledSlatrange[i]=0.5 } print(i) if (i==Ntaxa) {break} } #finish the loop for species i genusassignment[(starting+1):(starting+length(descendants))]=descendants print(zzz) } ################################################################################################# genusancmidpoint=tapply(ancmidpoint, genusassignment, function(x) {x=x[!is.na(x)];x[1]}) genusdescendant=tapply(descendant, genusassignment, function(x) {x=x[!is.na(x)]; max(x, na.rm=T)}) genusabsancmidpoint=tapply(absancmidpoint,genusassignment, function(x) {x=x[!is.na(x)];x[1]}) genusabsdescendant=tapply(absdescendant,genusassignment, function(x) {x=x[!is.na(x)];max(x)}) par(mfcol=c(2,2)) par(mar=c(4,4,2,1)) maxylim=150 hist(Slatrange*111.25, breaks=seq(0,16000,by=250), ylim=c(0,maxylim), col="gray", xlab="Species lat. range (km)", main="", ylab="Number of species") abline(v=median(Slatrange*111.25),lty=1,lwd=2) par(mfcol=c(2,2)) par(mar=c(4,4,2,1)) plot(absancmidpoint,absdescendant, pch=21, bg="gray", xlab="Absolute latitude of ancestor", ylab="Absolute latitude of descendant", main="Species", cex.axis=1,cex.lab=1,cex.main=0.9, ylim=c(0,90), xlim=c(0,90)) lines(c(0,90),c(0,90), lty=1) vp <- baseViewports() pushViewport(vp$inner,vp$figure,vp$plot) #pushViewport(viewport(x=0.75,y=0.35,width=0.5,height=0.5)) pushViewport(viewport(x=0.7,y=0.4,width=0.6,height=0.6)) par(plt=gridPLT(),new=T) hist(absdescendant-absancmidpoint, breaks=seq(-80,80,by=5),main="",yaxt="n", ylab="", xlab="", bty="n", col="gray", xaxt="n") axis(1,at=c(-30,-15,0,15,30),labels=c(-30,-15,0,15,30), tick=FALSE, cex.axis=0.7, mgp=c(0,0,0)) abline(v=0,lty=1,lwd=2) abline(v=median(absdescendant-absancmidpoint, na.rm=T),lty=1,lwd=3) popViewport(4) par(mar=c(4, 4, 2, 1) + 0.1) plot(genusabsancmidpoint,genusabsdescendant, pch=21, bg="gray", xlab="Absolute latitude of ancestor", ylab="Absolute latitude of descendant", main="Genera", cex.axis=1,cex.lab=1,cex.main=0.9,ylim=c(0,90), xlim=c(0,90)) lines(c(0,90),c(0,90), lty=1) vp <- baseViewports() pushViewport(vp$inner,vp$figure,vp$plot) pushViewport(viewport(x=0.7,y=0.4,width=0.6,height=0.6)) par(plt=gridPLT(),new=T) hist(genusabsdescendant-genusabsancmidpoint, breaks=seq(-80,80,by=5), main="", yaxt="n", ylab="", xlab="", bty="n", col="gray", xaxt="n") axis(1,at=c(-30,-15,0,15,30),labels=c(-30,-15,0,15,30), tick=FALSE, cex.axis=0.7, mgp=c(0,0,0)) abline(v=0,lty=1,lwd=2) abline(v=median(genusabsdescendant-genusabsancmidpoint, na.rm=T),lty=1,lwd=3) popViewport(4) par(mar=c(4, 4, 2, 1) + 0.1)