# R/mcen.r In mcen: Multivariate Cluster Elastic Net

#' Provides initial estimates for the mcen functionF
#'@param x the n x p design matrix
#'@param y the n x y matrix of responses
#'@param family type of likelihood used two options "mgaussian" or "mbinomial"
#'@param delta sparsity tuning parameter
#'@param gamma_y tuning parameter for clustering responses
#'@param intercept whether an intercept should be included in the model
#'@return matrix of coefficients
#'@importFrom glmnet glmnet
#'@importFrom stats coefficients
mcen.init <- function(x,y,family="mgaussian",delta=NULL,gamma_y=1,intercept=FALSE){
if(length(delta) > 1 | length(gamma_y) > 1){
stop('initialization only for single delta, gamma_y')
}
r <- dim(y)[2]
n <- dim(y)[1]
p <- dim(x)[2]
if(family=="mgaussian"){
init_family = "gaussian"
intercept_b=intercept
}
if(family=="mbinomial"){
init_family="binomial"
intercept_b=TRUE
x=x[,-1]
## probably need to force the intercept here to?
}
init_model <- apply(y,2,glmnet,family=init_family,x=x,alpha=delta/(delta+gamma_y),lambda=(delta+gamma_y)/(2*n),intercept=intercept_b)
if(intercept_b){
do.call(cbind,lapply(init_model,coefficients))
} else{
do.call(cbind,lapply(init_model,coefficients))[-1,]
}
}

#' Fits an MCEN model
#'@param x Matrix of predictors.
#'@param y Matrix of responses.
#'@param family Type of likelihood used two options "mgaussian" or "mbinomial".
#'@param ky Clusters for response.
#'@param delta L1 penalty.
#'@param gamma_y Penalty for with y clusters difference in predicted values.
#'@param ndelta Number of delta parameters.
#'@param delta.min.ratio Ratio between smallest and largest delta.
#'@param eps Convergence criteria.
#'@param scale_x Whether x matrix should be scaled, default is True.
#'@param scale_y Whether y matrix should be scaled, default is True.
#'@param clusterMethod K-means function used kmeans or kmeanspp.
#'@param clusterStartNum Number of random starting points for clustering.
#'@param clusterIterations Number of iterations for cluster convergence.
#'@param cluster_y An a priori definition of clusters. If clusters are provided they will remain fixed and are not estimated. Objective function is then convex.
#'@param max_iter Maximum number of iterations for coefficient estimates.
#'@param init_beta Clustering step requires an initial estimate, default is to use elastic net solution.
#'@param n.cores Number of cores used for calculation default is 1.
#'@return returns a MCEN object
#' \item{beta}{List of the coefficient estimates.}
#' \item{delta}{Value of delta.}
#' \item{gamma_y}{Value of gamma_y.}
#' \item{ky}{Value of ky.}
#' \item{y_clusters}{List of the clusters of y.}
#'@examples
#' x <- matrix(rnorm(400),ncol=4)
#' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4)
#' y <- x%*%beta + rnorm(400)
#' mcen_fit <- mcen(x,y,ky=2,delta=1)
#'@references Price, B.S. and Sherwood, B. (2018). A Cluster Elastic Net for Multivariate Regression. arXiv preprint arXiv:1707.03530. \url{http://arxiv-export-lb.library.cornell.edu/abs/1707.03530}.
#'@export
#'@importFrom stats cov

