引导自己构建的函数pvclust不起作用
我使用序列分析的方法来度量不同的“空间使用序列”之间的相似性，用字符串表示。这里有一个理论例子，分为三个类(A:城市，B:农业，C:山区)，分别是两个序列：
t1,t2,........tx Individual 1: A A A B B B C C Individual 2: A B B B A A C C 0 1 1 0 1 1 0 0 = \*\*4\*\*
我们用来度量序列间相似性的距离度量是hamming距离(即测量序列中字符需要替换多少次才能使序列相等，在上面的例子中，4字符需要替换才能等号)。根据计算汉明距离后得到的距离矩阵(给出每对可能序列的距离或不同)，利用Ward聚类方法(ward.D2)生成树状图。
现在，我还想介绍一种很好的集群健壮性度量，以便识别相关的集群。为此，我试图使用pvclust，它包含了几种计算引导值的方法，不过，它仅限于一些距离度量。使用未发布版本的pvclust，我尝试实现正确的距离度量(即hamming距离)，并尝试创建一个引导树。脚本正在运行，但结果不正确。应用在我的数据集上，使用1,000的nboot，"bp“值接近0，所有其他值"au”、"se.au“、"se.bp”、"v“、"c”、"pchi“都是0，这表明集群是人工制品。
这里我提供了一个示例脚本：
数据涉及非常同质的模拟序列(例如，继续使用1种特定状态)，因此每个集群都应该是显着的。我将引导的次数限制在10，以限制计算时间。
####################################################################
####Create the sequences####
dfr = data.frame()
a = list(dfr)
b = list(dfr)
c = list(dfr)
d = list(dfr)
data = list(dfr)
for (i in c(1:10)){
set.seed(i)
a[[i]] <- sample(c(rep('A',10),rep('B', 90)))
b[[i]] <- sample(c(rep('B',10),rep('A', 90)))
c[[i]] <- sample(c(rep('C',10),rep('D', 90)))
d[[i]] <- sample(c(rep('D',10),rep('C', 90)))
}
a = as.data.frame(a, header = FALSE)
b = as.data.frame(b, header = FALSE)
c = as.data.frame(c, header = FALSE)
d = as.data.frame(d, header = FALSE)
colnames(a) <- paste(rep('seq_urban'),rep(1:10), sep ='')
colnames(b) <- paste(rep('seq_agric'),rep(1:10), sep ='')
colnames(c) <- paste(rep('seq_mount'),rep(1:10), sep ='')
colnames(d) <- paste(rep('seq_sea'),rep(1:10), sep ='')
data = rbind(t(a),t(b),t(c),t(d))
#####################################################################
####Analysis####
## install packages if necessary
#install.packages(c("TraMineR", "devtools"))
library(TraMineR)
library(devtools)
source_url("https://www.dropbox.com/s/9znkgks1nuttlxy/pvclust.R?dl=0") # url to my dropbox for unreleased pvclust package
source_url("https://www.dropbox.com/s/8p6n5dlzjxmd6jj/pvclust-internal.R?dl=0") # url to my dropbox for unreleased pvclust package
dev.new()
par( mfrow = c(1,2))
## Color definitions and alphabet/labels/scodes for sequence definition
palet <- c(rgb(230, 26, 26, max = 255), rgb(230, 178, 77, max = 255), "blue", "deepskyblue2") # color palet used for the states
s.alphabet <- c("A", "B", "C", "D") # the alphabet of the sequence object
s.labels <- c("country-side", "urban", "sea", "mountains") # the labels of the sequence object
s.scodes <- c( "A", "U", "S", "M") # the states of the sequence object
## Sequence definition
seq_ <- seqdef(data, # data
1:100, # columns corresponding to the sequence data
id = rownames(data), # id of the sequences
alphabet = s.alphabet, states = s.scodes, labels = s.labels,
xtstep = 6,
cpal = palet) # color palet
##Substitution matrix used to calculate the hamming distance
Autocor <- seqsubm(seq_, method = "TRATE", with.missing = FALSE)
# Function with the hamming distance (i.e. counts how often a character needs to be substituted to equate two sequences to each other. Result is a distance matrix giving the distances for each pair of sequences)
hamming <- function(x,...) {
res <- seqdist(x, method = "HAM",sm = Autocor)
res <- as.dist(res)
attr(res, "method") <- "hamming"
return(res)
}
## Perform the bootstrapping using the distance method "hamming"
result <- pvclust(seq_, method.dist = hamming, nboot = 10, method.hclust = "ward")
result$hclust$labels <- rownames(test[,1])
plot(result)为了进行这种分析，我使用R包pvclust的未发布版本，它允许使用您自己的距离方法(在本例中是: hamming)。有人知道如何解决这个问题吗？
回答 1
Stack Overflow用户
发布于 2015-03-26 22:22:11
pvclust的目标是对变量(或属性)进行聚类，而不是在情况下。这就是为什么你有一些没有意义的结果。你可以试试
data(iris)
res <- pvclust(iris[, 1:4])
plot(res)要测试案例集群的稳定性，可以从包clusterboot中使用fpc。见我在这里的答案：Measuring reliability of tree/dendrogram (Traminer)
在您的示例中，您可以使用：
library(fpc)
ham <- seqdist(seq_, method="HAM",sm = Autocor)
cf2 <- clusterboot(as.dist(ham), clustermethod=disthclustCBI, k=4, cut="number", method="ward.D")例如，使用k=10，结果会很糟糕，因为您的数据实际上有4个集群(按构造计算)。
https://stackoverflow.com/questions/29284468
复制相似问题
腾讯云开发者
