# NOT RUN {
## Parallel setup
library(parallel)
length(cl <- makeCluster(detectCores()))
## script to create lung73.pvclust and lung73.sb
## multiscale bootstrap resampling of hierarchical clustering
library(pvclust)
data(lung)
### default pvclust scales
lung.pvclust <- pvclust(lung, nboot=10000, parallel=cl)
lung.sb <- sbfit(lung.pvclust,cluster=cl) # model fitting
### wider range of scales than pvclust default
sa <- 9^seq(-1,1,length=13)
lung73.pvclust <- pvclust(lung,r=1/sa,nboot=10000,parallel=cl)
lung73.sb <- sbfit(lung73.pvclust,cluster=cl) # model fitting
# }
# NOT RUN {
## replace si/au/bp entries in pvclust object
library(pvclust)
data(lung73) # loading the previously computed bootstrap
### the original pvclust result
plot(lung.pvclust, print.pv = c("si", "au", "bp"), cex=0.5, cex.pv=0.5)
pvrect(lung.pvclust, pv="si") # (defualt pvclust uses pv="au")
### default pvclust scales with p-values of k=2
lung.k2 <- sbpvclust(lung73.pvclust,lung73.sb, k=2)
plot(lung.k2, print.pv = c("si", "au", "bp"), cex=0.5, cex.pv=0.5)
pvrect(lung.k2, pv="si")
### wider scales with p-values of k=3 (default of scaleboot)
lung73.k3 <- sbpvclust(lung73.pvclust,lung73.sb)
plot(lung73.k3, print.pv = c("si", "au", "bp"), cex=0.5, cex.pv=0.5)
pvrect(lung73.k3, pv="si")
## diagnostics of fitting
### diagnose edges 61,...,69
lung73.sb[61:69] # print fitting details
plot(lung73.sb[61:69]) # plot curve fitting
summary(lung73.sb[61:69]) # print raw(=bp)/si/au p-values
### diagnose edge 67
lung73.sb[[67]] # print fitting
plot(lung73.sb[[67]],legend="topleft") # plot curve fitting
summary(lung73.sb[[67]]) # print au p-values
# }
Run the code above in your browser using DataLab