R/addpart.R

In stratigraph: Toolkit for the plotting and analysis of stratigraphic and palaeontological data

addpart <- function(x = NULL, file = 'rawcounts.txt', outfile = NULL){

# Read data in
#  expects raw tab separated count matrix with one header row giving column
#  lables and the first column giving row lables. Note that the header row
#  should have one fewer records than all the other rows.
if(is.null(x)){
  rawcounts <- read.table(file, header=TRUE, row.names=1)
}else{
  rawcounts <- x
}

# abort if there are empty rows or columns in input data
empty.rows <- (1:nrow(rawcounts))[rowSums(rawcounts) == 0]
if(length(empty.rows) > 0)
  stop(paste('empty row or rows in input data:', paste(empty.rows, collapse = ', ')))
empty.cols <- (1:ncol(rawcounts))[colSums(rawcounts) == 0]
if(length(empty.cols) > 0)
  stop(paste('empty col or cols in input data:', paste(empty.cols, collapse = ', ')))

# initialize proportion matrix
propcounts <- rawcounts

# divide each cell by its column sum to give proportional count matrix
for (i in 1:ncol(rawcounts))
  propcounts[,i] <- rawcounts[,i]/colSums(rawcounts)[i]
  
# lambda gets sums of squares of the proportional counts
lambda <- colSums(propcounts^2)

# Simpson's D gets 1-lambda
D <- 1 - lambda

# Hurlbert's PIE gets D*n/(n-1) where n=sample size
PIE <- D * (colSums(rawcounts)/(colSums(rawcounts) - 1))

# BELOW ARE CALCULATIONS RELATED TO S (richness)

# gives a vector of observed S for each sample
sampleSobs <- colSums(rawcounts != 0)

# total richness of pooled samples
St <- length(rowSums(rawcounts) != 0)

# gives a vector of differences between pooled richness (St) and sample richnesses (Sj)
Sdif <- colSums (rawcounts == 0)

# vectors of sample weights (qj) using n=sample size, ln(n) and equal weights
qjn <- colSums(rawcounts) / sum(rawcounts)
qjlnn <- log(colSums(rawcounts)) / sum(log(colSums(rawcounts)))
qjeq <- rep(1 / ncol(rawcounts), ncol(rawcounts))

# calculates weighted B-richness (Sbeta) for each sample (qj)*(St-Sj) using three weightings
siteSbetan <- qjn * Sdif
siteSbetalnn <- qjlnn * Sdif
siteSbetaeq <- qjeq * Sdif

# calculates weighted B-richness (Sbeta) for pooled sample using three weightings
Sbetan <- sum(siteSbetan)
Sbetalnn <- sum(siteSbetalnn)
Sbetaequal <- sum(siteSbetaeq)

# calculates Salpha (St-Sbeta) for pooled sample using three weightings
Salphan <- St-Sbetan
Salphalnn <- St-Sbetalnn
Salphaequal <- St-Sbetaequal

# BELOW ARE CALCULATIONS RELATING TO SIMPSONS D

# calculates pbari (sum of qj*pij for all samples) using three weightings
pbarin <- rep(0,nrow(propcounts)) 
pbarilnn <- rep(0,nrow(propcounts)) 
pbarieq <- rep(0,nrow(propcounts))

for (i in 1:nrow(propcounts))
  pbarin[i] <- sum(propcounts[i,] * qjn)
for (i in 1:nrow(propcounts))
  pbarilnn[i] <- sum(propcounts[i,] * qjlnn)
for (i in 1:nrow(propcounts))
  pbarieq[i] <- sum(propcounts[i,] * qjeq)

sumpbarin <- sum(pbarin^2)
sumpbarilnn <- sum(pbarilnn^2)
sumpbarieq <- sum(pbarieq^2)

# calculates Dbeta for three weightings
Dbetan <- sum(lambda * qjn) - sumpbarin
Dbetalnn <- sum(lambda * qjlnn) - sumpbarilnn
Dbetaeq <- sum(lambda * qjeq) - sumpbarieq

# calculates weighted mean Dalpha for three weightings
meanDalphan <- sum(D*qjn)
meanDalphalnn <- sum(D*qjlnn)
meanDalphaeq <- sum(D*qjeq)

# calculates Dtotal for three weightings
Dtotaln <- meanDalphan+Dbetan
Dtotallnn <- meanDalphalnn+Dbetalnn
Dtotaleq <- meanDalphaeq+Dbetaeq

# prints Salpha (St-Sbeta), Sbeta, percSalpha ((St-Sbeta)/St),
#  percSbeta (Sbeta/St) and St for each weighting which are Salphan,
#  Salphalnn, Salphaequal, Sbetan, Sbetalnn, Sbetaequal, St, percSalpha, percSbeta
counttable <- cbind(rawcounts,Sums=rowSums(rawcounts))
counttable <- rbind(counttable,Sums=colSums(counttable))

Sstats <- rbind(Sobs=c(sampleSobs,St),
                   SimpsonsD=c(D,Dtotaln),
				   PIE=c(PIE, sum(rawcounts)/(sum(rawcounts) - 1) * Dtotaln))

Stotal <- rep(St, 3)
pooledSbeta <- c(Sbetan, Sbetalnn, Sbetaequal)
meanSalphas <- Stotal-pooledSbeta
proportionSalpha <- meanSalphas/(Stotal)
proportionSbeta <- pooledSbeta/(Stotal)
Dtotals <- c(Dtotaln, Dtotallnn, Dtotaleq)
Dalphas <- c(meanDalphan, meanDalphalnn, meanDalphaeq)
Dbetas <- c(Dbetan, Dbetalnn, Dbetaeq)
proportionDalpha <- Dalphas/Dtotals
proportionDbeta <- Dbetas/Dtotals

foo <- rbind(Stotal,meanSalphas,pooledSbeta,proportionSalpha,proportionSbeta,
             Dtotals,Dalphas,Dbetas,proportionDalpha,proportionDbeta)
colnames(foo) <- c('nweight', 'lnnweight', 'eqweight')

returnlist <- list(counttable = counttable, Sstats = Sstats, Saddpart = foo)

if(!is.null(outfile)){
  write.table(returnlist$counttable, file = outfile, append = FALSE)
  sink(outfile, append = TRUE)
  cat('\n')
  sink()
  write.table(returnlist$Sstats, file = outfile, append = TRUE)
  sink(outfile, append = TRUE)
  cat('\n')
  sink()
  write.table(returnlist$Saddpart, file = outfile, append = TRUE)
}

return(returnlist)
}