# R/seqivardur.R In TraMineR: Trajectory Miner: a Sequence Analysis Toolkit

#### Documented in print.seqivardurseqivardur

```seqivardur <- function(seqdata, type=1, with.missing=FALSE){

if (!inherits(seqdata,"stslist"))
stop(" [!] data is NOT a sequence object, see seqdef function to create one")
if (!type %in% c(1,2))
stop(" [!] type should be 1 or 2!")

dss <- suppressMessages(seqdss(seqdata, with.missing=with.missing))

lgth <- seqlength(seqdata, with.missing=with.missing)
dlgth <- suppressMessages(seqlength(dss, with.missing=with.missing))
sdist <- suppressMessages(seqistatd(seqdata, with.missing=with.missing))
nnvisit <- rowSums(sdist==0)

## the var function in R gives the unbiased variance with the n-1 denominator
## but we need the real variance here
realvar <- function(x) {
n <- sum(!is.na(x))
var <- 1/n*sum((x - mean(x,na.rm=TRUE))^2,na.rm=TRUE)
return(var)
}

dur <- seqdur(seqdata, with.missing=with.missing)

if (type == 1) {
ret <- apply(dur,1,realvar)
meand <- apply(dur,1,mean,na.rm=TRUE)
var.max <- (dlgth-1)*(1-meand)^2
} else if (type == 2) {
meand <- rowSums(dur, na.rm=TRUE)/(dlgth + nnvisit)
ddur <- (dur - as.vector(meand))^2
ret <- (rowSums(ddur, na.rm=TRUE) + nnvisit*meand^2) /  (dlgth + nnvisit)
alph <- alphabet(seqdata, with.missing=with.missing)
alph.size <- length(alph)
if (alph.size<2) maxnnv <- 0
else
maxnnv <- ifelse(dlgth==1, alph.size-1, alph.size-2)
meand.max <- meand*(dlgth + nnvisit)/(dlgth + maxnnv)
##print(meand.max)
var.max <- ((dlgth-1)*(1-meand.max)^2 + (lgth-dlgth+1-meand.max)^2 + maxnnv*meand.max^2) / (dlgth + maxnnv)
}

ret <- as.vector(ret)
attr(ret,"vmax") <- as.vector(var.max)
attr(ret,"meand") <- as.vector(meand)

class(ret) <- c("seqivardur", class(ret))
return(ret)
}

print.seqivardur <- function(x, stat='var', ...) {
## Conversion for printing without attributes
statlist <- c('mean','std','var','vmax','all')
if (any(!stat %in% statlist)) msg.stop.in("stat",statlist)

dvar <- as.vector(x)
std <- sqrt(as.vector(x))
meand <- as.vector(attr(x,'meand'))
vmax <- as.vector(attr(x,'vmax'))

x <- cbind(meand,std,dvar,vmax)
a <- which(stat=='all')
if (length(a)>0)  stat <- stat[-a]
if (length(stat) >0) x <- x[,which(statlist %in% stat)]

NextMethod("print", ...)
}
```

## 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.