pvclust: Calculating P-values for Hierchical Clustering

Description

calculates $p$-values for hierarchical clustering via multiscale bootstrap resampling. Hierarchical clustering is done for given data and $p$-values are computed for each of the clusters.

Usage

pvclust(data, method.hclust="average",
        method.dist="correlation", use.cor="pairwise.complete.obs",
        nboot=1000, r=seq(.5,1.4,by=.1), store=FALSE, weight=FALSE)

parPvclust(cl=NULL, data, method.hclust="average",
           method.dist="correlation", use.cor="pairwise.complete.obs",
           nboot=1000, r=seq(.5,1.4,by=.1), store=FALSE, weight=FALSE,
           init.rand=TRUE, seed=NULL, iseed=NULL)

Arguments

data

numeric data matrix or data frame.

method.hclust

the agglomerative method used in hierarchical clustering. This should be (an abbreviation of) one of "average", "ward", "single", "complete", "mcquitty", "median"

method.dist

the distance measure to be used. This should be (an abbreviation of) one of "correlation", "uncentered", "abscor" or those which are allowed for method argument in dist functi

use.cor

character string which specifies the method for computing correlation with data including missing values. This should be (an abbreviation of) one of "all.obs", "complete.obs" or "pairwise.complete.obs"

nboot

the number of bootstrap replications. The default is 1000.

numeric vector which specifies the relative sample sizes of bootstrap replications. For original sample size $n$ and bootstrap sample size $n'$, this is defined as $r=n'/n$.

store

locical. If store=TRUE, all bootstrap replications are stored in the output object. The default is FALSE.

a cluster object created by package parallel or snow. If NULL, use the registered default cluster.

weight

logical. If weight=TRUE, resampling is made by weight vector instead of index vector. Useful for large r value (r>10). Currently, available only for distance "correlation" and "absc

init.rand

logical. If init.rand=TRUE, random number generators are initialized. Use iseed argument to achieve reproducible results.

seed

integer vector of random seeds. It should have the same length as cl. This argument is duplicated and will be unavailable in the future. Consider using iseed instead.

iseed

an integer to be passed to clusterSetRNGStream when init.rand=TRUE.

Value

hclusthierarchical clustering for original data generated by function hclust. See hclust for details.
edgesdata frame object which contains $p$-values and supporting informations such as standard errors.
countdata frame object which contains primitive information about the result of multiscale bootstrap resampling.
msfitlist whose elements are results of curve fitting for multiscale bootstrap resampling, of class msfit. See msfit for details.
nbootnumeric vector of number of bootstrap replications.
rnumeric vector of the relative sample size for bootstrap replications.
storelist contains bootstrap replications if store=TRUE was given for function pvclust or parPvclust.

Details

Function pvclust conducts multiscale bootstrap resampling to calculate $p$-values for each cluster in the result of hierarchical clustering. parPvclust is the parallel version of this procedure which depends on package parallel for parallel computation. For data expressed as $(n \times p)$ matrix or data frame, we assume that the data is $n$ observations of $p$ objects, which are to be clustered. The $i$'th row vector corresponds to the $i$'th observation of these objects and the $j$'th column vector corresponds to a sample of $j$'th object with size $n$. There are several methods to measure the dissimilarities between objects. For data matrix $X={x_{ij}}$, "correlation" method takes $$1 - \frac{ \sum_{i=1}^n (x_{ij} - \bar{x}_j) (x_{ik} - \bar{x}_k) } { \sqrt{\sum_{i=1}^n (x_{ij} - \bar{x}_j)^2} \sqrt{\sum_{i=1}^n (x_{ik} - \bar{x}_k)^2} }$$ for dissimilarity between $j$'th and $k$'th object, where $\bar{x}_j = \frac{1}{n} \sum_{i=1}^n x_{ij} \mbox{and} \bar{x}_k = \frac{1}{n} \sum_{i=1}^n x_{ik}$. "uncentered" takes uncentered sample correlation $$1 - \frac{ \sum_{i=1}^n x_{ij} x_{ik} } { \sqrt{\sum_{i=1}^n x_{ij}^2} \sqrt{\sum_{i=1}^n x_{ik}^2} }$$ and "abscor" takes the absolute value of sample correlation $$1 - \ \Biggl| \frac{ \sum_{i=1}^n (x_{ij} - \bar{x}_j) (x_{ik} - \bar{x}_k) } { \sqrt{\sum_{i=1}^n (x_{ij} - \bar{x}_j)^2} \sqrt{\sum_{i=1}^n (x_{ik} - \bar{x}_k)^2} } \Biggl|.$$

References

Shimodaira, H. (2004) "Approximately unbiased tests of regions using multistep-multiscale bootstrap resampling", Annals of Statistics, 32, 2616-2641. Shimodaira, H. (2002) "An approximately unbiased test of phylogenetic tree selection", Systematic Biology, 51, 492-508. Suzuki, R. and Shimodaira, H. (2004) "An application of multiscale bootstrap resampling to hierarchical clustering of microarray data: How accurate are these clusters?", The Fifteenth International Conference on Genome Informatics 2004, P034. http://www.is.titech.ac.jp/~shimo/prog/pvclust/

Examples

Run this code

## using Boston data in package MASS
library(MASS)
data(Boston)

## multiscale bootstrap resampling
boston.pv <- pvclust(Boston, nboot=100)

## CAUTION: nboot=100 may be too small for actual use.
##          We suggest nboot=1000 or larger.
##          plot/print functions will be useful for diagnostics.

## plot dendrogram with p-values
plot(boston.pv)

ask.bak <- par()$ask
par(ask=TRUE)

## highlight clusters with high au p-values
pvrect(boston.pv)

## print the result of multiscale bootstrap resampling
print(boston.pv, digits=3)

## plot diagnostic for curve fitting
msplot(boston.pv, edges=c(2,4,6,7))

par(ask=ask.bak)

## Print clusters with high p-values
boston.pp <- pvpick(boston.pv)
boston.pp

## parallel computation
library(parallel)
cl <- makeCluster(2, type = "PSOCK")

## parallel version of pvclust
boston.pv <- parPvclust(cl, Boston, nboot=1000)
stopCluster(cl)

Run the code above in your browser using DataLab