sd.resid.2.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 ##########################################################
#res2: EB estimates of normalized level 2 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 
	vm <- ev
	for(j in 1:n.class) {
		bv[[j]] <- as.matrix(bv[[j]][ind.r[[j]], ind.r[[j]]])
		vm[[j]] <- ev[[j]] + xmat[, ind.r[[j]]] %*% bv[[j]] %*% t(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 2 residual objects
	res2 <- as.list(c(1:n.class))
	for(j in 1:n.class) {
		res2[[j]] <- matrix(0, n.obs, len.beta)/0
	}
	
   # compute class-specific eb estimates for standardized level 2 residuals 
	for(i in 1:n.obs) {
		# indices for nonmissing components for subject i
       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) {
				
				# covariance of the eb estimates of level 2 residuals
				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(bv[[j]])
				invb <- ginverse(vbb)
				newv <- ((invb %*% (inv1 %*% u1)) %*% (vm[[j]][s, s]) %*%
					t(invb %*% (inv1 %*% u1)))
					
				# standardization	
				res2[[j]][i, ind.r[[j]]] <- (eb.beta[[j]][i, ind.r[[
					j]]] - beta[[j]][ind.r[[j]]])/(diag(newv)^0.5)
			}
		}
	}
	res2
#*********************************** End ****************************************#	
}