# Main function, computes both the SSAS (Sexual Segregation and # Aggregation Statistic) and the 95% limits of SSAS # under the assumption of random association of both sexes SSAS <- function(x, conf.int = 0.95, B = 10000) { x <- as.matrix(x) nr <- nrow(x) nc <- ncol(x) sr <- rowSums(x) sc <- colSums(x) n <- sum(x) E <- outer(sr, sc, "*")/n dimnames(E) <- dimnames(x) tmp <- .C("chisqsim", as.integer(nr), as.integer(nc),as.integer(sr), as.integer(sc), as.integer(n), as.integer(B), as.double(E), integer(nr * nc), double(n + 1), integer(nc), results = double(B), PACKAGE = "stats") obs <- sum(sort((x - E)^2/E, decreasing = TRUE))/n sim <- tmp\$results/n p0 <- (1 - conf.int)/2 return(c(obs, quantile(sim, p0), quantile(sim, 1 -p0))) } # This function formats data to be run with the SSAS function splitmfd <- function(mfd) { loca1 <- function(x) { x <- t(x[, 1:2]) dimnames(x) <- list(c("mal", "fem"), as.character(1:ncol(x))) x } l0 <- split(mfd, mfd\$mon) lapply(l0, loca1) } # Example 1: Isard rup <- read.table("http://pbil.univ-lyon1.fr/R/donnees/mfdrupicapra.txt", h = T) # Load data from the web plot1 <- function(w, titre = "") { plot(1:12, w[, 1], ylim = range(w), ax = F, pch = 19, type = "n", ylab = "IK", xlab = "") title(main = titre) box() axis(1, 1:12, c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) axis(2, pretty(range(w)), tck = 1) polyx <- c(1:12, 12:1) polyy <- c(w[, 3], rev(w[, 2])) polygon(polyx, polyy, col = grey(0.9)) points(w[, 1], pch = 19, type = "b") } # Function to plot data and 95\ levels of significance l1 <- splitmfd(rup) # Format data to be used with SSAS function w <- matrix(unlist(lapply(l1, SSAS)), ncol = 3, byrow = T) # "w" is a matrix having 3 columns and 12 rows. In the first columns are # the SSAS estimates for each month, and the lower and upper limits in columns 2 # and 3 respectively. plot1(w, "Isard") # Plot figure 3a # Example 2: Red deer cer <- read.table("http://pbil.univ-lyon1.fr/R/donnees/mfdcervus.txt", h = T) l1 <- splitmfd(cer) w <- matrix(unlist(lapply(l1, SSAS)), ncol = 3, byrow = T) plot1(w, "Red deer") # Plot figure 3c # Example 3: Roe deer cap <- read.table("http://pbil.univ-lyon1.fr/R/donnees/mfdcapreolus.txt", h = T) l1 <- splitmfd(cap) w <- matrix(unlist(lapply(l1, SSAS)), ncol = 3, byrow = T) plot1(w, "Roe deer") # Plot figure 3e