mcen <- function(x,y,family="mgaussian",ky=NULL,delta=NULL,gamma_y=1,ndelta=25,delta.min.ratio = NULL,
eps=.00001,scale_x=TRUE, scale_y=TRUE, clusterMethod="kmeans",clusterStartNum=30,clusterIterations=10,cluster_y=NULL,
max_iter=10,init_beta=NULL, n.cores=1){
if(is.null(cluster_y)==FALSE){
cluster_y_num <- length(unique(cluster_y))
if(is.null(ky)==FALSE){
if(cluster_y_num !=ky){
ky <- cluster_y_num
warning("ky does not match number of clusters in cluster_y, cluster_y then overrides ky")
}
}
}
if(scale_x){
x <- scale(x)
sigma_x <- attributes(x)$'scaled:scale' mean_x <- attributes(x)$'scaled:center'
}
if(scale_y && family!="mbinomial"){
y <- scale(y)
mean_y <- attributes(y)$'scaled:center' sigma_y <- attributes(y)$'scaled:scale'
}
p=dim(x)[2]
r=dim(y)[2]
n=dim(x)[1]
if(is.null(delta.min.ratio)==TRUE){
delta.min.ratio=ifelse(n < p, 0.01, 1e-04)
}
if(family=="mgaussian"){
Cxy=cov(x,y)*(n-1)
Cxx=cov(x)*(n-1)
}
if(is.null(delta)){
if(family=="mgaussian"){
delta_max = 2*max(abs(Cxy))+1
}
if(family=="mbinomial"){
## Need to add bens rule here
phat<-apply(y,2,mean)
zw<-log(phat/(1-phat))*phat*(1-phat)+y-phat
delta_max=2*max(abs(apply(zw,2,function(m){t(x)%*%zw})))+1
}
# plus one is just to make sure we get a full sparse matrix returned
delta_min = delta.min.ratio*delta_max
delta = exp(seq(log(delta_max),log(delta_min),length.out=ndelta))
}

ndelta <- length(delta)

delta_iter <- 1
beta_list <- list()
cluster_list <- list()
init_beta <- NULL
if(family=="mbinomial"){
x=cbind(rep(1,n),x)
}
for(i in 1:ndelta){
if(i == 1 || sum(init_beta)==0){
#first step get initial estimates, however we don't want to get stuck at initial estimates of all zeros for coefficients
if(length(delta)>1){
init_beta <- mcen.init(x,y,family,delta[i],gamma_y)
}else{
init_beta<-mcen.init(x,y,family,delta,gamma_y)
}
}
if(family=="mgaussian"){
beta_cluster <- mcen_workhorse(init_beta,delta[i],Cxx,Cxy,family,ky=ky,gamma_y=gamma_y,eps=eps,
clusterMethod=clusterMethod,clusterIterations=clusterIterations,
clusterStartNum=clusterStartNum,cluster_y=cluster_y, max_iter=max_iter, x=x)
}
if(family=="mbinomial"){

if(length(delta)>1){
beta_cluster<-mcen_bin_workhorse(init_beta,delta=delta[i],x=x,y=y,family="mbinomial",ky=ky,gamma_y=gamma_y,eps=eps,
clusterMethod=clusterMethod,clusterIterations=clusterIterations,
clusterStartNum=clusterStartNum,cluster_y=cluster_y, max_iter=max_iter)
}else{
beta_cluster<-mcen_bin_workhorse(init_beta,delta=delta,x=x,y=y,family="mbinomial",ky=ky,gamma_y=gamma_y,eps=eps,
clusterMethod=clusterMethod,clusterIterations=clusterIterations,
clusterStartNum=clusterStartNum,cluster_y=cluster_y, max_iter=max_iter)
}
}

init_beta <- beta_cluster$beta if(family=="mgaussian"){ beta_list[[i]] <- beta_adjust(init_beta,sigma_x,sigma_y,mean_x,mean_y) } if(family=="mbinomial"){ beta_list[[i]]<-beta_adjust_bin(init_beta,sigma_x) } cluster_list[[i]] <- beta_cluster$clusters
}
return_val <- list(beta=beta_list,delta=delta,gamma_y=gamma_y,ky=ky,y_clusters=cluster_list)
if(family=="mgaussian"){
class(return_val) <- c("mgauss_mcen", "mcen")
}
if(family=="mbinomial"){
class(return_val)<-c("mbinom_mcen","mcen")
}
return_val
}

