ma.resid.1.mix_function(ymat, xmat, ind.r, beta, eb.beta, bv, ev, mispat)
{
#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 ##############################################################
#ma1: EB estimates of level 1 Mahalanobis (Ma) residuals w.r.t. classes#
########################################################################

	# 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 
	n.class <- length(bv)
	ma <- vm <- ev
	for(j in 1:n.class) {
		bv[[j]] <- as.matrix(bv[[j]][ind.r[[j]], ind.r[[j]]])
		vm[[j]] <- ev[[j]] + as.matrix(xmat[, ind.r[[j]]]) %*% bv[[j]] %*% t(
			as.matrix(xmat[, ind.r[[j]]]))
	}
	
	# create the list of fixed effect indices
	betaind <- c(1:len.beta)
	ind.f <- ind.r
	for(j in 1:n.class) {
		ind.f[[j]] <- setdiff(betaind, ind.r[[j]])
	}
	
	# create the list of level 1 Ma residual objects
	ma1 <-res1 <- as.list(c(1:n.class))
	for(j in 1:n.class) {
		res1[[j]] <- matrix(0, n.obs, ncol(ymat))/0
		ma1[[j]] <- rep(0, n.obs)/0
	}
	
	# compute level 1 Ma residual estimates
	for(i in 1:n.obs) {
		s <- sort.list(mispat[i,  ])
		s <- s[1:(ncol(ymat) - sum(mispat[i,  ]))]
		l.s <- length(s)
		if(l.s > (len.beta - 1)) {
			for(j in 1:n.class) {
				v.e1 <- ev[[j]][s, s] %*% ginverse(vm[[j]][s, s]) %*% ev[[
					j]][s, s]
				if(length(ind.f[[j]]) < 1 && length(ind.r[[j]]) > 1) {
					res1[[j]][i, s] <- ( - xmat[s, ind.r[[j]]] %*% 
						eb.beta[[j]][i, ind.r[[j]]] + ymat[i,
						s])
				}
				if(length(ind.f[[j]]) == 1 && length(ind.r[[j]]) > 1) {
					res1[[j]][i, s] <- ( - xmat[s, ind.f[[j]]] * beta[[
						j]][ind.f[[j]]] - xmat[s, ind.r[[j]]] %*%
						eb.beta[[j]][i, ind.r[[j]]] + ymat[i,
						s])
				}
				if(length(ind.f[[j]]) > 1 && length(ind.r[[j]]) > 1) {
					res1[[j]][i, s] <- ( - xmat[s, ind.f[[j]]] %*% 
						beta[[j]][ind.f[[j]]] - xmat[s, ind.r[[
						j]]] %*% eb.beta[[j]][i, ind.r[[j]]] +
						ymat[i, s])
				}
				if(length(ind.f[[j]]) > 1 && length(ind.r[[j]]) == 1) {
					res1[[j]][i, s] <- ( - xmat[s, ind.f[[j]]] %*% 
						beta[[j]][ind.f[[j]]] - xmat[s, ind.r[[
						j]]] * eb.beta[[j]][i, ind.r[[j]]] + ymat[
						i, s])
				}
				if(length(ind.f[[j]]) < 1 && length(ind.r[[j]]) == 1) {
					res1[[j]][i, s] <- ( - xmat[s, ind.r[[j]]] * 
						eb.beta[[j]][i, ind.r[[j]]] + ymat[i,
						s])
				}
				if(length(ind.f[[j]]) == 1 && length(ind.r[[j]]) == 1) {
					res1[[j]][i, s] <- ( - xmat[s, ind.f[[j]]] * beta[[
						j]][ind.f[[j]]] - xmat[s, ind.r[[j]]] *
						eb.beta[[j]][i, ind.r[[j]]] + ymat[i,
						s])
				}
				ma1[[j]][i] <- t(res1[[j]][i, s]) %*% ginverse(v.e1) %*% 
					res1[[j]][i, s]
			}
		}
	}
	ma1
#*********************************** End ****************************************#	

}