Как вычислить различия между последовательностями, когда последовательности содержат пробелы? - PullRequest
0 голосов
/ 10 декабря 2018

Я хочу кластеризовать последовательности с оптимальным соответствием с TraMineR::seqdist() из данных, которые содержат пропуски, то есть последовательности, содержащие пропуски.

library(TraMineR)
data(ex1)
sum(is.na(ex1))

# [1] 38

sq <- seqdef(ex1[1:13])
sq

#    Sequence                 
# s1 *-*-*-A-A-A-A-A-A-A-A-A-A
# s2 D-D-D-B-B-B-B-B-B-B      
# s3 *-D-D-D-D-D-D-D-D-D-D    
# s4 A-A-*-*-B-B-B-B-D-D      
# s5 A-*-A-A-A-A-*-A-A-A      
# s6 *-*-*-C-C-C-C-C-C-C      
# s7 *-*-*-*-*-*-*-*-*-*-*-*-*

sm <- seqsubm(sq, method='TRATE')
round(sm,digits=3)

#      A-> B->   C-> D->
# A->   0 2.000   2 2.000
# B->   2 0.000   2 1.823
# C->   2 2.000   0 2.000
# D->   2 1.823   2 0.000

Когда я запускаю seqdist()

dist.om <- seqdist(sq, method="OM", indel=1, sm=sm)

Я получаю

Error: 'with.missing' must be TRUE when 'seqdata' or 'refseq' contains missing values

, но когда я устанавливаю опцию with.missing=TRUE, я получаю

 [>] including missing values as an additional state
 [>] 7 sequences with 5 distinct states
 [>] checking 'sm' (one value for each state, triangle inequality)
Error:  [!] size of substitution cost matrix must be 5x5

Итак, как мы можем вычислить различия между последовательностями, используя seqdist() свывод seqsubm() правильный путь, когда данные содержат пропуски, то есть последовательности содержат пропуски?

Примечание: Я не совсем уверен, имеет ли это смысл вообще.Пока я просто исключаю наблюдения с пропусками, но из-за своих данных я теряю много наблюдений по этому поводу.Поэтому было бы полезно узнать, есть ли такая опция.

1 Ответ

0 голосов
/ 10 декабря 2018

Существуют разные стратегии для вычисления расстояний, когда у вас есть пробелы.

1) Первое решение состоит в том, чтобы рассматривать отсутствующее состояние как дополнительное состояние.Это то, что seqdist делает, когда вы устанавливаете with.missing=TRUE.В этом случае матрица sm должна содержать затраты на замену любого состояния отсутствующим состоянием.Используя seqsubm, вам просто нужно указать with.missing=TRUE и для этой функции.По умолчанию затраты замещения для замены «отсутствующими» устанавливаются как фиксированное значение miss.cost (по умолчанию 2).

sm <- seqsubm(sq, method='TRATE', with.missing=TRUE)
round(sm,digits=3)

#     A->   B-> C->   D-> *->
# A->   0 2.000   2 2.000   2
# B->   2 0.000   2 1.823   2
# C->   2 2.000   0 2.000   2
# D->   2 1.823   2 0.000   2
# *->   2 2.000   2 2.000   0

Чтобы получить затраты на замещение для «отсутствующих» на основе вероятностей перехода

sm <- seqsubm(sq, method='TRATE', with.missing=TRUE, miss.cost.fixed=FALSE)
round(sm,digits=3)

#       A->   B->   C->   D->   *->
# A-> 0.000 2.000 2.000 2.000 1.703
# B-> 2.000 0.000 2.000 1.823 1.957
# C-> 2.000 2.000 0.000 2.000 1.957
# D-> 2.000 1.823 2.000 0.000 1.957
# *-> 1.703 1.957 1.957 1.957 0.000

Используя последний sm, мы получаем расстояния между последовательностями

dist.om <- seqdist(sq, method="OM", indel=1, sm=sm, with.missing=TRUE)
round(dist.om, digits=2)

#       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]
# [1,]  0.00 22.87 21.91 18.41  6.41 17.00 17.03
# [2,] 22.87  0.00 13.76 11.56 19.91 19.87 22.57
# [3,] 21.91 13.76  0.00 14.25 18.96 18.91 21.57
# [4,] 18.41 11.56 14.25  0.00 13.70 15.70 18.14
# [5,]  6.41 19.91 18.96 13.70  0.00 15.70 16.62
# [6,] 17.00 19.87 18.91 15.70 15.70  0.00 16.70
# [7,] 17.03 22.57 21.57 18.14 16.62 16.70  0.00

Конечно, последовательности тогда будут находиться близко друг от друга только потому, что они имеют много пропущенных состояний (*).Поэтому вы можете захотеть сохранить только те последовательности, в которых, например, пропущено менее 10% элементов.

2) Вторым решением является удаление пробелов, что вы делаете в seqdef.(Однако обратите внимание, что это меняет выравнивание.)

## Here, we drop seq 7 that contains only missing values

sq <- seqdef(ex1[-7,1:13], left='DEL', gaps='DEL')
sq

#    Sequence           
# s1 A-A-A-A-A-A-A-A-A-A
# s2 D-D-D-B-B-B-B-B-B-B
# s3 D-D-D-D-D-D-D-D-D-D
# s4 A-A-B-B-B-B-D-D    
# s5 A-A-A-A-A-A-A-A    
# s6 C-C-C-C-C-C-C  

sm <- seqsubm(sq, method='TRATE')
round(sm,digits=3)

#       A->   B-> C->   D->
# A-> 0.000 1.944   2 2.000
# B-> 1.944 0.000   2 1.823
# C-> 2.000 2.000   0 2.000
# D-> 2.000 1.823   2 0.000

dist.om <- seqdist(sq, method="OM", indel=1, sm=sm)
round(dist.om, digits=2)

#       [,1]  [,2]  [,3]  [,4]  [,5] [,6]
# [1,]  0.00 19.61 20.00 13.78  2.00   17
# [2,] 19.61  0.00 12.76  9.59 17.61   17
# [3,] 20.00 12.76  0.00 13.29 18.00   17
# [4,] 13.78  9.59 13.29  0.00 11.78   15
# [5,]  2.00 17.61 18.00 11.78  0.00   15
# [6,] 17.00 17.00 17.00 15.00 15.00    0
...