#' Estimates the clusters and provides the coefficients for an mcen object
#'@param beta The initial value of the coefficients
#'@param delta The sparsity (L1) tuning parameter
#'@param xx Matrix of transpose of x times x.
#'@param xy Matrix of transpose of x times y.
#'@param family Type of likelihood used two options "mgaussian" or "mbinomial"
#'@param ky Number of clusters for the response
#'@param gamma_y Penalty for the y clusters difference in predicted values
#'@param eps Convergence criteria
#'@param clusterMethod Which clustering method was used, currently support kmeans or kmeanspp
#'@param clusterIterations Number of iterations for cluster convergence
#'@param clusterStartNum Number of random starting points for clustering
#'@param cluster_y An a priori definition of clusters. If clusters are provided they will remain fixed and are not estimated. Objective function is then convex.
#'@param max_iter The maximum number of iterations for estimating the coefficients
#'@param x The design matrix
#'@importFrom methods as
#'@importFrom Matrix sparseMatrix
#'@useDynLib mcen CDU
mcen_workhorse<-function(beta,delta=NULL,xx,xy,family="mgaussian",ky=NULL,gamma_y=.5,eps=.00001,
clusterMethod="kmeans",clusterIterations=100,clusterStartNum=30,cluster_y=NULL,max_iter=10, x=x){
#really need to think about how to do this function if x and y are not scaled and there is an intercept. This maybe important for the glm case
beta_dim <- dim(beta)
p <- beta_dim[1]
r <- beta_dim[2]
if(sum(beta==0)==p*r){
#if initial estimate is totally sparse then mcen update is not needed. Philosophically, I'm ok with this because mcen idea is
# to group coefficient estimates, if all coefficients are zero there is nothing to group and nothing for mcen to do.
return_val <- list(beta = as(beta,"sparseMatrix"), clusters=rep(0,r))
return( return_val)
} else{
fitted_vals <- matrix(x%*%beta,ncol=r)
iter <- 0
fitting_done <- FALSE
while(!fitting_done){
#print(paste("fitting done is ", fitting_done, " and iter is ", iter))
num_potential_centers <- ncol( unique(fitted_vals,MARGIN=2))
if(num_potential_centers < ky){
fitting_done <- TRUE
warning(paste0("Algorithm stopped early for gamma_y ", gamma_y, ", delta ", delta, " and ky ", ky, " because number of unqiue fitted values is less than ", ky, " Common for large values of deltas that cause spares solution and little variability in fitted values.\n"))
clusters_used <- rep(0,r)
y_clusters <- NULL
} else if(is.null(cluster_y)==FALSE){
clusters_used <- cluster_y
} else{
y_clusters <- cluster(fitted_vals,ky,clusterMethod,clusterIterations,clusterStartNum)
clusters_used <- y_clusters$cluster } if(iter > 0 && fitting_done==FALSE){ if(is.null(cluster_y)){ fitting_done=SetEq(y_clusters2$cluster,y_clusters$cluster) #fitting_done=SetEq(y_clusters2$centers,y_clusters$centers) } else{ fitting_done <- TRUE } } if(fitting_done == FALSE){ beta <-matrix(.C("CDU",beta=as.double(as.vector(beta)),xx=as.double(as.vector(xx)), xy=as.double(as.vector(xy)), y_clusters=as.double(clusters_used),#y_clusters$cluster),
delta=as.double(delta), gamma_y=as.double(gamma_y), eps=as.double(eps),miter=as.integer(max_iter), r=as.integer(r), p=as.integer(p), beta0=as.double(as.vector(beta+1)), set=as.double(rep(0,r)), First=as.double(0),
Second=as.double(0),mine=as.double(rep(0,p*r)))$beta, nrow=p,ncol=r,byrow=FALSE) } if(iter == max_iter){ fitting_done <- TRUE } iter <- iter + 1 if(is.null(cluster_y)==TRUE){ y_clusters2 <- y_clusters } fitted_vals=matrix(x%*%beta,ncol=r) } if(is.null(cluster_y)==FALSE){ fitting_done <- TRUE #only one iteration of fitting done if clusters are prespecified } } return_val <- list(beta = as(beta,"sparseMatrix"), clusters=clusters_used) return( return_val) } #' Adjusts the value of the coefficients to account for the scaling of x and y. #'@param beta The estiamte of beta with scaled data. #'@param sigma_x Sample tandard deviations of the original predictors. #'@param sigma_y Sample standard deviations of the orignal responses. #'@param mean_x Sample means of the original predictors . #'@param mean_y Sample means of the original responses. #'@return Returns the adjusted coefficients #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> beta_adjust <- function(beta,sigma_x,sigma_y,mean_x,mean_y){ beta <- beta*sigma_x for(i in 1:dim(beta)[2]){ beta[,i] <- beta[,i]*sigma_y[i] } intercept <- mean_y - mean_x %*% beta beta <- rbind(intercept,beta) beta } #' Adjusts the value of the binomial coefficients to account for the scaling of x. #'@param beta The estiamte of beta with scaled data. #'@param sigma_x Sample tandard deviations of the original predictors. #'@return Returns the adjusted coefficients #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> beta_adjust_bin<-function(beta,sigma_x){ beta_out<-beta beta_out[-1,]-beta[-1,]*sigma_x return(beta_out) } #' SetEq test set equivalence of two clustering sets #'@param set1 is the cluster assignments of the previous iteration #'@param set2 is the cluster assignments of the current clusters #'@return Returns a logical saying if the two clusterings are equal #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> SetEq <-function(set1,set2){ s1<-vector("list",max(set1)) for(i in 1:max(set1)){ s1[[i]]=which(set1==i) } Sets=sapply(s1,function(z){ s2=vector("list",max(set2)) for(i in 1:max(set2)){ s2[[i]]=which(set2==i) } sum(sapply(s2,function(w){identical(z,w)})) }) sum(Sets)==max(set1) } #' Wrapper function for different clustering methods #'@param x data to be clustered. Clustering will be done on the columns. #'@param cNum number of cluster centers #'@param clusterMethod "kmean" for kmeans function, "kmeanspp" for kcca implementation of kmeans++ #'@param clusterIterations number of maximum iterations for clustering #'@param clusterStartNum random number of starting points used #'@return Returns cluster assignments #'@importFrom flexclust kcca #'@importFrom stats kmeans #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> cluster <- function(x, cNum, clusterMethod="kmeans", clusterIterations=100, clusterStartNum=30){ #get cluster results using method on data x using cluster tuning paramater c_tune if(clusterMethod == "kmeans"){ kmx <- kmeans(t(x), cNum, iter.max=clusterIterations, nstart=clusterStartNum) return_val <- list(cluster=kmx$cluster, centers=t(kmx$centers), av_dist=(kmx$withinss/kmx$size)/dim(x)[1], size=kmx$size)
}
else if(clusterMethod == "kmeanspp"){
kmx <- kcca(t(x),cNum,control=list(initcent="kmeanspp",iter.max=clusterIterations))
centers <- t(attributes(kmx)$centers) #center_trans <- centers[,attributes(kmx)$cluster]
#av_dist <- tapply(colMeans((x-center_trans)^2),attributes(kmx)$cluster,mean) return_val <- attributes(kmx)$cluster
}
return_val
}

