README.md

In laurasoul/dispeRse: Models For Paleobiogeography

dispeRse

This repository is an in development R package intended to generate new models for detecting phylogenetic signals of endemism and dispersal and is supported by two Short-Term Visiting Scholarships from NESCent (to Laura Soul and Graeme Lloyd).

You can install and load dispeRse into R using the following:

# Install the devtools package from CRAN:
install.packages("devtools")

# Load the devtools package into R:
library(devtools)

# Install the dispeRse package from github:
install_github("laurasoul/dispeRse")

# Load the dispeRse package into R:
library(dispeRse)

Script to perform (and plot) a random walk on a sphere:

# Load libraries:
library(dispeRse)
install.packages("sphereplot", dependencies = TRUE)
library(sphereplot)

# Plot spherical grid:
rgl.sphgrid()

# Start at North Pole:
lonlat <- c(0, 90)

# For 1000 random walk steps:
for(i in 1:1000) {

    # Take a random step and store new coordinates:
    lonlat2 <- unlist(EndPoint(lonlat[1], lonlat[2], runif(1, 0, 360), abs(rnorm(1, 0, 100)))[c(1, 2)])

    # Plot coordinates on sphere in rainbow colour order:
    rgl.sphpoints(lonlat2[1], lonlat2[2], 1, deg=TRUE, col=rainbow(1000)[i], cex=2)

    # Update lonlat ready for next step:
    lonlat <- lonlat2

}

Script to perform (and plot) a random walk on a tree on a sphere in continuous time:

# Add and load the Claddis package into R:
library(devtools)
install_github("graemetlloyd/Claddis")
library(Claddis)

# Load libraries:
library(dispeRse)
install.packages(c("ape", "maps"), dependencies = TRUE)
library(ape)
library(maps)

# Generate random 50-taxon tree:
tree <- rtree(10)

# Run function setting start point at equator-Greenwich Meridian intersection:
out <- TreeWalkerContinuous(tree, slon = 0, slat = 0, niter = 100, steplengthmean = 0, steplengthsd = 1000)

# Plot map:
map()

# For each branch:
for(i in 1:length(out)) {

    # Plot walk along branch:
    lines(out[[i]][2, ], out[[i]][3, ], col=rainbow(length(out))[i])

}

Script to perform (and plot) a random walk on a tree on a sphere in discrete time:

# For discrete function
run <- TreeWalkerDiscrete(b = 0.1, d = 0.05, steps = 50, slon = 0, slat = 0, steplengthsd = 1000)

# Plot map:
map()

# For discrete function
for (i in 1:nrow(run$latitudes)) {

    lines(run$longitudes[i, ], run$latitudes[i, ], col = "red")

}

Script to check how different responses to hitting the edge of a continent affect eveness of distribution.

# Load maps library:
library(maps)

# Set Earth radius in kilometres:
EarthRad <- 6367.4447

# Set number of time steps to run simulation for:
niter <- 100

# Set number of animals to run simulation for:
n_aminals <- 10000

# Set standard deviation of step size for dispersal in kilometres:
stepsize_sd <- 100

# Set centre of continent (long, lat):
continent_centre <- c(0, 0)

# Set radius of continent in kilometres:
continent_radius <- 2000

# Create empty matrix to store start positions (long, lat) of each animal:
start_positions <- matrix(nrow = 0, ncol = 2)

# Create vectors to store distances from continent centre at beginning and end of simlation:
end_distances <- start_distances <- vector(mode="numeric")

# Pick a response type (one of "bounce", "redraw", or "stop")
response <- "bounce"

# For each animal:
for(i in 1:n_aminals) {

    # Pick a random start position from the square in which the circular continent sits:
    start_position <- c(runif(1, -continent_radius / ((2 * pi * EarthRad) / 360), continent_radius / ((2 * pi * EarthRad) / 360)), runif(1, -continent_radius / ((2 * pi * EarthRad) / 360), continent_radius / ((2 * pi * EarthRad) / 360)))

    # Record distance from centre of continent:
    start_distances[i] <- GreatCircleDistanceFromLongLat(continent_centre[1], continent_centre[2], start_position[1], start_position[2])

    # If the starting point is not on the circular continent:
    if(start_distances[i] > continent_radius) {

        # Whilst the starting point is not on the circular continent:
        while(start_distances[i] > continent_radius) {

            # Pick a new starting point:
            start_position <- c(runif(1, -continent_radius / ((2 * pi * EarthRad) / 360), continent_radius / ((2 * pi * EarthRad) / 360)), runif(1, -continent_radius / ((2 * pi * EarthRad) / 360), continent_radius / ((2 * pi * EarthRad) / 360)))

            # Record distance from centre of continent:
            start_distances[i] <- GreatCircleDistanceFromLongLat(continent_centre[1], continent_centre[2], start_position[1], start_position[2])

        }

    }

    # Store point as start position (now we know it is safely on the continent:
    start_positions <- rbind(start_positions, start_position)

}

# Ossify start positions now as they will get overwritten later:
ossified_start_positions <- start_positions

