R/as.seqtree.R

In WeightedCluster: Clustering of Weighted Data

###########################
## disstree main function
###########################
as.seqtree <- function(object, seqdata, diss, weighted=TRUE,...){
	UseMethod("as.seqtree")
}
as.seqtree.twins <- function(object, seqdata, diss, weighted=TRUE, ncluster, ...) {
	return(as.seqtree.hclust(object, seqdata=seqdata, diss=diss, weighted=weighted, ncluster=ncluster,...))
}
as.seqtree.hclust <- function(object, seqdata, diss, weighted=TRUE, ncluster, ...) {
	pred <- data.frame(Split2=factor(cutree(object, 2)))
	for(p in 3:ncluster){
		pred[, paste("Split", p, sep="")] <- factor(cutree(object, p))
	}
	object <- pred
	return(as.seqtree.default(object, seqdata=seqdata, diss=diss, weighted=weighted, ...))
}

as.seqtree.default <- function(object, seqdata, diss, weighted=TRUE, ...) {
	predictor <- object
	ncluster <- ncol(object)+1
	if(ncluster<2){
		stop(" [!] ncluster should be bigger than 2")
	}
	if (inherits(diss, "dist")) {
		diss <- as.matrix(diss)
 	}
	## Model matrix from forumla
	nobs= nrow(diss)

	## Allow integer weights for replicates
	if(weighted & !is.null(attr(seqdata,"weights") )){
		weights <- attr(seqdata,"weights")
	}
	else {
		weights <- as.double(rep(1,nobs))
	}
	pop <- sum(weights)
#	TraMineR:::.localstuffDissTree$DTNnodeCounter <- as.integer(1)

	vardis <- dissvar(diss, weights=weights)

	as.seqtreeDTNBuildNode <- function(ind, vardis, depth, current) {
		node <- TraMineRInternalNodeInit(ind=ind, vardis=vardis, depth=depth, dmat=diss, weights=weights)
		node$info$splitschedule <- depth
		SCtot <- vardis*node$info$n
		SCres <- SCtot
		if(current>ncluster ||length(ind)==1){
			return(node)
		}
		## print(SCtot)
		#varnames <- colnames(pred)

		for (p in current:ncluster) {
			clust <- predictor[ind, p-1]
			SplitLabels <- unique(clust)
			if(length(SplitLabels)==2){
				#print(p)
				#print(ncluster)
				bestSpl <- list()
				bestSpl$variable <- clust == SplitLabels[1]
				# print(ind)
				# print(as.integer(ind[bestSpl$variable]))
				# print(diss)
				lSCres <- TraMineRInternalWeightedInertiaDist(diss, as.integer(nrow(diss)),
					as.integer(FALSE), as.integer(ind[bestSpl$variable]), as.double(weights),
					as.integer(FALSE))
				rSCres <- TraMineRInternalWeightedInertiaDist(diss, as.integer(nrow(diss)),
					as.integer(FALSE), as.integer(ind[!bestSpl$variable]), as.double(weights),
					as.integer(FALSE))
				info <- list(
						lpop=sum(weights[ind[bestSpl$variable]]),
						rpop=sum(weights[ind[!bestSpl$variable]]),
						SCres=lSCres+rSCres
				)
				info$lvar=lSCres/info$lpop
				info$rvar=rSCres/info$rpop
				bestSpl$spl <- TraMineRInternalSplitInit(p-1, index = 1:2,
					 prob = c(info$lpop, info$rpop)/node$info$n, info = info, labels=SplitLabels)
				SCres <- bestSpl$spl$info$SCres
				node$split <- bestSpl$spl
				node$split$info$R2 <- 1-(SCres/SCtot)
				node$surrogates <- bestSpl$sur

				#print(bestSpl)
				left <- as.seqtreeDTNBuildNode(ind=ind[bestSpl$variable], vardis=bestSpl$spl$info$lvar, depth=p, current=p+1)
				right <- as.seqtreeDTNBuildNode(ind=ind[!bestSpl$variable], vardis=bestSpl$spl$info$rvar, depth=p, current=p+1)
				node$kids <- list(left, right)
				## We have found the split, so leave the loop
				return(node)
			}
		}
		## Maximum depth reached
		return(node)
	}

	root <- as.seqtreeDTNBuildNode(ind=1:nobs, vardis=vardis, depth=1, current=2)
	#print(root)
	tree <- list()
	tree$fitted <- data.frame(disstreeleaf(root))
	names(tree$fitted) <- "(fitted)"
	tree$info <- list(method="disstree", n=pop, parameters= list(minSize=1, maxdepth=ncluster, R=0, pval=1), object=seqdata, weight.permutation="diss")
	if(!weighted) {
		tree$info$adjustment <- dissassoc(diss, tree$fitted[,1], R=0, weights=NULL)
	}
	else {
		tree$info$adjustment <- dissassoc(diss, tree$fitted[,1], R=0, weights=weights, weight.permutation="diss")
	}
	tree$data <- predictor
	tree$terms <- NULL
	tree$weights <- weights
	##tree <- party(root, data=predictor, fitted =fitted, terms = terms(formula.call),  info = info)
	tree$root <- root

	class(tree) <- c("seqtree", "disstree", class(tree))
	return(tree)
}