使用作为 TraMineR
的一部分的生物燃料数据集:
library(TraMineR)
data(biofam)
lab <- c("P","L","M","LM","C","LC","LMC","D")
biofam.seq <- seqdef(biofam[,10:25], states=lab)
head(biofam.seq)
Sequence
1167 P-P-P-P-P-P-P-P-P-LM-LMC-LMC-LMC-LMC-LMC-LMC
514 P-L-L-L-L-L-L-L-L-L-L-LM-LMC-LMC-LMC-LMC
1013 P-P-P-P-P-P-P-L-L-L-L-L-LM-LMC-LMC-LMC
275 P-P-P-P-P-L-L-L-L-L-L-L-L-L-L-L
2580 P-P-P-P-P-L-L-L-L-L-L-L-L-LMC-LMC-LMC
773 P-P-P-P-P-P-P-P-P-P-P-P-P-P-P-P
我可以执行聚类分析:
library(cluster)
couts <- seqsubm(biofam.seq, method = "TRATE")
biofam.om <- seqdist(biofam.seq, method = "OM", indel = 3, sm = couts)
clusterward <- agnes(biofam.om, diss = TRUE, method = "ward")
cluster3 <- cutree(clusterward, k = 3)
cluster3 <- factor(cluster3, labels = c("Type 1", "Type 2", "Type 3"))
但是,在此过程中,来自biofam.seq的唯一ID已被数字1到N的列表所取代:
head(cluster3, 10)
[1] Type 1 Type 2 Type 2 Type 2 Type 2 Type 3 Type 3 Type 2 Type 1
[10] Type 2
Levels: Type 1 Type 2 Type 3
现在,我想知道每个簇中哪些序列,以便我可以应用其他函数来获得每个簇内的平均长度,熵,子序列,不相似性等 . 我需要做的是:
-
将旧ID映射到新ID
-
将每个簇中的序列插入到单独的序列对象中
-
在每个新序列对象上运行我想要的统计信息
如何在上面的列表中完成2和3?
2 回答
我想这会回答你的问题 . 我使用了我在这里找到的代码http://www.bristol.ac.uk/cmm/software/support/workshops/materials/solutions-to-r.pdf来创建
biofam.seq
,因为你所建议的都没有为我工作 .首先,我使用
split
为每个集群创建索引列表,然后我在lapply
循环中使用它来创建biofam.seq
的子序列列表:最后,您可以使用
lapply
或sapply
对每个子序列运行分析其中
FUN
可以是您通常在单个序列上运行的任何函数 .例如,可以简单地获得第一簇的状态序列对象