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

#### Defines functions print.TraMineRPermuthist.TraMineRPermutsummary.TraMineRPermutTraMineR.permutationweightTraMineR.permutation

```TraMineR.permutation <- function(data, R, statistic, ...) {
if (is.vector(data)) {
n <- length(data)
}
else {
n<- nrow(data)
}
perms <- list()
perms\$R <- sum(R)
perms\$t0 <- statistic(data,1:n,...)
ntest <- length(perms\$t0)
if (perms\$R>1) {
perms\$t <- matrix(NA, perms\$R,length(perms\$t0))
perms\$t[1,] <- perms\$t0
for(i in 2:(perms\$R)){
perms\$t[i,] <- statistic(data,sample(n), ...)
}
perms\$pval <- numeric(length(perms\$t0))

for(i in 1:length(perms\$pval)){
perms\$pval[i] <- (sum(perms\$t[, i]>=perms\$t0[i]))/perms\$R
}
}
else {
perms\$t <- NA
perms\$pval <- as.numeric(rep(NA, ntest))
}
class(perms) <- "TraMineRPermut"
return(perms)
}

TraMineR.permutationweight <- function(data, R, statistic, samplesize, sampleprob, t0, ...) {
if (is.vector(data)) {
n <- length(data)
}
else {
n<- nrow(data)
}
perms <- list()
perms\$R <- sum(R)
allind <- 1:n
perms\$t0 <- t0
ntest <- length(perms\$t0)
# spop <- 1:samplesize
#print(data.frame(allind,sampleprob))
if (perms\$R>1) {
#perms\$t0boot <- matrix(NA, perms\$R,length(perms\$t0))
perms\$t <- matrix(NA, perms\$R,length(perms\$t0))
perms\$t[1,] <- perms\$t0
#perms\$t0boot[1,] <- perms\$t0
for(i in 2:perms\$R){
pop <- sort.int(sample(allind, size=samplesize, replace=TRUE, prob=sampleprob), method="quick")
pop2 <- sample(allind, size=samplesize, replace=TRUE, prob=sampleprob)
#bt <- statistic(pop,pop2, ...)
perms\$t[i,] <- statistic(pop,pop2, ...)
#perms\$t0boot[i,] <- bt[-(1:ntest)]
}
perms\$pval <- numeric(length(perms\$t0))

for(i in 1:length(perms\$pval)){
perms\$pval[i] <- sum(perms\$t[, i]>=perms\$t0[i])/perms\$R
}
}
else {
perms\$t <- NA
perms\$pval <- as.numeric(rep(NA, ntest))
}
class(perms) <- c("TraMineRPermut")
return(perms)
}

summary.TraMineRPermut <- function(object, ...){
perms <- object
resume <- data.frame(t0 = perms\$t0, "p-value"=perms\$pval)
rownames(resume) <- paste("t[,",1:length(perms\$t0),"]")
return(resume)
}
hist.TraMineRPermut <- function(x, index=1, breaks="FD",main=paste("Distribution of test statistic number ", index), xlab=index, pvalue.limit=NULL, freq=FALSE, type="density", txtcex=1, ...) {
ntest <- length(x\$t0)
if(index >ntest){
stop("index should be in range 1, ", ntest)
}
if (x\$R==0) {
stop("Cannot plot permutation test (none detected)")
}
if(type=="hist") {
hist(x\$t[, index], main=main, xlab=xlab, breaks=breaks, freq=freq, ...)
if (!is.null(pvalue.limit)) {
testbootorder <- order(x\$t[, index], decreasing=TRUE)
abline(v=x\$t[testbootorder[round(pvalue.limit*x\$R)], index], col="blue")
}
abline(v=x\$t0[index], col="red")
}
else if(type=="density"){
dens <- density(x\$t[,index])
##hist(x\$t[,index], xlim = range(dens\$x), xlab = "x", ylab = "density", freq = FALSE)
plot(dens, main=main, xlab=xlab, ...)
nullhypregion <- dens\$x>=x\$t0[index]
polygon(x=c(min(dens\$x[nullhypregion]),dens\$x[nullhypregion]),y=c(0,dens\$y[nullhypregion]),  col="red")
if(txtcex>0) {
text(dens\$x[nullhypregion][1], dens\$y[nullhypregion][1],labels=paste("  Statistic:",  format(x\$t0[index], digits =3), "(P-value:",x\$pval[index],")"), adj=c(0,0), cex= txtcex)
}
}

}

print.TraMineRPermut <- function(x,...){
cat("\nPermutation test\n\n")
print(summary(x),...)
cat("P-value approximated using", x\$R, "permutations\n")
}
```

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