brsim {brsim}  R Documentation 
BrainerdRobinson similarity coefficient matrix calculation, with permutationbased pvalues and optional clustering
Description
This function calculates the BrainerdRobinson (BR) similarity coefficient for each pair of row
of the input table (RobinsonBrainerd 1951, 1952).
It also performs a permutation test to assess the significance of each BR coefficient (DeBoerKintighRostoker 1996), and allows to carry
out a hierarchical agglomerative clustering. An optimal cluster solution can be established using the silhouette method
(see details provided below). The function produces a correlation matrix in tabular form, which is also visually plotted as an heatmap.
In the heatmap (which is built using the corrplot
package), the size and the color of the
squares are proportional to the BrainerdRobinson coefficients. Optionally, the heatmap can be reordered on the basis
of the hierachical clustering, with clusters enclosed by red rectangles.
Visit this LINK to access the package's vignette.
Usage
brsim(
df,
num.perm = 1000,
clust = FALSE,
part = NULL,
aggl.meth = "ward.D2",
sort.map = FALSE,
number.cex = 0.7,
cex.dndr.lab = 0.7,
cex.sil.lab = 0.7,
cex.dot.plt.lab = 0.7,
oneplot = TRUE
)
Arguments
df 
A table (dataframe format) where each row represents an assemblage and each column represents an item. 
num.perm 
A numeric value indicating the number of permutations to perform in each test (default is 1000). 
clust 
TRUE (default) or FALSE if the user does or does not want a agglomerative hierarchical clustering to be performed. 
part 
Desired number of clusters; if NULL (default), an optimal partition is calculated (see Details). 
aggl.meth 
Agglomeration method ("ward.D2" by default) to be used; the selected method is internally used for the reordering
of the heatmap if 
sort.map 
TRUE or FALSE (default) if the user does or does not want the rendered heatmap to be ordered on the basis of the selected hierachical clustering. 
number.cex 
Numeric. Set the size of the labels used for the coefficients displayed in the rendered heatmap. 
cex.dndr.lab 
Numeric. Set the size of the labels used in the dendrogram. 
cex.sil.lab 
Numeric. Set the size of the labels used in the silhouette plot. 
cex.dot.plt.lab 
Numeric. Set the size of the labels used in the Cleveland's dotplots representing bycluster proportions. 
oneplot 
TRUE (default) or FALSE if the user wants or does not want the plots to be visualized in a single window. 
Details
Permutationbased pvalues
The rationale behind the pvalue calculation is as follows: for each pair of assemblages in the input data, the function first calculates the observed BrainerdRobinson (BR) coefficient. This is a measure of the similarity between the two assemblages.
The function then performs a certain number of permutations (default is 1000). In each permutation, it generates two new assemblages (each featuring a sample size corresponding to the size of each assemblage being compared) by randomly sampling from the global pool (the combined data of all assemblages), and calculates the BR coefficient for this new pair of assemblages (see DeBoerKintighRostoker 1996). This creates a distribution of BR coefficients that we would expect to see by chance.
The pvalue is then calculated as the proportion of the permuted BR coefficients that are less than or equal
to the observed BR coefficient. A small pvalue (typically < 0.05) suggests that the observed similarity
between the two assemblages is statistically significant; it is unlikely to have occurred just by chance.
In simple terms, the pvalue calculation is trying to answer the question: if there were no real
similarity between these two assemblages, what is the probability that I would observe a similarity as extreme as
(or more extreme than) the one I actually observed, just by chance?
The pvalues are returned in two matrices: in the first, the pvalues are reported as they are, whereas in the second
they are classified as <0.05, <0.01, <0.001, or not significant.
Hierarchical agglomerative clustering
By setting the parameter clust
to TRUE
, the units (rows) for which the BR coefficients have been
calculated will be clustered. Note that the clustering is based on a dissimilarity matrix which
is internally calculated as the maximum value of the BR coefficient (200) minus the observed BR coefficient.
This allows a simpler reading of the dendrogram which is produced by the function, where the less dissimilar (i.e., more similar) units will be placed at lower levels, while more dissimilar (i.e., less similar) units will be placed at higher levels within the dendrogram. The latter depicts the hierarchical clustering based (by default) on the Ward's agglomeration method; rectangles identify the selected cluster partition.
Optionally, by setting the sort.map
to TRUE
, the heatmap can be reordered on the basis of the hierarchical clustering,
with clusters indicated by red rectangles. The number of clusters indicated depends on what requested by the user (see
the next section). Note that, internally, the reordering is based on the same agglomeration method
selected by the user via the aggl.method
parameter, which is set to ward.D2
by default.
Number of clusters and silhouette method
Besides the dendrogram, a silhouette plot is produced, which allows to measure how 'good' is the selected cluster solution.
If the parameter part
is left empty (default), an optimal cluster solution is
obtained.
The optimal partition is selected via an iterative procedure which identifies at which cluster solution the highest average silhouette width is achieved. The cluster solution ranges from a minimum of 2 to a maximum which is equal to the number of units (i.e., the rows of the input dataset) minus 1. The number of units essentially represents the maximum number of clusters that could potentially be formed if each row were its own cluster. However, since a cluster solution requires at least two groups, the maximum number of meaningful clusters is one less than the number of rows.
If a userdefined partition is needed, the user can input the desired number of clusters using the parameter part
.
In either case, an additional plot is returned besides the cluster dendrogram and the silhouette plot; it displays a scatterplot in which the cluster solution (xaxis) is plotted against the average silhouette width (yaxis). A black dot represents the partition selected either by the iterative procedure or by the user.
Note that in the silhouette plot, the labels on the lefthand side of the
chart show the units' names and the cluster number to which each unit is closer.
The silhouette plot is obtained from the silhouette()
function out from the cluster
package.
For a detailed description of the silhouette plot, its rationale, and its interpretation, see
Rousseeuw 1987.
Descriptive bycluster dotplots
The function also provides a Cleveland's dotplots that represent bycluster proportions. The
clustered units are grouped according to their cluster membership, the frequencies are summed, and
then expressed as percentages. The latter are represented by the dotplots, along with the average
percentage. The latter provides a frame of reference to understand which percentage is below,
above, or close to the average. The raw data on which the plots are based are stored in the
list returned by the function (see below).
Value
The function returns a list storing the following components