#' randomly assign n samples to k groups
#' @param n number of samples
#' @param k number of groups
#' @return Returns assignments of n into k groups
randomly_assign <- function(n,k){
#randomly assign n samples into k groups
small_set <- floor(n/k)
group_assign <- NULL
if(n %% k == 0){
group_assign <-  rep(seq(1,k),n/k)
} else{
remainder <- n %% k
for(i in 1:remainder){
group_assign <- c(group_assign, rep(i,small_set+1))
}
group_assign <- c(group_assign, rep(seq((i+1),k),small_set))
}
sample(group_assign)
}

#'matrix multiply
#'@param beta Matrix of coefficients.
#'@param x Design matrix.
#'@return Returns x times beta
matrix_multiply <- function(beta,x){
x%*%beta
}

#' predictions from a mcen model
#'@method predict mcen
#'@param object The mcen object.
#'@param newx A matrix of new observations.
#'@param ... Additional variables to be sent to predict.
#'@return Returns predictions for each beta of an mcen object
#'@examples
#' x <- matrix(rnorm(400),ncol=4)
#' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4)
#' y <- x%*%beta + rnorm(400)
#' mcen_fit <- mcen(x,y,ky=2,delta=1)
#' new_x <- matrix(rnorm(12),ncol=4)
#' mcen_preds <- predict(mcen_fit, new_x)
#'@export
predict.mcen <- function(object,newx,...){
newx <- cbind(1,newx)
preds <- lapply(object$beta,matrix_multiply,newx) preds } #' Calculates the out of sample likelihood for an mcen object #' @param obj The mcen object. #' @param test_x The matrix of test predictors. #' @param test_y The matrix of test responses. #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> pred_eval <- function(obj,test_x,test_y) { UseMethod("pred_eval") } #' Calculates sum of squared error between two vectors or matrices #' @param pred the predictions #' @param test_y the testing response values #' @return returns the sum of the squared differences between pred and test_y #' @author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> squared_error <- function(pred,test_y){ sum((pred-test_y)^2) } #' Calculates the prediction error for a mgauss_mcen object. #'@param obj The mgauss_mcen object. #'@param test_x The matrix of test predictors. #'@param test_y The matrix of test responses. #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> #'@importFrom stats predict pred_eval.mgauss_mcen <- function(obj,test_x,test_y){ predicted_vals <- predict(obj,test_x) pred_error <- sapply(predicted_vals,squared_error,test_y) pred_error } #' Calculates out of sample error on the binomial likelihood #'@param pred The predicted values. #'@param test_y The test response values. #'@importFrom faraway ilogit #'@author Brad Price <brad.price@@mail.wvu.edu> vl_binom<-function(pred,test_y){ pred=as.matrix(pred) out=which(ilogit(pred)=="NaN") pred[out]=1-10^-10 pred[which(ilogit(pred)>1-10^-10)]=1-10^-10 pred[which(ilogit(pred)<10^-10)]=10^-10 foo=sum(log(ilogit(pred)^test_y)+log((1-ilogit(pred))^(1-test_y))) return(-foo) } #' Evaluates prediction error for multiple binomial responses. #'@param obj The mbinom_mcen object. #'@param test_x A matrix of the test predictors. #'@param test_y A matrix of the test responses. #' @author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> #'@importFrom stats predict pred_eval.mbinom_mcen<-function(obj,test_x,test_y){ predicted_vals<-predict(obj,test_x) vl<-sapply(predicted_vals,vl_binom,test_y) return(vl) } #' Cross validation for mcen function #' @param x Matrix set of predictors. #' @param y Matrix set of responses. #' @param family The exponential family the response corresponds to. #' @param ky A vector with the number of possible clusters for y. #' @param gamma_y Set of tuning parameter for clustering penalty in response categories. #' @param nfolds Number of folds used in the cross-validation. #' @param folds A vector of length n, where this identifies what fold of the kfold cross validation each observation belongs to. #' @param cluster_y a priori definition of clusters. If clusters are provided they will remain fixed and are not estimated. Objective function is then convex. #' @param n.cores Number of cores used for parallel processing. #'@param ... The variables passed to mcen #' @return Returns a cv.mcen object. #' \item{models}{A list of mcen objects.} #' \item{cv}{Cross validation results.} #' \item{ky}{The same value as the input ky.} #' \item{gamma_y}{The same value as the input gamma_y.} #'@examples #'\donttest{ #' x <- matrix(rnorm(400),ncol=4) #' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4) #' y <- x%*%beta + rnorm(400) #' cv_fit <- cv.mcen(x,y,ky=2) #'} #'@references Price, B.S. and Sherwood, B. (2018). A Cluster Elastic Net for Multivariate Regression. arXiv preprint arXiv:1707.03530. \url{http://arxiv-export-lb.library.cornell.edu/abs/1707.03530}. #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> #'@export #'@importFrom parallel mclapply cv.mcen <-function(x,y,family="mgaussian",ky=seq(2,4),gamma_y=seq(.1,5.1,.5),nfolds=10,folds=NULL,cluster_y=NULL,delta=NULL, n.cores=1,...){ if(is.null(cluster_y)==FALSE){ ky = length(unique(cluster_y)) warning("ky does not match number of clusters in cluster_y, cluster_y then overrides ky") } clus=length(ky) gamma_yl=length(gamma_y) ResOut <- NULL#vector("list",clus) if(is.null(folds)){ folds <- randomly_assign(dim(x)[1],nfolds) } Eval<-function(del,gammay,kyn,fam,...){ #print(paste("working on eval for delta", del, "and gamma", gammay)) pse <- 0 for(s in unique(folds)){ #print(paste("working on fold", s)) Set=which(folds==s) #bet=mcen(x[-Set,],y[-Set,],delta=del,gamma_y=gammay,ky=kyn,...) if(fam=="mbinomial"){ train_model <- mcen(x[-Set,],y[-Set,],family="mbinomial",delta=del,gamma_y=gammay,ky=kyn,cluster_y=cluster_y,...) }else{ train_model <- mcen(x[-Set,],y[-Set,],delta=del,gamma_y=gammay,ky=kyn,cluster_y=cluster_y,...) } pse <- pse + pred_eval(train_model,x[Set,],y[Set,]) } pse } mc_Eval <- function(mcen_obj,family,...){ delta_vals <- mcen_obj$delta
gamma_y <- mcen_obj$gamma_y ky <- mcen_obj$ky
#print(paste("Evaluating model of ky working on ", ky, " gammay is ", gamma_y))
pe <- Eval(delta_vals,gamma_y,ky,fam=family,...)
cbind(ky,gamma_y,delta_vals,pe)
}

