1

文字列として表される異なる「空間使用のシーケンス」間の類似性を測定するために、シーケンス分析メソッドを使用しています。以下は、2 つのシーケンスに対して 3 つのクラス (A: 都市、B: 農業、C: 山) を使用した理論上の例です。

                t­1,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**

シーケンス間の類似性を測定するために使用する距離尺度は、ハミング距離です (つまり、シーケンスを同一視するためにシーケンス内の文字を置換する必要がある頻度を測定します。上記の例では、順番に4文字を置換する必要があります)。シーケンスを同一視します)。ハミング距離の計算後に得られた距離行列 (可能なすべてのシーケンスのペアの距離または非類似度を与える) に基づいて、Ward (ward.D2) のクラスタリング方法を使用して樹状図が作成されました。

ここで、関連するクラスターを識別するために、クラスターの堅牢性の適切な尺度も含めたいと思います。このために、ブートストラップ値を計算するためのいくつかの方法を含む pvclust を使用しようとしましたが、距離測定の数に制限されていました。リリースされていないバージョンの pvclust を使用して、適切な距離測定 (つまり、ハミング距離) を実装しようとし、ブートストラップ ツリーを作成しようとしました。スクリプトは機能していますが、結果が正しくありません。1000 の 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 の未リリース バージョンを使用しています。これにより、独自の距離法 (この場合はハミング) を使用できます。誰かがこの問題を解決する方法を知っていますか?

4

1 に答える 1

1

の目的は、ケースではなく変数(または属性)pvclustをクラスター化することです。これが、意味をなさない結果が得られる理由です。あなたが試すことができます

data(iris)
res <- pvclust(iris[, 1:4])
plot(res)

ケースのクラスタリングの安定性をテストするには、clusterbootfrom packageを使用できますfpc。ここで私の答えを参照してください:ツリー/デンドログラムの信頼性を測定する (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 つのクラスターがあるため (構造上)、悪い結果が得られます。

于 2015-03-26T22:22:11.753 に答える