Nothing
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.