mc_mcen <- function(tune_vals,family,...){
#print(paste("Fitting full model ky",tune_vals[1],"gammay",tune_vals[2]))
mcen(x,y,ky=tune_vals[1],gamma_y=tune_vals[2],cluster_y=cluster_y,family=family,...)
}

mcen_list <- list()
mcen_num <- 1
if(n.cores == 1){
for(k in 1:clus){
k_gamma_results <- NULL
for(g in 1:gamma_yl){
#print(paste0("working on k ", k, " and g ", g))
mcen_full <- mcen(x,y,family,ky=ky[k], gamma_y=gamma_y[g],cluster_y=cluster_y,delta=delta,...)
mcen_list[[mcen_num]] <- mcen_full
train_delta <- mcen_full$delta eval_results <- Eval(train_delta,gamma_y[g],ky[k],fam=family,...) ResOut <- rbind(ResOut, cbind(ky[k],gamma_y[g],train_delta,eval_results,mcen_num)) mcen_num <- mcen_num + 1 } } } else{ ky_vals <- sort(rep(ky,gamma_yl)) gamma_y_vals <- rep(gamma_y,clus) tune_mat <- rbind(ky_vals,gamma_y_vals) tune_list <- split(tune_mat, rep(1:ncol(tune_mat), each = nrow(tune_mat))) mcen_list <- mclapply(tune_list,mc_mcen,mc.cores=n.cores,family=family,...) ResOut <- mclapply(mcen_list,mc_Eval,mc.cores=n.cores,family=family,...) ResOut <- do.call(rbind,ResOut) model_num <- rep(seq(1,length(mcen_list)),each=length(mcen_list[[1]]$delta))
ResOut <- cbind(ResOut,model_num)

}

