R/summarizeFarmsLaplaceExact.R

In cn.farms: cn.FARMS - factor analysis for copy number estimation

#' Summarization Laplacian approach with exact computation
#'
#' This function implements an exact Laplace FARMS algorithm. 
#' 
#' @param probes A matrix with numeric values.
#' @param mu Hyperparameter value which allows to quantify different aspects of potential prior knowledge. Values near zero assumes that most positions do not contain a signal, and introduces a bias for loading matrix elements near zero. Default value is 0 and it's recommended not to change it.
#' @param weight Hyperparameter value which determines the influence of the Gaussian prior of the loadings
#' @param weightSignal Hyperparameter value on the signal.
#' @param weightZ Hyperparameter value which determines how strong the Laplace prior of the factor should be at 0. Users should be aware, that a change of weightZ in comparison to the default parameter might also entail a need to change other parameters. Unexperienced users should not change weightZ.
#' @param weightProbes Parameter (TRUE/FALSE), that determines, if the number of probes should additionally be considered in weight. If TRUE, weight will be modified.
#' @param cyc Number of cycles. If the length is two, it is assumed, that a minimum and a maximum number of cycles is given. If the length is one, the value is interpreted as the exact number of cycles to be executed (minimum == maximum).
#' @param tol States the termination tolerance if cyc[1]!=cyc[2]. Default is 0.00001.
#' @param weightType Flag, that is used to summarize the probes of a sample. 
#' @param centering States how the data should be centered ("mean", "median"). Default is median.
#' @param rescale Parameter (TRUE/FALSE), that determines, if moments in exact Laplace FARMS are rescaled in each iteration. Default is FALSE.
#' @param backscaleComputation Parameter (TRUE/FALSE), that determines if the moments of hidden variables should be reestimated after rescaling the parameters.
#' @param maxIntensity Parameter (TRUE/FALSE), that determines if the expectation value (=FALSE) or the maximum value (=TRUE) of p(z|x_i) should be used for an estimation of the hidden varaible. 
#' @param refIdx index or indices which are used for computation of the centering
#' @param ... Further parameters for expert users.
#' 
#' @return A list including:
#' the found parameters: lambda0, lambda1, Psi
#' 
#' the estimated factors: z (expectation), maxZ (maximum)
#' 
#' p: log-likelihood of the data given the found lambda0, lambda1, 
#' Psi (not the posterior likelihood that is optimized)
#' 
#' varzx: variances of the hidden variables given the data
#' 
#' KL: Kullback Leibler divergences between between posterior and prior 
#' distribution of the hidden variables
#' 
#' IC: Information Content considering the hidden variables and data
#' 
#' ICtransform: transformed Information Content
#' 
#' Case: Case for computation of a sample point (non-exception, special exception)
#' 
#' L1median: Median of the lambda vector components
#' 
#' intensity: back-computed summarized probeset values with mean correction
#' 
#' L_z: back-computed summarized probeset values without mean correction
#' 
#' rawCN: transformed values of L_z
#' 
#' SNR: some additional signal to noise ratio value
#' 
#' 
#' @author Andreas Mayr \email{mayr@@bioinf.jku.at} and Djork-Arne Clevert \email{okko@@clevert.de} and Andreas Mitterecker \email{mitterecker@@bioinf.jku.at}
#' @examples 
#' x <- matrix(rnorm(100, 11), 20, 5)
#' summarizeFarmsExact(x)
#' @export

