R/vg_mixture.R

In mixture: Mixture Models for Clustering and Classification

# MAIN VGPCM FUNCTION
vgpcm <- function(data=NULL,  G=1:3, mnames=NULL, # main inputs with mnames being the Model Name 
				start=2, label=NULL, # starting inputs , start = 0: random soft, start = 2, random hard. start = 3 mkmeans. 
				veo=FALSE, da=c(1.0), # veo (variables exceed observations), da is deterministic annealing  
				nmax=1000, atol=1e-8, mtol=1e-8, mmax=10, burn=5, # convergence settings for matrix and loglik
				pprogress=FALSE, pwarning=FALSE, stochastic = FALSE, latent_method = "standard")  # progress settings 
{

	# Do some sanity checks. 
	if (is.null(data)) stop('Hey, we need some data, please! data is null')
	if (!is.matrix(data)) stop('The data needs to be in matrix form')
	if (!is.numeric(data)) stop('The data is required to be numeric')
	if (nrow(data) == 1) stop('nrow(data) is equal to 1')
	if (ncol(data) == 1) stop('ncol(data) is equal to 1; This function currently only works with multivariate data p > 1')
	# check for full NAs vector. as mixture version 1.6+ can handle missing data.  
	apply(data,1,checkNA)

	# some more sanity checks. 
	if (is.null(G)) stop('G is NULL')
	G = as.integer(ceiling(G))
	if (!is.integer(G)) stop('G is not a integer')
	if ( any(G < 1)) stop('G is not a positive integer')
	n <- dim(data)[1]
	p <- dim(data)[2]

	row_tags <- c()
	# grab row_tags if na 
	if(any(is.na(data))){
		for(i in 1:n)
		{
			if(any(is.na(data[i,])))
			{
				row_tags <- append(i,row_tags)
			}
		}
	}

	# IF MODEL NAME NOT SPECIFIED DO EM ON ALL POSSIBLE MODELS 
	if (is.null(mnames) )  mnames = c("EII", "VII", "EEI", "VEI", "EVI", "VVI", "EEE", "EEV", "VEV", "VVV", "EVE", "VVE", "VEE", "EVV")

	# CREATE MULTIDIMENSIONAL ARRAY CONTAINING THE FOLLOWING:
	# G x MODEL_LENGTH x INFO (loglik, npar, BIC)
	info_BIC <- c() # BIC for each model 
	info_loglik <- c() # logliklihoods for each model
	info_npar <- c() # gives number of parameters vector
	info_model_lexicon <- c() # gives header tag for the full model, useful for print statements and summary
	info_model <- c() # stores full model data 

	# SET START OBJECT
	# handeling z_ig matrices. 
	if ( is.matrix(start))
	{
		if(length(G) != 1){
			stop("Initialization z matrix should only be for a single G")
		}
		if(dim(start)[1] != dim(data)[1])
		{
			stop("Initialization z matrix should have the same number of rows as data")
		}
		if(dim(start)[2] != G)
		{
			stop("Initialization z matrix should have G number of columns")
		}
		startobject <- "matrix"
	}
	else if ( start == 0){
		startobject <- "random_soft" # check for validity of start object. 
	} 
	else if ( start == 1){
		startobject <- "random_hard" 
	}
	else if ( start == 2){
	  if ( any(is.na(data)) ) { stop("You cannot use kmeans on missing values, try soft initialization first, then after a kmeans start.") }
		startobject <- "kmeans"
	} 
	else if (start > 2) {
		startobject <- "multi"
	}
	else{
		stop("start setting is not valid")
	}

	# handle latent methods. 
	if( !( latent_method %in% c("random", "standard") ) )
	{
		stop("latent method is invalid, can only be either random, or standard")
	}

	# deterministic annealing sanity checks
	if(any(is.na(da))){ stop("deterministic annealing should contain no NAs") } # check NAs
	else if(!is.vector(da)) { stop("deterministic annealing should be a vector")} # check vector 
	else if(!is.numeric(da)) { stop("deterministic annealing should be numeric")} # check numeric 

	# burn setting check 
	if(!is.numeric(burn)){ stop("burn-in setting has to be a number")}
	if(round(burn) <= 0) { stop("burn-in setting has to be a positive round number")}
	if(burn != round(burn)){
		if(!pwarning){
			warning("Warning: rounding burn-in setting number")
		}
		burn <- round(burn) 
	}

	# label sanity checks 
	if(!is.null(label)) { # if there is no null
		if(!is.vector(label)){ stop("label has to be a vector")} #check that its a vector. 
		else if(any(is.na(label))) { stop("No NAs allowed in label vector.")} # check for NAs
		else if(any(label < 0)) { stop("labels have to be greater or equal to 0")} # make sure they are greater than 0
		else if(any(label > G)) { stop("labels cant be bigger than G groups")}  # make sure they are greater than G
		else if(all(label == 0)) { stop("There is no need a vector with all 0!, just don't pass in anything")} # check for stupidity. 
		else if(length(label) != n) { stop("label vector has to be the same size as the number of observtions. ")} # check length. 
		else if(length(G) != 1) { stop("Labels can only be inputted in for a single G") } # check for single input G. 
	}

	# go through all G's 
	for(G_i in G)
	{
		# go through all models
		for(model_name in mnames)
		{
			if(pprogress) {	cat("Running VGPCM Model:",model_name,"G=",G_i,"\n")}
			# calculate number of parameters once. 
			number_of_params <- npar.model.skew(model_name,G=G_i,p=p,family_name = "VG")

			check_veo <- TRUE

			if(number_of_params > n){
				if(veo){
					if(!pwarning){
					warning("Model: ",model_name , " G: " , G_i ," ","Number of Parameters exceed number of observations.\n")
					}
				}
				else{
					check_veo <- FALSE
				}
			}

			# check to see if number of parameters exceed observations. 
			if(check_veo){

				model_id <- model_to_id(model_name)
				model_id_stochastic_check <- model_id

				if(G_i > 1){
					# set up intialization matrix depending on choice. 
					in_zigs <- switch(startobject,
					"random_soft"=z_ig_random_soft(n,G_i),
					"random_hard"=z_ig_random_hard(n,G_i),
					"kmeans"=z_ig_kmeans(data,G_i),
					"matrix"=start,
					"multi"=z_ig_random_soft(n, G_i))


					# handle labels within z_ig matrix.  
					if(!is.null(label)){
						model_id_stochastic_check <- 2

						if(stochastic)
						{
							stop('Under current version, you cannot have semi-supervision and stochastic EM')
						}
						
						# start observation count
						i <- 1
						# go through the entire label system 
						for(label_i in label){
							# if its not zero, replace the one entry in z_ig with just 0s and 1 one. 
							if(label_i != 0){
								# create vector of zeros and put a 1 based on the label i 
								classif_vector <- rep(0,G_i) # rep 0. 
								classif_vector[label_i] <- 5 # put one 
								in_zigs[i,] <- classif_vector # replace
							}
							i <- i + 1
						}
					}
				}
				else{
					in_zigs <- as.matrix(rep(1.0,n))
				}

				
				# Handle multiple starts 
				if(startobject == "multi" && G_i > 1){

					# run first model with settings. 
					model_results_i <- main_loop_vg(X = t(data),
															G = G_i, in_zigs = in_zigs, model_id = model_id_stochastic_check,
															model_type = model_id + stochastic*20, in_nmax = nmax, in_l_tol = atol,
															in_m_iter_max = mmax, in_m_tol = mtol, anneals = da,
															t_burn = burn, latent_step = latent_method)

					# handle errors just in case. 
					if(!is.null(model_results_i$Error)){
						best_loglik <- -Inf 
						best_z_start <- z_ig_random_soft(n, G_i)
					}
					else{
						best_loglik <- tail(model_results_i$logliks, 1) 
						best_z_start <- model_results_i$zigs
					}

					# attempt multiple starts. 
					for(i in 2:start){

						# run rest model with settings. 
						model_results_i <- main_loop_vg(X = t(data),
																G = G_i, in_zigs = in_zigs, model_id = model_id_stochastic_check,
																model_type = model_id + stochastic*20, in_nmax = nmax, in_l_tol = atol,
																in_m_iter_max = mmax, in_m_tol = mtol, anneals = da,
																t_burn = burn, latent_step = latent_method)

						# handle errors. 
						if(!is.null(model_results_i$Error)){
							current_loglik <- -Inf 
							current_z_start <- z_ig_random_soft(n, G_i)
						}
						else{
							current_loglik <- tail(model_results_i$logliks, 1) 
							current_z_start <- model_results_i$zigs
						}

						# check if the start was better. 
						if(current_loglik > best_loglik){
							best_loglik <- current_loglik
							best_z_start <- current_z_start
						}

					}

					# give a good start. 
					in_zigs <- as.matrix(best_z_start)

				}
			  # RUN MODEL 
			  model_results_i <- main_loop_vg(X = t(data),
			                                  G = G_i, in_zigs = in_zigs, model_id = model_id_stochastic_check,
			                                  model_type = model_id + stochastic*20, in_nmax = nmax, in_l_tol = atol,
			                                  in_m_iter_max = mmax, in_m_tol = mtol, anneals = da,
			                                  t_burn = burn, latent_step = latent_method)
			  status <- "Failed Aitken's Convergence Criterion"
			  if (nmax > length(model_results_i$logliks)) {
			    status <- "Converged according to Aitken's Convergence Criterion"
			  }

				if(pwarning && !is.null(model_results_i$Error) ){
					cat(paste(model_name,"| G =",G_i,":",model_results_i$Error),"\n")
				}
			
			  # AQUIRE STATUS IF ANY 
			  model_results_i$status <- status
			  if(is.null(model_results_i$Error)){
			    info_loglik <- append(tail(model_results_i$logliks,1),info_loglik) # first one is loglik 
			    info_npar <- append(number_of_params,info_npar) # number of paramaters. 
			    info_BIC <- append(2*info_loglik[1] - log(n)*info_npar[1],info_BIC) # append BIC value 
			    info_model_lexicon <- append(paste("Model:", model_name, "G: ",G_i),info_model_lexicon) # lexicon for summary
			    model_results_i$sigs <- convert_matrices(model_results_i,G_i,p) # convert matrices into proper form
			    info_model <- append(list(model_results_i),info_model) # store all model objects. 
			  }
			  else {
					info_loglik <- append(NA,info_loglik) # first one is loglik 
					info_npar <- append(number_of_params,info_npar) # number of paramaters. 
					info_BIC <- append(NA,info_BIC) # append BIC value 
					info_model_lexicon <- append(paste("Model:", model_name, "G: ",G_i),info_model_lexicon) # lexicon for summary
					info_model <- append(list(NA),info_model) # store all model objects. 
				}
			}
		}

	}

	if(start > 2){
		startobject <- paste(start, "random soft initializations.")
	} 

	info_matrix <- list(startobject=startobject, # gives starting object information. 
						info_loglik=info_loglik,info_npar=info_npar,info_BIC=info_BIC, # logliks, params, BICs, 
						lexicon=info_model_lexicon,model_objs=info_model,Gs=G,row_tags=(row_tags+1)) # lexicon (tags), and model_objs 
	if(length(info_matrix$model_obj) < 1){                                       # plus 1 for row_tags due to c++
		stop("No models estimated")
	}
	info_matrix$best_model <- vg_get_best_model(info_matrix)
	info_matrix$BIC <- construct_BIC_3D(info_matrix)
  info_matrix$map <- MAP(info_matrix$best_model$model_obj[[1]]$zigs)
  info_matrix$best_model$map <- info_matrix$map 


    gpar= list()
	for (k in 1:info_matrix$best_model$G ) {
		gpar[[k]] = list()		
		gpar[[k]]$mu       = info_matrix$best_model$model_obj[[1]]$mus[[k]]
		gpar[[k]]$sigma    = info_matrix$best_model$model_obj[[1]]$sigs[[k]]
		gpar[[k]]$alpha    = info_matrix$best_model$model_obj[[1]]$alphas[[k]]
		gpar[[k]]$invSigma = try(solve(info_matrix$best_model$model_obj[[1]]$sigs[[k]]))
		gpar[[k]]$logdet   = info_matrix$best_model$model_obj[[1]]$log_dets[k]

	}

	gpar$pi = info_matrix$best_model$model_obj[[1]]$pi_gs

	info_matrix$gpar <- gpar
	info_matrix$z <- info_matrix$best_model$model_obj[[1]]$zigs



	class(info_matrix) <- "vgpcm"

	return(info_matrix)

}

