R/other.R

In ManlyMix: Manly Mixture Modeling and Model-Based Clustering

Gmat <- function(p){

    M <- p * (p + 1) / 2

    G <- matrix(c(rep(NA, M * p^2)), ncol = M)

    n <- 1

    Z <- rep(0, M)

    for (a in 1:p){

        for (b in 1:p){

            Zn <- Z

            i1 <- max(a, b)            
            i2 <- min(a, b)

            ind <- M - (p - i2 + 1) * (p - i2 + 2) / 2 + i1 - i2 + 1

            Zn[ind] <- 1

            G[n,] <- Zn

            n <- n + 1
            
        }

    }

    return(G)

}
Manly.plot <- function(X, var1 = NULL, var2 = NULL, model = NULL, x.slice = 100, y.slice = 100, x.mar = 1, y.mar = 1, col = "lightgrey", ...){

	if(is.null(var2)){
		if(is.null(var1)){
			stop("Please specify which variable(s) to plot...\n")
		}
		else{
			cat("Univariate density plot is provided...\n")
			Manly.density(X, var = var1, model = model, x.slice = x.slice, x.mar = x.mar, col = col, ...)
		}
	}
	else{
		cat("Bivariate contour plot is provided...\n")
		Manly.contour(X, var1 = var1, var2 = var2, model = model, x.slice = x.slice, y.slice = y.slice, x.mar = x.mar, y.mar = y.mar, col = col, ...)
	}
}

Manly.density <- function(X, var = 1, model = NULL, x.slice = 100, x.mar = 1, col = "lightgrey", ...){	
	if(!is.matrix(X)){
		if(is.vector(X)){
			n <- length(X)
			X <- matrix(X, n, 1)
		}
	}
	n <- dim(X)[1]
	K <- length(model$tau)
	hist(X[,var], freq = FALSE, ...)
	x.seq <- seq(min(X[,var]) - x.mar, max(X[,var]) + x.mar, length = x.slice)	
	z <- rep(0, length(x.seq))

	for (i in 1:length(x.seq)){
		z[i] <- Manly.mix(x.seq[i], la = matrix(model$la[,var], K, 1), tau = model$tau, Mu = matrix(model$Mu[,var], K, 1), S = array(model$S[var,var,], dim = c(1,1,K)))
	}	
	points(x.seq, z, type = "l", col = col, ...)
	box()


}

Manly.contour <- function(X, var1 = 1, var2 = 2, model = NULL, x.slice = 100, y.slice = 100, x.mar = 1, y.mar = 1, col = "lightgrey", ...){
	n <- dim(X)[1]
	p <- dim(X)[2]
 	if(p != 2) warning("Only two variables are used for the contour plot...\n")

	if(n < 1) stop("Wrong number of observations n...\n")
	if(p < 1) stop("Wrong dimensionality p...\n")


	x.seq <- seq(min(X[,var1]) - x.mar, max(X[,var1]) + x.mar, length = x.slice)
	y.seq <- seq(min(X[,var2]) - y.mar, max(X[,var2]) + y.mar, length = y.slice)
	z <- matrix(0, length(x.seq), length(y.seq))

	for (i in 1:length(x.seq)){
		for (j in 1:length(y.seq)){
			z[i,j] <- Manly.mix(c(x.seq[i], y.seq[j]), la = model$la[,c(var1, var2)], tau = model$tau, Mu = model$Mu[,c(var1, var2)], S = model$S[c(var1, var2),c(var1, var2),])
		}
	}	

	contour(x.seq, y.seq, z, col = col, ...)
	points(X[,c(var1, var2)], col=model$id, cex = 0.5, pch = 19)

}