#returns list of all models and cross-validation erorr
#later functions will handle coefficients and predictions from this object
Out = list(models=mcen_list,cv=ResOut,ky=ky,gamma_y=gamma_y)
class(Out) <- c(paste0(family,"_cv.mcen"),"cv.mcen")
return(Out)

}

#' Prints nice output for an mcen object.
#'@param x The mcen object.
#'@return Prints out some basic information about the mcen object.
#'@export
print.mcen <- function(x,...){
cat("\n mcen object \n")
cat("class:")
print(class(x))
cat("\n gamma_y value")
print(x$gamma_y) cat("\n Number of groups for response:") print(x$ky)
}

#' Prints nice output for a cv.mcen object.
#'@param x The cv.mcen object.
#'@return Prints out information about where the cv.mcen object was minimized.
#'@export
print.cv.mcen <- function(x,...){
cat("\nCross Validation Error minimized at:\n")
min_spot <- which.min(x$cv[,4]) print(x$cv[min_spot,])
}

#' Makes predictions from the model with the smallest cross-validation error.
#'@method predict cv.mcen
#'@param object The cv.mcen object.
#'@param newx The X matrix of predictors.
#'@param ... Additional parameters to be sent to predict.
#'@return Returns the predicted values from the model with the smallest cross-validation error.
#'@examples
#'\donttest{
#' x <- matrix(rnorm(400),ncol=4)
#' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4)
#' y <- x%*%beta + rnorm(400)
#' mcen_fit <- cv.mcen(x,y,ky=2,gamma_y=3)
#' new_x <- matrix(rnorm(12),ncol=4)
#' mcen_preds <- predict(mcen_fit, new_x)
#'}
#'@export
#'@importFrom stats coef
predict.cv.mcen <- function(object,newx,...){
newx <- cbind(1,newx)
preds <- newx %*% coef(object)
if("mbinomial_cv.mcen" %in% class(object)){
ilogit(preds)
} else{
preds
}
}

#' Returns the cluster values from a cv.mcen object.
#'@param obj The cv.mcen object.
#'@return Returns the clusters from the model with the smallest cross-validation error.
#'@examples
#'\donttest{
#' x <- matrix(rnorm(400),ncol=4)
#' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4)
#' y <- x%*%beta + rnorm(400)
#' mcen_fit <- cv.mcen(x,y,ky=2,gamma_y=3)
#' mcen_cluster <- cluster.vals(mcen_fit)
#' }
#'@export
cluster.vals <- function(obj){
min_spot <- which.min(obj$cv[,4]) min_delta <- obj$cv[min_spot,3]
min_model <- obj$cv[min_spot,5] best_model <- obj$models[[min_model]]
delta_spot <- which(best_model$delta==min_delta) best_model$y_clusters[[delta_spot]]
}

#' Returns the coefficients from the cv.mcen object with the smallest cross-validation error.
#'@method coef cv.mcen
#'@param object The cv.mcen object.
#'@param ... Additional values to be passed.
#'@return The matrix of coefficients for the best MCEN model as determined by cross-validation.
#'@examples
#'\donttest{
#' x <- matrix(rnorm(400),ncol=4)
#' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4)
#' y <- x%*%beta + rnorm(400)
#' mcen_fit <- cv.mcen(x,y,ky=2,gamma_y=3)
#' best_coef <- coefficients(mcen_fit)
#' }
#'@export
#'@importFrom stats coef
coef.cv.mcen <- function(object,...){
min_spot <- which.min(object$cv[,4]) min_delta <- object$cv[min_spot,3]
min_model <- object$cv[min_spot,5] best_model <- object$models[[min_model]]
best_beta <- coef(best_model,min_delta)
best_beta
}