# For each time step:
for(i in 1:niter) {

    # For each organism:
    for(j in 1:n_aminals) {

        # Draw a random bearing for the next move:
        bearing <- runif(1, 0, 360)

        # Draw a random distance for the next move:
        distance <- abs(rnorm(1, 0, stepsize_sd))

        # Calculate the new position of the organism following the move:
        new_position <- EndPoint(start_positions[j, 1], start_positions[j, 2], bearing, distance)[c("long", "lat")]

        # If draw means leaving continent:
        if(GreatCircleDistanceFromLongLat(new_position$long, new_position$lat, continent_centre[1], continent_centre[2]) > continent_radius) {

            # If response to leaving continent is to redraw:
            if(response == "redraw") {

                # As long as the draw removes the taxon from the continent:
                while(GreatCircleDistanceFromLongLat(new_position$long, new_position$lat, continent_centre[1], continent_centre[2]) > continent_radius) {

                    # Redraw the bearing:
                    bearing <- runif(1, 0, 360)

                    # Redraw the distance:
                    distance <- abs(rnorm(1, 0, stepsize_sd))

                    # Calculate the new position:
                    new_position <- EndPoint(start_positions[j, 1], start_positions[j, 2], bearing, distance)[c("long", "lat")]

                }

                # Overwrite start positions:
                start_positions[j, ] <- unlist(new_position)

                # Record distance from centre of continent:
                end_distances[j] <- GreatCircleDistanceFromLongLat(start_positions[j, 1], start_positions[j, 2], continent_centre[1], continent_centre[2])

            # If response to leaving continent is not to redraw (i.e., is to bounce off of, or stop at, continent edge):
            } else {

# FIND COLLISION POINT STARTS HERE:

                # Scalar which describes the proportion of the disatnce the organism has to travel before colliding with the continent edge:
                distance_modifier <- 0.5

                # Starting stepsize (will shrink as answer is honed in on more precisely):
                stepsize <- 0.1

                # Get starting new position (that will eb overwritten until the collision point is found):
                new_position <- EndPoint(start_positions[j, 1], start_positions[j, 2], bearing, distance_modifier * distance)[c("long", "lat")]

                # Whilst the collision point has not been found:
                while(all.equal(as.vector(GreatCircleDistanceFromLongLat(continent_centre[1], continent_centre[2], new_position$long, new_position$lat, EarthRad = EarthRad, Warn = FALSE)), continent_radius) != TRUE) {

                    # Record current distance of new position from centre of continent:
                    current_distance_from_centre <- as.vector(GreatCircleDistanceFromLongLat(continent_centre[1], continent_centre[2], new_position$long, new_position$lat, EarthRad = EarthRad, Warn = FALSE))

                    # Do we need to increase the distance modifier value?:
                    if(current_distance_from_centre < continent_radius) {

                        # Create potential better value by adding step size to modifier:
                        limit <- distance_modifier + stepsize

                    # Or do we need to decrease the degree modifier value?:
                    } else {

                        # Create potential better value by subtracting step size from modifier:
                        limit <- distance_modifier - stepsize

                    }

                    # Store the new new position (so we can ask if this is better):
                    new_new_position <- EndPoint(start_positions[j, 1], start_positions[j, 2], bearing, limit * distance)[c("long", "lat")]

                    # Get new distance from continent centre based on potential better modifer:
                    new_distance_from_centre <- as.vector(GreatCircleDistanceFromLongLat(continent_centre[1], continent_centre[2], new_new_position$long, new_new_position$lat, EarthRad = EarthRad, Warn = FALSE))

                    # If new distance is closer to continent radius:
                    if(abs(current_distance_from_centre - continent_radius) > abs(new_distance_from_centre - continent_radius)) {

                        # Update new position:
                        new_position <- new_new_position

                        # Update distance modifier itself:
                        distance_modifier <- limit

                    # If current distance is stil our best estimate:
                    } else {

                        # Shrink the step size so we can hone in closer:
                        stepsize <- stepsize * 0.1

                    }

                }

# FIND COLLISION POINT ENDS HERE:

                # If response to hitting a continent edge is to bounce off it:
                if(response == "bounce") {

                    # Get the remaining distance of the step:
                    remaining_distance <- as.vector(distance - GreatCircleDistanceFromLongLat(start_positions[j, 1], start_positions[j, 2], new_position$long, new_position$lat))

                    # Get the bearing to the centre of the continent from the colision point:
                    bearing_to_centre <- as.vector(BearingBetweenTwoLongLatPoints(new_position$long, new_position$lat, continent_centre[1], continent_centre[2]))

                    # Get the bearing to the start of the dispersal step from the colision point:
                    bearing_to_start <- as.vector(BearingBetweenTwoLongLatPoints(new_position$long, new_position$lat, start_positions[j, 1], start_positions[j, 2]))

                    # Get the difference between these two bearings:
                    bearing_difference <- min(c(abs(bearing_to_start - bearing_to_centre), 360 - abs(bearing_to_start - bearing_to_centre)))

                    # If bearing difference should be subtracted to get new "bounce" bearing:
                    if(all.equal((bearing_to_centre + bearing_difference) %% 360, bearing_to_start) == TRUE) {

                        # Get new "bounce" bearing:
                        bearing_to_end <- (bearing_to_centre - bearing_difference) %% 360

                    # If bearing difference should be added to get new "bounce" bearing:
                    } else {

                        # Get new "bounce" bearing:
                        bearing_to_end <- (bearing_to_centre + bearing_difference) %% 360

                    }

                    # Record actual new position (after bounce):
                    new_position <- EndPoint(new_position$long, new_position$lat, bearing_to_end, remaining_distance)[c("long", "lat")]

                    # Overwrite start positions:
                    start_positions[j, ] <- unlist(new_position)

                    # Record distance from centre of continent:
                    end_distances[j] <- GreatCircleDistanceFromLongLat(start_positions[j, 1], start_positions[j, 2], continent_centre[1], continent_centre[2])

                    # Little conditional to catch a second collision (does not do what it should, this is just a test to see if it happens):
                    if(end_distances[j] > continent_radius) stop("Fallen off after bouncing!")

                }

                # If response to hitting a continent edge is to stick to it:
                if(response == "stop") {

                    # Overwrite start positions:
                    start_positions[j, ] <- unlist(new_position)

                    # Record distance from centre of continent:
                    end_distances[j] <- GreatCircleDistanceFromLongLat(start_positions[j, 1], start_positions[j, 2], continent_centre[1], continent_centre[2])

                }

            }

        # If draw means staying on continent:
        } else {

            # Overwrite start positions:
            start_positions[j, ] <- unlist(new_position)

            # Record distance from centre of continent:
            end_distances[j] <- GreatCircleDistanceFromLongLat(start_positions[j, 1], start_positions[j, 2], continent_centre[1], continent_centre[2])

        }

    }

}