# PRINT SUMMARY AND PLOT STATEMENTS
print.vgpcm<-function(x, ...){
  # split strings and parse
	splitted_strings <- strsplit(x$best_model$model_type," ")[[1]]
  # print to line parsed strings. 
  cat("The model chosen by applying the BIC criteria has", trimws(splitted_strings[5]), "component(s) and the", trimws(splitted_strings[2]), "covariance structure\n using", x$startobject, "\n"  )
}
# just prints the compare BIC matrix. 
summary.vgpcm<- function(object, ...){
    cat("BIC for each model, number of components (rows), and covariance structure (columns).\n")
	print(object$BIC[,,3])
}

# plots a line graph of BICs. 
plot.vgpcm<- function(x, ...) {
	bicl = x$BIC[,,3]
	# you need to wrap this plot up in a print statement otherwise it doesnt work. 
	print(levelplot(bicl,
          col.regions = colorRampPalette(c("black","brown","red","orange","gold","yellow","white"))(prod(dim(bicl)) + 10),
          xlab = "G",ylab = "Covariance Type"))
	trellis.focus("legend", side="right", clipp.off=TRUE, highlight=FALSE)
	trellis.unfocus()	
}



# Print output for the gpcm_best class. 
print.vgpcm_best <-function(x, ...){

	splitted_strings  <- strsplit(x$model_type," ")[[1]]
	cov_string <- splitted_strings[2]
	component_string <- splitted_strings[5]

	cat("============================\n")
    cat("Best VG Model According To BIC \n")
	cat("============================\n")
	cat("Status: ", strsplit(x$status," ")[[1]][1], "\n")
	cat("Covariance Model Type: ",cov_string,"\n")
	cat("Number of Components: ",component_string,"\n")
	cat("Initalization: ",x$startobject,"\n")
	cat("BIC: ", x$BIC, "\n")
	cat("============================\n")
	
}