Manly.mix <- function(X, la = NULL, tau = NULL, Mu = NULL, S = NULL){


	if(is.null(dim(X))){
		n <- 1
		p <- dim(la)[2]

	}
	else{
		n <- dim(X)[1]
		p <- dim(X)[2]
		
	}
	K <- length(tau)

	if( is.null(Mu) || is.null(la) || is.null(tau) || is.null(S)) stop("Must provide the parameters to calculate mixture density from...\n")

	equal.K <- c(dim(la)[1], dim(Mu)[1], dim(S)[3])
	equal.p <- c(dim(la)[2], dim(Mu)[2], dim(S)[1], dim(S)[2])


	if(K < 1) stop("Wrong number of mixture components K...\n")
	if ((K != equal.K[1]) || (K != equal.K[2]) || (K != equal.K[3])) stop("Inconsistent number of mixture components K...\n")
	if ((p != equal.p[1]) || (p != equal.p[2]) || (p != equal.p[3]) || (p != equal.p[4])) stop("Inconsistent number of dimensionaltiy p...\n")



	x1 <- as.vector(t(X))
	la1 <- as.vector(t(la))
	Mu1 <- as.vector(t(Mu))
	S1 <- as.vector(S)


	misc_int <- c(p, n, K)
	mix1 <- rep(0, n)


	result <- .C("run_Manly_mix", x1 = as.double(x1), misc_int = as.integer(misc_int), tau = as.double(tau), Mu1 = as.double(Mu1), S1 = as.double(S1), la1 = as.double(la1), mix = as.double(mix1), PACKAGE = "ManlyMix")	
	
	return(result$mix)



}



Manly.var <- function(X, model = NULL, conf.CI = NULL){

	n <- dim(X)[1]
	p <- dim(X)[2]
	K <- length(model$tau)
	
	Inf.mat <- NULL

	tau <- model$tau
	gamma <- model$gamma
	Mu <- model$Mu
	S <- model$S
	la <- model$la

	x1 <- as.vector(t(X))


	misc_int <- c(p, n, K)
	
	Y <- array(NA, c(n, p, K))

	for(k in 1:K){

		y1 <- rep(0, n*p)
		result <- .C("run_Manly_transX", x1 = as.double(x1), misc_int = as.integer(misc_int), la1 = as.double(la[k,]), y1 = as.double(y1), PACKAGE = "ManlyMix")
		Y[,,k] <- matrix(result$y1, nrow = n, byrow = TRUE)

	}


	for (i in 1: n){

		grad1 <- NULL

		for (k in 1:(K-1)){

			grad1 <- c(grad1, gamma[i,k] / tau[k] - gamma[i,K] / tau[K])

		}


		grad2 <- NULL


		for (k in 1:K){
			grad2 <- c(grad2, gamma[i,k] * solve(S[,,k]) %*% (Y[i,,k] - Mu[k,]))

		}


		grad3 <- NULL

		I <- matrix(0,p,p)
		diag(I) <- rep(1,p)

		for (k in 1:K){

			grad3 <- c(grad3, t(Gmat(p)) %*% as.vector(gamma[i,k] / 2 * solve(S[,,k]) %*% ((Y[i,,k] - Mu[k,]) %*% t(Y[i,,k] - Mu[k,]) %*% solve(S[,,k]) - I)))


		}

		grad4 <- NULL


		for (k in 1:K){

			Dk <- matrix(0,p,p)
			for (j in 1:p){
				if(la[k,j] != 0){
					Dk[j,j] <- (1 + exp(la[k,j] * X[i,j]) * (X[i,j] * la[k,j] - 1)) / la[k,j]^2
				}
			}

			vect.temp <- -gamma[i,k] * Dk %*% solve(S[,,k]) %*% (Y[i,,k] - Mu[k,]) + gamma[i,k] * X[i,]
			index <- la[k,] == 0
			grad4 <- c(grad4, vect.temp[!index])


		}
		q <-  c(grad1, grad2, grad3, grad4)
	
		
		Inf.mat <- cbind(Inf.mat, q)

	}
	Inf.mat <- Inf.mat %*% t(Inf.mat)

	V <- solve(Inf.mat)

	if(!is.null(conf.CI)){

		if((conf.CI >1) || (conf.CI <=0)) stop("The confidence level has to be in between 0 and 1...\n")
		st.err <- sqrt(diag(V))
		Estimates <- c(model$tau[-K], as.vector(t(model$Mu)), unique(as.vector(model$S)), as.vector(t(model$la)))
		quantile <- (1- conf.CI)/2 +conf.CI
		
		Lower <- Estimates - qnorm(quantile) * st.err
		Upper <- Estimates + qnorm(quantile) * st.err
		CI <- cbind(Estimates, Lower, Upper)
	
	}
	else{
		CI <- NULL

	}

	return(list(V = V, CI = CI))
}