# Now make plots to se if there is a bias in animal distribution based on response choice

# Plot map at a size just accommodating the circular continent:
map(xlim=c(-continent_radius / ((2 * pi * EarthRad) / 360), continent_radius / ((2 * pi * EarthRad) / 360)), ylim=c(-continent_radius / ((2 * pi * EarthRad) / 360), continent_radius / ((2 * pi * EarthRad) / 360)))

# Plot start positions in "red":
points(ossified_start_positions[, 1], ossified_start_positions[, 2], pch=19, col="red", cex=0.5)

# Plot end positions in blue:
points(start_positions[, 1], start_positions[, 2], pch=19, col="blue", cex=0.5)

# Create a set of 10 concentric rings to use in determining evenness of animal distribution:
concentric_circles <- seq(0, continent_radius, continent_radius / 10)

# Empty vectors to store N animals in each ring at beginning and end and the area of each ring:
start_N_animals_per_ring <- end_N_animals_per_ring <- sphere_areas <- vector(mode="numeric")

# For each ring:
for(i in 2:length(concentric_circles)) {

    # Get the area:
    sphere_areas[(i - 1)] <- SphericalCapArea(concentric_circles[i]) - SphericalCapArea(concentric_circles[(i - 1)])

    # Get the number of animals in the ring at the start of the simulation:
    start_N_animals_per_ring[(i - 1)] <- length(intersect(which(start_distances > concentric_circles[(i - 1)]), which(start_distances <= concentric_circles[i])))

    # Get the number of animals in the ring at the end of the simulation:
    end_N_animals_per_ring[(i - 1)] <- length(intersect(which(end_distances > concentric_circles[(i - 1)]), which(end_distances <= concentric_circles[i])))

}

# Make a barplot of N animals per area in each ring at the start of the simulation:
barplot(start_N_animals_per_ring / sphere_areas, main="Start")

# Make a barplot of N animals per area in each ring at the end of the simulation:
barplot(end_N_animals_per_ring / sphere_areas, main="End")

Animated results of a dispersal simulation run:

# Load libraries (must have ImageMagick installed too):
library(dispeRse)
library(animation)

# Run a simulation (note, will crash if complete extinction occursso so re-run until it works):
out <- DispersalSimulator(N_steps = 1000, organism_multiplier = 1, radius = 2000, continent_speed_mean = 50, continent_speed_sd = 20, organism_step_sd = 1000, start_configuration = "supercontinent", stickiness = 0)

# Start outputting to GIF:
saveGIF({

    # For each time step:
    for(i in 1:1001) {

        # Plot continents:
        MapContinents(cbind(out$continent_positions[, i, 1], out$continent_positions[, i, 2]), radius = 2000)

        # Plot organisms:
        points(x = out$organism_longitudes[, i], y = out$organism_latitudes[, i], pch=19, col="red")

    }

# Additional output information:
}, movie.name = "test.gif", interval = 0.1, nmax = 1001, ani.width = 600, ani.height = 600)

More to follow.

laurasoul/dispeRse documentation built on May 25, 2021, 4:50 a.m.