# Function: calculates the best model from a full list of models and their objects. 
# Get best model according to BIC 
vg_get_best_model <- function(gpcm_model)
{	
  # sanity checks for input 
  if(!is.list(gpcm_model)) { stop("Error: Input is not a gpcm_model") }
  if(!is.list(gpcm_model$model_objs)) { stop("Error: model_objs missing... input is not a gpcm_model ") }
  
	# replace infinity and NaNs with NA 
	gpcm_model$info_BIC[!is.finite(gpcm_model$info_BIC)] <- NA
	
	if(all(is.na(gpcm_model$info_BIC))){
		stop("error: no models estimated.")
	}

	# get best_model index. 
	bm_index <- try(match(try(max(gpcm_model$info_BIC,na.rm = TRUE)), gpcm_model$info_BIC))

  # get G 
  lexicon_best <- strsplit(gpcm_model$lexicon[bm_index]," ")[[1]]
  G_best <- as.numeric(lexicon_best[5])
  Cov_type <- lexicon_best[2]
  
  status <- gpcm_model$model_objs[[bm_index]]$status
  # construct best model 
  best_model <- list(model_type=gpcm_model$lexicon[bm_index], # what the call is. 
                     model_obj=gpcm_model$model_objs[bm_index], # return model objects for the best model. 
                     BIC=gpcm_model$info_BIC[bm_index], loglik=gpcm_model$info_loglik[bm_index], # get best BIC and loglik. 
                     nparam=gpcm_model$info_npar[bm_index],startobject=gpcm_model$startobject, G=G_best,cov_type=Cov_type ,# number of parameters and G. 
                     status=status) # convergence status
  
  



  best_model$map <- MAP(best_model$model_obj[[1]]$zigs)
  
  class(best_model) <- "vgpcm_best"
  
  return(best_model)
  
}