|
The following R functions calculate the degree of anthropocentricism in cladograms. The application of the functions is explained and discussed in the paper "Anthropocentricisms in cladograms" (Sandvik 2009). These functions are not part of any R package which you can install. However, simply pasting the functions into your R console should work without problems. There won’t be any help functions available. Instead, each functions starts with some comment lines that explain its purpose and parameters. Simply contact me in case of any trouble.
"%=%" <- function(arg1, arg2) identical(all.equal(arg1, arg2), TRUE)
# This auxiliary function is used by some of the others. It tests
# whether the variables on each side of the equality sign are
# identical (within meaningful accuracy bonds). The difference
# from the default function "==" is that "%=%" compares vectors
# as a whole, not element-wise, and that it ignores rounding
# errors.
Pascal <- function(row, i)
{# This auxiliary function returns a part of Pascal’s triangle.
# The part is specified by the variables "row" and "i":
# - "row" specifies a single row of Pascal’s triangle, where the
# first row is 1, the second row is c(1, 1), the third is
# c(1, 2, 1), etc.
# - "i" specifies the element(s) within this row that shall be
# returned. It may be a single integer, or a vector of integers.
if (missing(row) || !is.numeric(row) || length(row) != 1
|| is.na(row) || row <= 0 || missing(i) || !length(i)
|| !is.numeric(i) || any(is.na(i)) || any(i < 0))
{P <- NA
cat("Sorry, the function Pascal requires positive integers ",
"as input.\n", sep="")}
else
{if (!is.vector(i))
{i <- as.vector(i)}
if (i %=% 1 | i %=% row)
{P <- 1}
else
{if (all(i > row))
{P <- rep(0, length(i))}
else
{len <- min(row, max(i)) + 1
P <- rep(0, len)
P[2] <- 1
if (row > 1)
{for (j in 2:row)
{P <- c(0, P[1:(len-1)] + P[2:len])}}
P <- P[i + 1]
P[which(is.na(P))] <- 0}}}
P}
branches <- function(tips, side)
{# This auxiliary function returns the number of possible tree
# topologies for a given number of tips and a given number of
# side branches to the focal taxon. In counting different tree
# topologies, trees are regarded as rooted, unlabeled,
# bifurcating and non-mirrorable. The parameter "tips" has to
# be a single integer, while "side" (the number of side branches)
# can be single integer or a vector. Its values have to be in the
# range [1; tips - 1].
if (missing(tips) || !is.numeric(tips) || length(tips) != 1
|| is.na(tips) || tips <= 0 || missing(side) || !length(side)
|| !is.numeric(side) || any(is.na(side)) || any(side < 0))
{b <- NA
cat("Sorry, the function \"branches\" requires positive ",
"integers as input.\n", sep="")}
else
{if (!is.vector(side))
{side <- as.vector(side)}
if (all(side >= tips))
{b <- rep(0, length(side))}
else
{if (side %=% (tips - 1))
{b <- 1}
else
{len <- tips + 2
b <- rep(0, len)
b[2] <- 1
if (tips > 2)
{for (j in 3:tips)
{for (k in (len-1):2)
{b[k] <- b[k-1] + b[k+1]}}}
b <- b[side + 1]
b[which(is.na(b))] <- 0}}}
b}
attention <- function(t, s)
{# This function calculates the Human Attention Score of a
# cladogram. Its arguments are:
# - "t", the number of tips of the cladogram
# - "s", the number of branches that are side branches to the
# focal taxon. (Consequently, "s" has to be in the range
# [1; t-1].)
abs(s - log2(t)) /
(ifelse(s > log2(t), t - 1, 1) - log2(t)) / 2 + 0.5}
p.att <- p.attention <- function(t, s)
{# This function calculates the probability accompanying the Human
# Attention Score of a cladogram. The parameters are the same as
# for the function "attention".
g <- branches(t, 1:t)
sum(g[s:length(g)]) / sum(g)}
bs.att <- bs.attention <- function(t, s, nboot=12000, quiet=F)
{# This function bootstraps the Human Attention Scores for a number
# of cladograms and returns the mean Human Attention Score and the
# accompanying probability, both weighted by the number of tips of
# the cladograms and unweighted. The parameters "t" and "s" are
# the same as for the function "attention", but are vectors rather
# than single values. The two vectors have to be of equal length.
# Parameter "nboot" controls the number of bootstrap replicates.
# If "quiet" is TRUE, output showing the progress is suppressed.
f <- length(t)
im <- 0
for (j in 1:f)
{im[j] <- attention(t[j], s[j])}
um <- mean(im)
wm <- weighted.mean(im, t)
ubs <- wbs <- 0
for (i in 1:nboot)
{g <- rep(0, f)
for (j in 1:f)
{x <- runif(1, 1, sum(branches(t[j], 1:t[j])))
l <- 1
while (x > branches(t[j], l))
{x <- x - branches(t[j], l)
l <- l + 1}
g[j] <- attention(t[j], l)}
if (mean(g) >= um)
{ubs <- ubs + 1}
if (weighted.mean(g, t) >= wm)
{wbs <- wbs + 1}
if (!quiet)
{cat("p(unweighted) = ", rund(ubs/i, 4),
", p(weighted) = ", rund(wbs/i, 4),
" (after ", floor(100 * i/nboot),
"% of iterations)\n", sep="")}}
res <- c(um, wm, ubs/nboot, wbs/nboot)
names(res) <- c("unweighted.mean", "weighted.mean",
"unweighted.probability", "weighted.probability")
res}
right <- function(n, r)
{# This function calculates the Human Rightness Score of a
# cladogram. Its arguments are:
# - "n", the number of nodes below the focal taxon (including the
# basalmost and excluding the focal tip itself)
# - "r", the number of these nodes were the focal taxon lies on
# the right-hand side. (Consequently, "r" has to be in the
# range [0; n].)
r / n}
p.right <- function(n, r, tails=1)
{# This function calculates the probability accompanying the Human
# Rightness Score of a cladogram. The parameters "n" and "r" are
# the same as for the function "right". The parameter "tails"
# controls whether the test is one-tailed (default) or two-tailed
# (if tails=2). In the latter case, the probability accompanies
# the Human Utterness Score.
n <- n + 1
Pa <- Pascal(n, 1:n)
if (tails %=% 2)
{p <- 2 * Pa / sum(Pa)}
else
{p <- 0
for (i in 1:n)
{p[i] <- sum((Pa / sum(Pa))[n:i])}}
p[r + 1]}
|