Manly.sim <- function(n, la, tau, Mu, S){

	if(n < 1)  stop("Wrong sample size n...\n")
	if(sum(tau) != 1)  stop("Mixing proportions should sum up to 1...\n")

    	K <- length(tau)
    	p <- dim(Mu)[2]

	equal.K <- c(dim(la)[1], dim(Mu)[1], dim(S)[3])
	equal.p <- c(dim(la)[2], dim(S)[1], dim(S)[2])

	if (K < 1) stop("Wrong number of mixture components K...\n")
	if ((K != equal.K[1]) || (K != equal.K[2]) || (K != equal.K[3])) stop("Inconsistent number of mixture components K...\n")
	if ((p != equal.p[1]) || (p != equal.p[2]) || (p != equal.p[3])) stop("Inconsistent number of dimensionality p...\n")

	Y <- NULL
	Y <- matrix(rnorm(n * p), n, p)
        Nk <- rep(1, K) + drop(rmultinom(1, n - K, tau))
   	
    	id <- NULL
    	for (k in 1:K) {id <- c(id, rep(k, Nk[k]))}
	for (k in 1:K) {
	   	eS <- eigen(S[,,k], symmetric = TRUE)
   		ev <- eS$values	
		temp <- eS$vectors %*% diag(sqrt(ev), p) %*% t(Y[id == k, ]) + Mu[k,]
        	Y[id == k, ] <- t(temp)

	}


	X <- NULL	
	for (k in 1:K){
		ind <- id == k
		Z <- sweep(Y[ind,], 2, STATS = la[k,], FUN = "*")
		Z <- log(Z + 1)
		Z <- sweep(Z, 2, STATS = la[k,], FUN = "/")
		X <- rbind(X, Z)
	}

	return(list(X = X, id = id))
}	













prop.M <- function(Z, tau1, tau2, mu1, mu2, S1, S2, la1, la2, sqr.S1, inv.S1, inv.S2, det.S1, det.S2){
	Yk <- sqr.S1 %*% Z + mu1
	index <- la1!=0

	X <- Yk
	if(sum(index)>0){
		X[index] <- log(la1[index] * Yk[index] + 1) / la1[index]
	}
	index2 <- la2!=0	
	Ykprime <- X
	if(sum(index2)>0){
		Ykprime[index2] <- (exp(la2[index2] * X[index2]) - 1) / la2[index2]
	}
		
	a <- 1 / 2 * t(Ykprime - mu2) %*% inv.S2 %*% (Ykprime - mu2) - t(la2) %*% X -
	1 / 2 * t(Yk - mu1) %*% inv.S1 %*% (Yk - mu1) + t(la1) %*% X
	b <- log(tau2 / tau1 * det.S2^(-1/2) / det.S1^(-1/2))
	re <- ifelse(a < b, 1, 0)
	return(re)
}