summarizeFarmsExact <- function(
		probes, 
		mu = 1.0, 
		weight = 0.001, 
		weightSignal = 1.0,
		weightZ = 1.0, 
		weightProbes = TRUE,
		cyc = c(10, 10), 
		tol = 0.00001, 
		weightType = "mean", 
		centering = "median",
		rescale = FALSE,
		backscaleComputation = FALSE,
		maxIntensity = TRUE,
		refIdx,
		...) {
    
	eps1 <- 0.01
	eps2 <- 10^-10
	probes <- as.matrix(probes)
	sigmaZ <- 1 / sqrt(2 * weightZ)
	if(length(cyc) == 1)
		cyc <- c(cyc, cyc)
	additional <- list(...)
	
	if("initPsi" %in% names(additional)) {
        initPsi <- additional$initPsi
    } else {
        initPsi <- 0.1
    } 
	
	if("init" %in% names(additional)) {
        init <- additional$init
    } else {
        init <- "mean"
    }
	
	if("algorithm" %in% names(additional)) {
        algorithm <- additional$algorithm
    } else {
        algorithm <- "exact"
    }
	
	if(algorithm == "v" || algorithm == "ard") {
		if("boundedLapla" %in% names(additional)) {
            boundedLapla <- additional$boundedLapla
        } else {
            boundedLapla <- TRUE
        }
		
		if("spuriousCorrelation" %in% names(additional)) {
            spuriousCorrelation <- additional$spuriousCorrelation
        } else {
            spuriousCorrelation <- 1.0
        }
	}
	
	if(algorithm == "ard") {
		if("ard2Max" %in% names(additional))
			ard2Max <- additional$ard2Max
		else
			ard2Max <- -1
	}
	
	
	## Initialize variables
	Data <- probes
	DataT <- t(probes)
	n <- ncol(Data)
	dimension <- nrow(Data)
	
	if(algorithm == "mean") {
		allMedian <- median(Data) 
		twoMedian <- median(apply(Data, 1, median))
		allMean <- mean(Data) 
		twoMean <- mean(apply(Data, 1, mean))
		if (centering == "median") { mean.Data <- apply(Data, 1, median) }
		if (centering == "mean") { mean.Data <- rowMeans(Data) } 
		NData <- Data - mean.Data
		sd.Data <- sqrt(rowSums(NData^2) / n)
		NData <- NData / (sd.Data + 10^-3)
		return(list(colMean = colMeans(NData), colMedian = apply(NData, 2, median), 
                        allMedian = allMedian, twoMedian = twoMedian, 
                        allMean = allMean, twoMean = twoMean))
	}
	
	if (missing(refIdx)) {
		refIdx <- 1:n 
	}
	if (centering == "median") { 
		if (length(refIdx) == 1) {
			mean.Data <- Data[, refIdx]
		}
		else {
			mean.Data <- Biobase::rowMedians(Data[, refIdx])
		}
	}
	if (centering == "mean") { 
		if (length(refIdx) == 1) {
			mean.Data <- Data[, refIdx]
		}
		else {
			mean.Data <- rowMeans(Data[, refIdx])
		}
	}
	NData <- Data - mean.Data
	
	sd.Data <- sqrt(rowSums(NData^2) / n)
	NData <- NData / (sd.Data + 10^-3)
	NDataT <- t(NData)
	DataCov <- NData %*% NDataT / n
	DataCov[DataCov <= 0] <- 10^-3
	DataCov <- (DataCov + t(DataCov)) / 2
	
	myLambda <- rep(mu, dimension)
	PsiLambda <- rep(1.0 / weight, dimension)
	if(weightProbes) {
        PsiLambda <- rep(1.0 / (weight*dimension), dimension)
    }
	sigmaS <- 1/sqrt(weightSignal)
	
	s <- 1.0
	Psi <- initPsi * diag(DataCov)
	lambda <- sqrt(diag(DataCov) - Psi)
	lapla <- rep(1, n)
	if(init == "maxAbsSample") {
		sampleIdx <- which.max(abs(apply(NData, 2, mean)))
		lambda <- (abs(NData[, sampleIdx])/max(abs(NData[, sampleIdx])))*sqrt(1-initPsi)
		Psi <- diag(DataCov)-lambda^2
	} else if(init == "maxAbs") {
		maxAcrossSamples <- apply(NData, 1, function(x) max(abs(x)))
		lambda <- (maxAcrossSamples/max(maxAcrossSamples))*sqrt(1-initPsi)
		Psi <- diag(DataCov)-lambda^2
	} else if(init == "absEqual") {
		sumSamples <- abs(rowSums(NData))
		lambda <- (sumSamples/max(sumSamples))*sqrt(1-initPsi)
		Psi <- diag(DataCov)-lambda^2
	}
	
	NData2 <- NData^2
	
	#EM algorithm
	nrCyc <- cyc[2]
	PsiOld <- Psi
	
	results <- list()
	
	log_p_ges <- 0
	lambda_old <- lambda
	Psi_old <- Psi
	
	for (i in 1:cyc[2]) {
		## E-Step
		
		if (algorithm == "exact") {
			InvPsi <- 1 / Psi
			av <- rep(-0.5 * (t(lambda) %*% (InvPsi * lambda))[1], n)
			bv <- as.vector((lambda * InvPsi) %*% NData)
			cv <- -0.5 * colSums(NData2 * InvPsi)
			nv <- rep(1 / ((2 * pi)^(dimension / 2) * prod(Psi)^(1/2) * 2 * sigmaZ), n)
			moments <- .Call("momentsGauss", i, eps1, eps2, av, bv, cv, sigmaZ, 
                    nv, 1, 0, PACKAGE = "cn.farms")
			if(rescale) {
				sdmom <- sqrt(1 / n * sum(moments$moment2)) / (sqrt(2.0) * sigmaZ) + 10^-3
				moments$moment1 <- moments$moment1 / sdmom 
				moments$moment2 <- moments$moment2 / sdmom^2
			}
			avg_xEz <- as.vector(NData %*% moments$moment1 / n)
			avg_Ez2 <- mean(moments$moment2)
		} else if (algorithm == "v") {
			InvPsi <- 1 / Psi
			sigma2 <- 1 / (lapla + (t(lambda) %*% (InvPsi * lambda))[1])
			Ez <- as.vector((lambda * InvPsi) %*% NData) * sigma2
			Ez2 <- Ez^2 + sigma2
			avg_xEz <- as.vector(NData %*% Ez / n)
			avg_Ez2 <- mean(Ez2)
			
			## M-Step only for v
			lapla <- 1 / sqrt(Ez2)
			if(boundedLapla)
				lapla[lapla < spuriousCorrelation] <- spuriousCorrelation 
		} else if (algorithm == "ard") {
			InvPsi <- 1 / Psi
			sigma2 <- 1 / (lapla + (t(lambda) %*% (InvPsi * lambda))[1])
			Ez <- as.vector((lambda * InvPsi) %*% NData) * sigma2
			Ez2 <- Ez^2 + sigma2
			avg_xEz <- as.vector(NData %*% Ez / n)
			avg_Ez2 <- mean(Ez2)
			
			## M-Step only for ard
			lapla <- 1 / Ez2
			if(boundedLapla)
				lapla[lapla < spuriousCorrelation] <- spuriousCorrelation 
		} else if (algorithm == "g") {
			PsiM1_Lambda <- lambda / Psi
			sigmaZ2 <- 1 / (1 + sum(lambda * PsiM1_Lambda))
			avg_xEz <- (as.vector(DataCov %*% PsiM1_Lambda) * sigmaZ2)
			avg_Ez2 <- sum(PsiM1_Lambda * avg_xEz) * sigmaZ2 + sigmaZ2
		}
		
		lambda <- lambda_old
		
		## M-Step
		Psi_PsiLambdaM1 <- Psi / PsiLambda
		
		s <- sum(lambda*avg_xEz/Psi)/(1/sigmaS^2+sum(lambda^2/Psi)*avg_Ez2)
		lambda <- (avg_xEz * s + Psi_PsiLambdaM1 * myLambda) / (avg_Ez2 * s^2 * rep(1, dimension) + Psi_PsiLambdaM1)
		#lambda <- (avg_xEz + Psi_PsiLambdaM1 * myLambda) / (avg_Ez2 * rep(1, dimension) + Psi_PsiLambdaM1)
		lambda <- ifelse(lambda < 0.0, rep(0, length(lambda)), lambda)
		Psi <- diag(DataCov) - avg_xEz * s * lambda + Psi_PsiLambdaM1 * (myLambda - lambda) * lambda
		#Psi <- diag(DataCov) - avg_xEz * lambda + Psi_PsiLambdaM1 * (myLambda - lambda) * lambda
		Psi[Psi < 10^-3] <- 10^-3
		
		lambda_old <- lambda
		lambda <- lambda*s
		
		if(i > cyc[1] && max(abs(Psi - PsiOld)) / max(abs(PsiOld)) < tol) {
			nrCyc <- i + 1
			break
		}
		
		PsiOld <- Psi
	}
    
	InvPsi <- 1 / Psi
	av <- rep(-0.5 * (t(lambda) %*% (InvPsi * lambda))[1], n)
	bv <- as.vector((lambda * InvPsi) %*% NData)
	cv <- -0.5 * colSums(NData2 * InvPsi)
	nv <- rep(1 / ((2 * pi)^(dimension/2) * prod(Psi)^(1 / 2) * 2 * sigmaZ), n)
	moments <- .Call("momentsGauss", i, eps1, eps2, av, bv, cv, sigmaZ, nv, 1, 0, PACKAGE="cn.farms")
	
	if(algorithm == "exact") {
		z <- moments$moment1
		maxZ <- moments$max
		varzx <- moments$moment2 - moments$moment1^2
		KL <- moments$CrossEntropy - moments$Entropy
		IC <- log(2 * exp(1.0) * sigmaZ) - moments$Entropy  
		sdz <- sqrt(1 / n * sum(moments$moment2)) / (sqrt(2.0) * sigmaZ)
	} else if ((algorithm == "v") || (algorithm == "ard")) {
		sigma2 <- 1 / (lapla + (t(lambda) %*% (InvPsi * lambda))[1])
		z <- as.vector((lambda * InvPsi) %*% NData) * sigma2
		maxZ <- z
		if(algorithm == "ard" && (ard2Max > 0)) {
			disc <- 8 * av + bv^2 * ard2Max
			maxZ[] <- 0.0
			maxZP <- maxZ
			maxZM <- maxZ
			d0 <- disc > 0.0
			maxZP[d0] <- 2.0 / (bv[d0] * ard2Max + sqrt(ard2Max * disc[d0]))
			maxZM[d0] <- 2.0 / (bv[d0] * ard2Max - sqrt(ard2Max * disc[d0]))
			d0P <- (disc > 0.0) & (maxZP > 0.0)
			maxZ[d0P] <- maxZM[d0P]
			d0M <- (disc > 0.0) & (maxZM < 0.0)
			maxZ[d0M] <- maxZP[d0M]
		}
		varzx <- sigma2
		KL <- (0.5 * log(2 * pi) - 0.5 * log(lapla) + ((z^2 + sigma2) * lapla) / 2.0) - 0.5 * (log(2 * sigma2 * pi * exp(1)))
		IC <- 0.5 * log(1.0 + (t(lambda) %*% (InvPsi * lambda))[1] / lapla) 
   		sdz <- sqrt(1 / n * sum(z^2 + varzx))
	} else if (algorithm == "g") {
		PsiM1_Lambda <- lambda / Psi
		sigma2 <- 1 / (1 + sum(lambda * PsiM1_Lambda))
		z <- as.vector(NDataT %*% PsiM1_Lambda) * sigma2
		maxZ <- z
		varzx <- rep(sigma2, ncol(Data))    
		KL <- (0.5 * log(2 * pi) + ((z^2 + sigma2)) / 2.0) - 0.5*(log(2 * sigma2 * pi * exp(1)))
		IC <- 0.5 * log(1.0 + (t(lambda) %*% (InvPsi * lambda))[1])   
   		sdz <- sqrt(1 / n * sum(z^2 + varzx))
	}
	ICtransform <- 1 / exp(IC * 2.0)
	
	if(sdz == 0.0)
		sdz <- 1
	
	z <- z / sdz
	maxZ <- maxZ / sdz
	lambda <- lambda * sdz
	if(algorithm == "v") {
        lapla <- lapla / sdz
    } else if (algorithm == "ard") {
        lapla <- lapla / sdz^2
    }
	
	avRET <- rep(-0.5 * (t(lambda) %*% (InvPsi * lambda))[1], n)
	bvRET <- as.vector((lambda * InvPsi) %*% NData)
	cvRET <- -0.5 * colSums(NData2 * InvPsi)
	nvRET <- rep(1 / ((2 * pi)^(dimension/2) * prod(Psi)^(1 / 2) * 2 * sigmaZ), n)
	momentsRET <- .Call("momentsGauss", i, eps1, eps2, avRET, bvRET, cvRET, sigmaZ, nvRET, 1, 0, PACKAGE="cn.farms")
	cvCOMP <- -0.5 * colSums(NData2 * rep(1.0, dimension))
	nvCOMP <- rep(1 / ((2 * pi)^(dimension/2)), n)
	MLQ <- pchisq(2 * (sum(log(momentsRET$normConst) - (cvCOMP + log(nvCOMP)))), dimension)
	log_p_ges <- sum(log(momentsRET$normConst))
	
	if(algorithm == "exact") {
		zScaled <- momentsRET$moment1
		maxZScaled <- momentsRET$max
		varzxScaled <- momentsRET$moment2 - momentsRET$moment1^2
		KLScaled <- momentsRET$CrossEntropy - momentsRET$Entropy
		ICScaled <- log(2 * exp(1.0) * sigmaZ) - momentsRET$Entropy    
	} else if ((algorithm == "v")||(algorithm == "ard")) {
		sigma2Scaled <- 1 / (lapla + (t(lambda) %*% (InvPsi * lambda))[1])
		zScaled <- as.vector((lambda * InvPsi) %*% NData) * sigma2Scaled
		maxZScaled <- zScaled
		if(algorithm == "ard" && (ard2Max > 0)) {
			disc <- 8 * avRET + bvRET^2 * ard2Max
			maxZScaled[] <- 0.0
			maxZScaledP <- maxZScaled
			maxZScaledM <- maxZScaled
			d0 <- disc > 0.0
			maxZScaledP[d0] <- 2.0 / (bvRET[d0] * ard2Max + sqrt(ard2Max * disc[d0]))
			maxZScaledM[d0] <- 2.0 / (bvRET[d0] * ard2Max - sqrt(ard2Max * disc[d0]))
			d0P <- (disc > 0.0) & (maxZScaledP > 0.0)
			maxZScaled[d0P] <- maxZScaledM[d0P]
			d0M <- (disc > 0.0) & (maxZScaledM < 0.0)
			maxZScaled[d0M] <- maxZScaledP[d0M]
		}
		varzxScaled <- sigma2Scaled
		KLScaled <- (0.5 * log(2 * pi) - 0.5 * log(lapla) + ((zScaled^2 + sigma2Scaled) * lapla) / 2.0) - 0.5*(log(2 * sigma2Scaled * pi * exp(1)))
		ICScaled <- 0.5 * log(1.0 + (t(lambda) %*% (InvPsi * lambda))[1] / lapla)    
	} else if (algorithm == "g") {
		PsiM1_Lambda <- lambda / Psi
		sigma2Scaled <- 1 / (1 + sum(lambda * PsiM1_Lambda))
		zScaled <- as.vector(NDataT %*% PsiM1_Lambda) * sigma2Scaled
		maxZScaled <- zScaled
		varzxScaled <- rep(sigma2Scaled, ncol(Data))    
		KLScaled <- (0.5 * log(2 * pi) + ((zScaled^2 + sigma2Scaled)) / 2.0) - 0.5*(log(2 * sigma2Scaled * pi * exp(1)))
		ICScaled <- 0.5 * log(1.0 + (t(lambda) %*% (InvPsi * lambda))[1])        
	}
	ICtransformScaled <- 1 / exp(ICScaled * 2.0)
    
	if(backscaleComputation == TRUE) {
		z <- zScaled
		maxZ <- maxZScaled
	}
	
	if(maxIntensity) {
        zint <- maxZ
    } else {
        zint <- z
    }
	
	lambdaNormalized <- lambda
	PsiNormalized <- Psi
	lambda0 <- mean.Data
	lambda1 <- lambda * sd.Data
	Psi <- Psi * sd.Data^2
	rm("lambda")
	
	names(lambdaNormalized) <- rownames(probes)
	names(PsiNormalized) <- rownames(probes)
	names(z) <- colnames(probes)
	names(maxZ) <- colnames(probes)
	names(lambda0) <- rownames(probes)
	names(lambda1) <- rownames(probes)
	names(Psi) <- rownames(probes)
	
	lambdaMeanNormalized <- mean(lambdaNormalized)
	lambdaMedianNormalized <- median(lambdaNormalized)
	lambda1Mean <- mean(lambda1)
	lambda1Median <- median(lambda1)
	SNR <- 1/(1+sum(lambda1^2/Psi))
	
	#linear
	PsiLL <- (lambda1^2 / Psi)
	sumPsiLL <- sum(PsiLL)
	if (sumPsiLL == 0) { 
		sumPsiLL <- 1
	}
	propPsiLL <- PsiLL / sumPsiLL
	L_c_Linear <- as.vector(crossprod(lambda1, propPsiLL)) * zint
	L_c_Linear_Normalized <- as.vector(crossprod(lambdaNormalized, propPsiLL)) * zint
	expressLinear <- mean(lambda0) + L_c_Linear
	
	#square
	PsiLL <- (lambda1^2 / Psi)^2
	sumPsiLL <- sum(PsiLL)
	if (sumPsiLL == 0) { 
		sumPsiLL <- 1
	}
	propPsiLL <- PsiLL / sumPsiLL
	L_c_Square <- as.vector(crossprod(lambda1, propPsiLL)) * zint
	L_c_Square_Normalized <- as.vector(crossprod(lambdaNormalized, propPsiLL)) * zint
	expressSquare <- mean(lambda0) + L_c_Square
	
	#softmax
	PsiLL <- exp(lambda1)
	sumPsiLL <- sum(PsiLL)
	if (sumPsiLL == 0) { 
		sumPsiLL <- 1
	}
	propPsiLL <- PsiLL / sumPsiLL
	L_c_Softmax <- as.vector(crossprod(lambda1, propPsiLL)) * zint
	L_c_Softmax_Normalized <- as.vector(crossprod(lambdaNormalized, propPsiLL)) * zint
	expressSoftmax <- mean(lambda0) + L_c_Softmax
	
	#mean
	L_c_Mean <- mean(lambda1)*zint
	L_c_Mean_Normalized <- mean(lambdaNormalized)*zint
	expressMeanBC <- mean(lambda0)+mean(lambda1)*zint
	expressBCMean <- colMeans(lambda0+lambda1%o%zint)
	
	#median
	L_c_Median <- median(lambda1)*zint
	L_c_Median_Normalized <- median(lambdaNormalized)*zint
	expressMedianBC <- median(lambda0)+median(lambda1)*zint
	expressBCMedian <- apply(lambda0+lambda1%o%zint, 2, median)
	
	if (weightType == "linear") {
		L_c <- L_c_Linear
		L_c_Normalized <- L_c_Linear_Normalized
		express <- expressLinear
	} else if(weightType == "square") {
		L_c <- L_c_Square
		L_c_Normalized <- L_c_Square_Normalized
		express <- expressSquare
	} else if (weightType == "softmax") {
		L_c <- L_c_Softmax
		L_c_Normalized <- L_c_Softmax_Normalized
		express <- expressSoftmax
	} else if ((weightType == "medianBC")||(weightType == "median")) {
		L_c <- L_c_Median
		L_c_Normalized <- L_c_Median_Normalized
		express <- expressMedianBC
	} else if ((weightType == "meanBC")||(weightType == "mean")) {
		L_c <- L_c_Mean
		L_c_Normalized <- L_c_Mean_Normalized
		express <- expressMeanBC
	} else if (weightType == "BCMedian") { 
		L_c <- L_c_Median
		L_c_Normalized <- L_c_Median_Normalized
		express <- expressBCMedian
	} else if (weightType == "BCMean") {
		L_c <- L_c_Mean
		L_c_Normalized <- L_c_Mean_Normalized
		express <- expressBCMean
	}
	
	cnINIMean <- mean(lambda1)*zint
	cnINIMedian <- median(lambda1)*zint
	cnINIMeanSignal <- cnINIMean*((sum(lambda1^2/Psi))/(1+sum(lambda1^2/Psi)))
	cnINIMedianSignal <- cnINIMedian*((sum(lambda1^2/Psi))/(1+sum(lambda1^2/Psi)))
	
	return(list(
                    sdData=sd.Data,
			        lambdaNormalized=lambdaNormalized,
			        PsiNormalized=PsiNormalized,
			        lambda0=lambda0, 
			        lambda1=lambda1, 
			        Psi=Psi,
			        z=z, 
			        maxZ=maxZ, 
			        ExactMaxZ=moments$max,
			        ExactMaxZScaled=momentsRET$max,
			        ExactP=log_p_ges, 
			        ExactMLQ=MLQ,
			        ExactCase=momentsRET$Case,
			        IndividualExactKL=momentsRET$CrossEntropy - momentsRET$Entropy,
			        SummaryExactKL=mean(momentsRET$CrossEntropy - momentsRET$Entropy),
			        IndividualExactKLScaled=moments$CrossEntropy - moments$Entropy,
			        SummaryExactKLScaled=mean(moments$CrossEntropy - moments$Entropy),
			        IndividualKL=KL, 
			        SummaryKL=mean(KL),
			        IndividualKLScaled=KLScaled,
			        SummaryKLScaled=mean(KLScaled),
			        IndividualIC=IC, 
			        SummaryIC=mean(IC),
			        IndividualICScaled=ICScaled,
			        SummaryICScaled=mean(ICScaled),
			        ICtransform=ICtransform, 
			        ICtransformScaled=ICtransformScaled,
			        IndividualINICall=ICtransformScaled,
			        SummaryINICall=min(ICtransformScaled),
			        varzx=varzx, 
			        varzxScaled=varzxScaled,
			        lambdaMeanNormalized=lambdaMeanNormalized,
			        lambdaMedianNormalized=lambdaMedianNormalized,
			        lambda1Mean=lambda1Mean,
			        lambda1Median=lambda1Median,
			        SNR=SNR,
			        cnINIMean=cnINIMean,
			        cnINIMedian=cnINIMedian,
			        cnINIMeanSignal=cnINIMeanSignal,
			        cnINIMedianSignal=cnINIMedianSignal,
			        L_z=as.numeric(L_c),
			        L_z_Normalized=as.numeric(L_c_Normalized),
			        intensity=as.numeric(express), 
			        L_c_Linear=L_c_Linear,
			        L_c_Linear_Normalized=L_c_Linear_Normalized,
			        expressLinear=expressLinear,
			        L_c_Square=L_c_Square,
			        L_c_Square_Normalized=L_c_Square_Normalized,
			        expressSquare=expressSquare,
			        L_c_Softmax=L_c_Softmax,
			        L_c_Softmax_Normalized=L_c_Softmax_Normalized,
			        expressSoftmax=expressSoftmax,
			        L_c_Mean=L_c_Mean,
			        L_c_Mean_Normalized=L_c_Mean_Normalized,
			        expressMeanBC=expressMeanBC,
			        expressMedianBC=expressMedianBC,
			        L_c_Median=L_c_Median,
			        L_c_Median_Normalized=L_c_Median_Normalized,
			        expressBCMean=expressBCMean,
			        expressBCMedian=expressBCMedian))
}