# R/seqdist-OMstran.R In TraMineR: Trajectory Miner: a Sequence Analysis Toolkit

#### Defines functions OMstran

```# Should only be used through seqdist()

OMstran <- function(seqdata, indel, sm, full.matrix, transindel, otto, previous,
add.column, with.missing, weighted, refseq, norm) {

alph <- alphabet(seqdata, with.missing=with.missing)
void <- attr(seqdata,"void")
dimnames(sm) <- list(alph, alph)
if(length(indel)==1){
indel <- rep(indel, length(alph))
}
names(indel) <- alph
if(transindel=="prob"){
tr <- seqtrate(seqdata, with.missing=with.missing, weighted=weighted)
dimnames(tr) <- list(alph, alph)
}
sl <- seqlength(seqdata)
seqdata <- as.matrix(seqdata)
maxcol <- ncol(seqdata)
newseqdata <-matrix("", nrow=nrow(seqdata), ncol=maxcol)
sep <- "@@@@TraMineRSep@@@@"

mypastefunc <- function(i){
minmax <- function(i){
return(max(1, min(i, maxcol)))
}
tostate <- seqdata[,minmax(i+1)]
# voids can only be at end of sequences
# when to state is void we set it as the previous state when add.column
tostate[tostate==void] <- seqdata[tostate==void,minmax(i)]
}
ret <- paste( seqdata[, minmax(i)], tostate, sep=sep)
# set transition as NA when from state is void
ret[seqdata[,minmax(i)]==void] <- NA
}else{
# set transition as NA when to state is void
ret[tostate==void] <- NA
}
if(previous){
ret.na <- is.na(ret)
ret <- paste(seqdata[, minmax(i-1)], ret, sep=sep)
ret[ret.na] <- NA
## ret[seqdata[,minmax(i-1)]==void] <- NA  ## cannot happen for !ret.na
}
return(ret)
}

for (i in 1:(maxcol)) {
newseqdata[, i] <- mypastefunc(i)
}
newseqdata <- newseqdata[, -maxcol]
if(previous){
newseqdata <- newseqdata[, -1]
}
}
## transition to void states have been set as NA  (to/from nr will be considered as a distinct transition)

newalph <- unique(as.character(newseqdata))
newalph <- newalph[!is.na(newalph)] ## deleting NAs
alphabet_size <- length(newalph)
suppressMessages(newseqdata <- seqdef(newseqdata, cpal=rep("blue", alphabet_size),
left='DEL', gaps='DEL', right='DEL',
id = rownames(seqdata)) )
transweight <- 1 - otto
if(previous){
transweight <- transweight/2
}
indelrate <- max(indel)/max(sm)
transweight <- transweight* indelrate
indel <- indelrate*indel/max(indel)

sm <- sm/max(sm)
## =========================================
## Building the new substitution cost matrix
## =========================================

## Build substitution matrix and new alphabet
alphabet <- alphabet(newseqdata)
alphabet_size <- length(alphabet)
msg("Creating", alphabet_size, "distinct transition states")
## Recomputing the subsitution matrix
indels <- numeric(alphabet_size)
names(indels) <- alphabet
newsm <- matrix(0, nrow=alphabet_size, ncol=alphabet_size)
stateindel <- indel* otto
## indel costs change according to the previous parameters
if(!previous){
for(i in 1:alphabet_size){
statesi <- strsplit(alphabet[i], sep)[[1]]
indels[i] <- stateindel[statesi[1]]
if(statesi[1]!= statesi[2]){
if(transindel=="constant"){
rawtransindel <- 1
}else if(transindel=="prob"){
rawtransindel <- (1-tr[statesi[1], statesi[2]])
}else if(transindel=="subcost"){
rawtransindel <- (sm[statesi[1], statesi[2]])
}
indels[i] <- indels[i] +  transweight*rawtransindel
}
}

}else{
for(i in 1:alphabet_size){
statesi <- strsplit(alphabet[i], sep)[[1]]
rawtransindel <- 0
if(transindel=="constant"){
rawtransindel <- (statesi[1]!= statesi[2])+(statesi[2]!= statesi[3])
}else if(transindel=="prob"){
rawtransindel <- (2-tr[statesi[1], statesi[2]]-tr[statesi[2], statesi[3]])
}else if(transindel=="subcost"){
rawtransindel <- (sm[statesi[1], statesi[2]]+sm[statesi[2], statesi[3]])
}
indels[i] <- stateindel[statesi[2]] +  transweight*rawtransindel
}

}
## Just fix the state we are comparing according to 'previous'
compare_state <- ifelse(previous, 2, 1)
for (i in 1:(alphabet_size-1)) {
statesi <- strsplit(alphabet[i], sep)[[1]]
for (j in (i+1):alphabet_size) {
statesj <- strsplit(alphabet[j], sep)[[1]]
cost <- sm[statesi[compare_state], statesj[compare_state]]
if(transindel %in% c("constant", "prob", "subcost")){
cost <- otto*cost + (indels[alphabet[i]] + indels[alphabet[j]] -stateindel[statesi[compare_state]] -stateindel[statesj[compare_state]])
}
newsm[i, j] <- cost
newsm[j, i] <- cost
}
}
suppressMessages(return(seqdist(newseqdata, method = "OM", indel = indels, sm = newsm, full.matrix = full.matrix,
weighted = weighted, refseq = refseq, norm = norm)))
}
```

## Try the TraMineR package in your browser

Any scripts or data that you put into this service are public.

TraMineR documentation built on May 29, 2024, 5 a.m.