overlap.M2 <- function(N, tau, Mu, S, la){
	p <- dim(Mu)[2]
	simu <- matrix(rnorm(N * p), N, p)
	inv.S1 <- solve(S[,,1])
	inv.S2 <- solve(S[,,2])
	det.S1 <- det(S[,,1])
	det.S2 <- det(S[,,2])
	a.eig <- eigen(S[,,1])
	a.sqrt <- a.eig$vectors %*% diag(sqrt(a.eig$values)) %*% solve(a.eig$vectors)

	b.eig <- eigen(S[,,2])
	b.sqrt <- b.eig$vectors %*% diag(sqrt(b.eig$values)) %*% solve(b.eig$vectors)

	options(warn=-1)
	Omega12 <- sum(apply(simu, 1, prop.M, tau1 = tau[1], tau2 = tau[2], mu1 = Mu[1,], mu2 = 
	Mu[2,], S1 = S[,,1], S2 = S[,,2], la1 = la[1,], la2 = la[2,], sqr.S1 = 
	a.sqrt, inv.S1 = inv.S1, inv.S2 = inv.S2, det.S1 = det.S1, det.S2 = det.S2), na.rm = TRUE) / N

	Omega21 <- sum(apply(simu, 1, prop.M, tau1 = tau[2], tau2 = tau[1], mu1 = Mu[2,], mu2 = 
	Mu[1,], S1 = S[,,2], S2 = S[,,1], la1 = la[2,], la2 = la[1,], sqr.S1 = 
	b.sqrt, inv.S1 = inv.S2, inv.S2 = inv.S1, det.S1 = det.S2, det.S2 = det.S1), na.rm = TRUE) / N
	return(list(Omega12 = Omega12, Omega21 = Omega21))
}


Manly.overlap <- function(tau, Mu, S, la, N = 1000){

	p <- length(tau)
	K <- dim(Mu)[1]
	mat <- combn(1:K, 2)
	OmegaMap <- matrix(0, K, K)
	Components <- NULL
	Overlap <- rep(NA, dim(mat)[2])
	OverlapMap <- matrix(NA, 3, dim(mat)[2])

	for(i in 1:dim(mat)[2]){
		whichtwo <- mat[,i]
		result <- overlap.M2(N, tau = tau[whichtwo], Mu = Mu[whichtwo,], S = S[,,whichtwo], la = la[whichtwo,])
		OmegaMap[whichtwo[1], whichtwo[2]] <- result$Omega12
		OmegaMap[whichtwo[2], whichtwo[1]] <- result$Omega21
		Overlap[i] <- result$Omega12 + result$Omega21
		Components <- c(Components, paste("(", mat[1,i], ", ", mat[2,i], ")", sep=""))

		
	}
	OverlapMap <- data.frame(Components, Overlap)

	for(k in 1:K){
		OmegaMap[k, k] <- 1-sum(OmegaMap[k,])
	}
	

	BarOmega <- mean(Overlap)
	MaxOmega <- max(Overlap)
	return(list(OmegaMap = OmegaMap, OverlapMap = OverlapMap, BarOmega = BarOmega, MaxOmega = MaxOmega))

}






ClassAgree <- function(est.id, trueid){
	if (length(est.id) != length(est.id)){
		stop("Lengths of the estimated id and true id do not match...\n")
	}
	n <- length(est.id)
	A <- as.factor(est.id)
	B <- as.factor(trueid)
	K1 <- nlevels(A)
	K2 <- nlevels(B)
	if(K1 < K2){
		combination <- rep(0, K2)
	}
	else{
		combination <- rep(0, K1)
	}
	for (i in 1:K1){
		ind <- A == levels(A)[i]
		est.id[ind] <- i 
	}
	for (i in 1:K2) {
        	ind <- B == levels(B)[i]
       		trueid[ind] <- i 
    	}


	if (min(K1, K2) == 1) {

		ClassificationTable <- table(trueid, est.id)
		MisclassificationNum <- n - max(ClassificationTable)
		
	}
	else{
		
		Q <- .C("runProAgree", n = as.integer(n), K1 = as.integer(K1), 
        K2 = as.integer(K2), id1 = as.integer(est.id-1), id2 = as.integer(trueid-1), 
        maxPro = as.double(0), combination = as.integer(combination), PACKAGE = "ManlyMix")
		label <- Q$combination + 1

		if(K1 < K2){
			for(i in 1:n){
				est.id[i] <- label[est.id[i]]
			}
		}
		else{
			for(i in 1:n){
				trueid[i] <- label[trueid[i]]

			}
		}
		MisclassificationNum <- (1- Q$maxPro) * n
		ClassificationTable <- table(trueid, est.id)
	}
	
	return(list(ClassificationTable = ClassificationTable, MisclassificationNum = MisclassificationNum))
}