#' Returns the coefficients from an mcen object.
#'@method coef mcen
#'@param object The mcen object.
#'@param delta The L1 tuning parameter
#'@param ... Additional values to pass on.
#'@return The matrix of coefficients.
#'@examples
#' x <- matrix(rnorm(400),ncol=4)
#' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4)
#' y <- x%*%beta + rnorm(400)
#' mcen_fit <- mcen(x,y,ky=2,gamma_y=3,delta=c(1,2))
#' best_coef <- coefficients(mcen_fit,delta=1)
#'@export
coef.mcen <- function(object,delta=NULL,...){
if(is.null(delta)){
object$beta } else{ delta_pos <- which(object$delta == delta)
object$beta[[delta_pos]] } } #' Gets the index position for the model with the smallest cross-validation error. #'@param model The cv.mcen object. #'@return Returns the index for the model with the smallest cross-validation error. #'@author Ben Sherwood <ben.sherwood@@ku.edu>, Brad Price <brad.price@@mail.wvu.edu> #'@examples #'\donttest{ #' x <- matrix(rnorm(400),ncol=4) #' beta <- beta <- matrix(c(1,1,0,0,0,0,-1,-1,0,0,-1,-1,1,1,0,0),ncol=4) #' y <- x%*%beta + rnorm(400) #' mcen_fit <- cv.mcen(x,y,ky=2,gamma_y=3) #' get_best_cvm(mcen_fit) #' } #'@export get_best_cvm <- function(model){ min(model$cvm)
}

#'Creates the probabilities and working response for the glmnet update for a given response with a binomial family
#' @param X the matrix of predictors.
#' @param Beta current iteration of the regression coefficients
#' @param Y is the matrix of responses
#' @param r the response of interest
#' result is a list of things needed for the working response in glmnet
CalcHorseEBin<-function(X,Beta,Y,r){
Probs=ilogit(X%*%Beta[,r])
set=which(Probs<=10^-5 | Probs>=(1-10^-5))
W=Probs*(1-Probs)
W[set]=10^-5
Z=X%*%Beta[,r]+(Y[,r]-Probs)/W
V=sqrt(W)*Z
Tm=apply(X,2,function(x){sqrt(W)*x})
R=t(X)%*%X
return(list(P=Probs,W=W,Z=Z,V=V,Tm=Tm,R=R))
}