BR.similarity.matrix:
similarity matrix reporting the BR coefficients. 
Pvalue.matrix:
matrix reporting the permuted pvalues. 
classified.Pvalues.matrix:
matrix reporting the permuted pvalue classified as <0.05, <0.01, <0.001, or not significant. 
BR.distance_matrix:
distance matrix on which the hierarchical clustering is performed (returned if clustering is selected). 
avr.silh.width.by.n.of.clusters:
average silhouette width by number of clusters (returned if clustering is selected). 
partition.silh.data:
silhouette data for the selected partition (returned if clustering is selected). 
data.with.cluster.membership:
copy of the input data table with an additional column storing the cluster membership for each row (returned if clustering is selected). 
by.cluster.proportion:
table reporting the proportion of column categories across each cluster; rows sum to 100 percent (returned if clustering is selected).
References
Robinson, W. S. (1951). A Method for Chronologically Ordering Archaeological Deposits. In American Antiquity (Vol. 16, Issue 4, pp. 293–301). Cambridge University Press.
Robinson, W. S., & Brainerd, G. W. (1952). Robinson’s Coefficient of Agreement – A Rejoinder. In American Antiquity (Vol. 18, Issue 1, pp. 60–61). Cambridge University Press.
Rousseeuw P J. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. In Journal of Computational and Applied Mathematics 20, 5365.
DeBoer, W. R., Kintigh, K., & Rostoker, A. G. (1996). Ceramic Seriation and Site Reoccupation in Lowland South America. In Latin American Antiquity (Vol. 7, Issue 3, pp. 263–278). Cambridge University Press.
See Also
Examples
# build a toy dataset (subset of the 'Nelson' dataset out of the 'archdata' package )
mytable < structure(list(Biscuit = c(10, 17, 2, 10, 2, 1),
Type_I = c(2,2, 10, 40, 118, 107),
Type_II_Red = c(24, 64, 68, 91, 45, 3),
Type_II_Yellow = c(23, 90, 18, 20, 1, 0),
Type_II_Gray = c(34,76, 48, 15, 5, 0)),
row.names = c("1", "2", "3", "7", "8", "9"),
class = "data.frame")
# run the function and store the results in the 'result' object
result < brsim(mytable, clust=TRUE)
# same as above, but with an userdefined cluster partition
result < brsim(mytable, clust=TRUE, part=3)
# same as above, but rendering with a reordered heatmap
result < brsim(mytable, clust=TRUE, part=3, sort.map=TRUE)