ma.resid.2.mix_function(ymat, xmat, ind.r, beta, eb.beta, bv, ev, misp)
{
#input##############################################################	
# ymat: outcome matrix & ymat[i,] for ith subject                  #
# xmat: common time covariates s.t. xmat*beta=y                    #
# beta: the list of coefficients estimates for each class          #
# eb.beta: the list of eb estimates for latent variables           #
# ind.r: the list of column indices corresponding to random effects#
# bv  : the list of covariance matrix for latent variables         #
# ev  : the list of within subject covariance matrix               #
#mispat: missing pattern ( 1 for missing)                          # 
####################################################################

## output ##############################################################
#ma2: EB estimates of level 2 Mahalanobis (Ma) residuals w.r.t. classes#
########################################################################

#*************************** Main ***********************************#
	# class number specified in GGMM
	n.class <- length(bv)
	
	# sample size 
	n.obs <- nrow(ymat)

	# the number of growth factors
	len.beta <- ncol(xmat)
				
	# create the list of marginal covariance matrices &
	# the list of level 2 ma residual objects 
	betaind<-c(1:len.beta)
	newv <- ma2 <- vm <- ev
	for(j in 1:n.class) {
		vm[[j]] <- ev[[j]] + xmat[, ind.r[[j]]] %*% as.matrix(bv[[j]][ind.r[[j]],
			ind.r[[j]]]) %*% t(xmat[, ind.r[[j]]])
		ma2[[j]] <- rep(0, n.obs)/0
	}
	
	# compute eb estimates of level 2 Ma residuals  
	if(sum(misp) == 0) {
		for(j in 1:n.class) {
			len.beta <- length(beta[[j]])
			u1 <- ginverse(t(xmat[, ind.r[[j]]]) %*% xmat[, ind.r[[j]]]) %*%
				t(xmat[, ind.r[[j]]])
			inv1 <- ginverse(u1 %*% ev[[j]] %*% t(u1))
			vbb <- inv1 + ginverse(as.matrix(bv[[j]][ind.r[[j]], ind.r[[j]]]))
			invb <- ginverse(vbb)
			newv <- ((invb %*% (inv1 %*% u1)) %*% (vm[[j]]) %*% t(invb %*% (
				inv1 %*% u1)))
			ma2[[j]] <- diag(t(t(eb.beta[[j]][, ind.r[[j]]]) - beta[[j]][
				ind.r[[j]]]) %*% ginverse(newv) %*% (t(eb.beta[[j]][,
				ind.r[[j]]]) - beta[[j]][ind.r[[j]]]))
		}
	}
	if(sum(misp) > 0) {
		for(j in 1:n.class) {
			len.beta <- length(beta[[j]])
			for(i in 1:n.obs) {
				s <- 1:ncol(ymat)
				s <- s[misp[i,  ] == 0]
				l.s <- length(s)
				if((l.s > (len.beta - 1)) & (len.beta > 1)) {
					u1 <- ginverse(t(xmat[s, ind.r[[j]]]) %*% xmat[
						s, ind.r[[j]]]) %*% t(xmat[s, ind.r[[j]]])
					inv1 <- ginverse(u1 %*% ev[[j]][s, s] %*% t(u1))
					vbb <- inv1 + ginverse(as.matrix(bv[[j]][ind.r[[
						j]], ind.r[[j]]]))
					invb <- ginverse(vbb)
					newv <- ((invb %*% (inv1 %*% u1)) %*% (vm[[j]][
						s, s]) %*% t(invb %*% (inv1 %*% u1)))
					ma2[[j]][i] <- c(t(eb.beta[[j]][i, ind.r[[
						j]]] - beta[[j]][ind.r[[j]]]) %*% ginverse(
						newv) %*% (eb.beta[[j]][i, ind.r[[j]]] -
						beta[[j]][ind.r[[j]]]))
				}
				if((l.s > (len.beta - 1)) & (len.beta == 1)) {
					u1 <- ginverse(t(xmat[s, ind.r[[j]]]) %*% xmat[
						s, ind.r[[j]]]) %*% t(xmat[s, ind.r[[j]]])
					inv1 <- ginverse(u1 %*% ev[[j]][s, s] %*% t(u1))
					vbb <- inv1 + ginverse(as.matrix(bv[[j]][ind.r[[
						j]], ind.r[[j]]]))
					invb <- ginverse(vbb)
					newv <- ((invb %*% (inv1 %*% u1)) %*% (vm[[j]][
						s, s]) %*% t(invb %*% (inv1 %*% u1)))
					ma2[[j]][i] <- ((eb.beta[[j]][i, ind.r[[j]]] -
						beta[[j]][ind.r[[j]]])/newv * (eb.beta[[
						j]][i, ind.r[[j]]] - beta[[j]][ind.r[[
						j]]]))
				}
			}
		}
	}
	ma2
}