#' Creates the the working response for all responses for glmnet binomial family
#' @param X the matrix of predictors.
#' @param Beta current iteration of the regression coefficients
#' @param Y is the matrix of responses
#' result is the list of vectors needed for the working responses in glmnet
CalcHorseBin<-function(Y,X,Beta){
TMm=NULL
TmV=NULL
Cm=NULL
Vm=NULL
R=NULL
for(l in 1:dim(Y)[2]){
out<-CalcHorseEBin(X,Beta,Y,l)
TMm=cbind(TMm,out$Tm) Cm=cbind(Cm,t(out$Tm)%*%out$Tm) Vm=cbind(Vm,out$V)
TmV=cbind(TmV,t(out$Tm)%*%out$V)
##R=cbind(R,out$R) } R=out$R
return(list(Tm=TMm,C=Cm,V=Vm,R=R,TmV=TmV))
}

#' The workhorse function for the binomial updates in mcen.  It uses IRWLS glmnet updates to solve the regression problem.
#' @param Y the matrix of responses
#' @param X the matrix of predictors with the intercept included
#' @param delta the tuning parameter for the lasso penalty
#' @param gamma_y the tuning parameter for the ridge fusion penalty
#' @param y_clusters the cluster assignments from the provided clustering algorithm
#' @param set_length the size of each cluster corresponding to a given response.  r dimensions with each element containing the cluster size of that responses cluster.
#' @param eps the tolerance for conversion normally 1e-5
#' @param maxiter the maximum number of iterations
#'@return Returns a matrix of coefficients
#'@useDynLib mcen BinUp
bin_horse<-function(Y,X,delta,gamma_y,y_clusters,set_length,eps,maxiter){
Beta01=matrix(0,nrow=dim(X)[2],ncol=dim(Y)[2])
BetaI2=Beta01+1
iter=0
while(sum((Beta01-BetaI2)^2)/length(Beta01)>eps && iter<maxiter){
BetaI2=Beta01
Mine<-CalcHorseBin(Y,X,BetaI2)
Out<-.C("BinUp",TmV=as.double((as.vector(Mine$TmV))),Tm=as.double(as.vector(Mine$Tm)),V=as.double(as.vector(Mine$V)),C=as.double(as.vector(Mine$C)),R=as.double(as.vector(Mine$R)), Beta=as.double(as.vector(BetaI2)),r=as.integer(dim(Y)[2]),p=as.integer(dim(X)[2]),n=as.integer(dim(X)[1]), y_clusters=as.double(as.vector(y_clusters)), sl=as.double(as.vector(set_length)), gamma_y=as.double(gamma_y), delta=as.double(delta), iter=as.double(maxiter),eps=as.double(eps),third=as.double(as.vector(Beta01)),fourth=as.double(as.vector(Beta01)), sixth=as.double(as.vector(Beta01)),first=as.double(as.vector(Beta01)),second=as.double(as.vector(Beta01)),fifth=as.double(as.vector(Beta01))) Beta01=matrix(Out$Beta,byrow = FALSE,ncol=dim(Y)[2],nrow=dim(X)[2])
##print(sum((Beta01-BetaI2)^2)/length(Beta01))
iter=iter+1
}
return(list(Beta=Beta01))
}

#' Calculates cluster assignment and coefficient estimates for a binomial mcen.
#'@param beta Initial estimate of coefficients.
#'@param delta Tuning parameter for L1 penalty.
#'@param y Matrix of responses.
#'@param x Matrix of predictors.
#'@param family type of likelihood used two options "mgaussian" or "mbinomial"
#'@param ky Number of clusters used for grouping response variables.
#'@param gamma_y Tuning parameter for the penalty between fitted values for responses in the same group.
#'@param eps Convergence criteria
#'@param clusterMethod Which clustering method was used, currently support kmeans or kmeanspp
#'@param clusterIterations Number of iterations for cluster convergence
#'@param clusterStartNum Number of random starting points for clustering
#'@param cluster_y An a priori definition of clusters. If clusters are provided they will remain fixed and are not estimated. Objective function is then convex.
#'@param max_iter The maximum number of iterations for estimating the coefficients
mcen_bin_workhorse<-function(beta,delta=NULL,y,x,family="mbinomial",ky=NULL,gamma_y=1,eps=.00001,
clusterMethod="kmeans",clusterIterations=100,clusterStartNum=30,cluster_y=NULL,max_iter=10){
beta_dim <- dim(beta)
p <- beta_dim[1]
r <- beta_dim[2]
if(sum(beta==0)==p*r){
#if initial estimate is totally sparse then mcen update is not needed. Philosophically, I'm ok with this because mcen idea is
# to group coefficient estimates, if all coefficients are zero there is nothing to group and nothing for mcen to do.
return_val <- list(beta = as(beta,"sparseMatrix"), clusters=rep(0,r))
return( return_val)
} else{
fitted_vals <- x%*%beta
fitted_vals<-as.matrix(fitted_vals)
iter <- 0
fitting_done <- FALSE
while(!fitting_done){
#print(paste("fitting done is ", fitting_done, " and iter is ", iter))
num_potential_centers <- ncol(unique(fitted_vals,MARGIN = 2))
if(num_potential_centers < ky){
fitting_done <- TRUE
warning(paste0("Algorithm stopped early for gamma_y ", gamma_y, ", delta ", delta, " and ky ", ky, " because number of unqiue fitted values is less than ", ky, " Common for large values of deltas that cause spares solution and little variability in fitted values.\n"))
clusters_used <- rep(0,r)
} else if(is.null(cluster_y)==FALSE){
clusters_used <- cluster_y
} else{
y_clusters <- cluster(fitted_vals,ky,clusterMethod,clusterIterations,clusterStartNum)
clusters_used <- y_clusters$cluster } if(iter > 0){ if(is.null(cluster_y) && fitting_done==FALSE){ #fitting_done=SetEq(y_clusters2$centers,y_clusters$centers) fitting_done=SetEq(y_clusters2$cluster,y_clusters$cluster) } else{ fitting_done <- TRUE } } if(fitting_done == FALSE){ sl=NULL for(m in 1:length(clusters_used)){ sl[m]=length(which(clusters_used==clusters_used[m])) } beta<-bin_horse(Y=y,X=x,delta = delta,gamma_y = gamma_y,y_clusters =clusters_used ,set_length = sl,eps=eps,maxiter = max_iter)$Beta
}
if(iter == max_iter){
fitting_done <- TRUE
}
iter <- iter + 1
if(is.null(cluster_y)==TRUE && fitting_done==FALSE){
y_clusters2 <- y_clusters
}
fitted_vals=x%*%beta
}
if(is.null(cluster_y)==FALSE){
fitting_done <- TRUE #only one iteration of fitting done if clusters are prespecified
}
}
return_val <- list(beta = as(beta,"sparseMatrix"), clusters=clusters_used)
return( return_val)
}


## Try the mcen package in your browser

Any scripts or data that you put into this service are public.

mcen documentation built on July 23, 2020, 5:08 p.m.