sd.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 ##########################################################
#res1: EB estimates of normalized level 1 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 for marginal covariance matrices 
	n.class <- length(bv)
	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 1 residual objects
	res1 <- as.list(c(1:n.class))
	for(j in 1:n.class) {
		res1[[j]] <- matrix(0, nrow(ymat), ncol(ymat))/0
	}
	
	# compute class-specific eb estimates of standardized level 1 residual 
	# v.e1: covariance of the eb estimates of level 1 residuals
	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.r[[j]]) < 2 && length(ind.f[[j]]) == 0) {
					res1[[j]][i, s] <- ( - xmat[s, ind.r[[j]]] * 
						eb.beta[[j]][i, ind.r[[j]]] + ymat[i,
						s])/diag(v.e1)^0.5
				}
				if(length(ind.r[[j]]) < 2 && length(ind.f[[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])/diag(v.e1)^0.5
				}
				if(length(ind.r[[j]]) < 2 && length(ind.f[[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])/diag(v.e1)^0.5
				}
				if(length(ind.r[[j]]) > 1 && length(ind.f[[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])/diag(v.e1)^0.5
				}
				if(length(ind.r[[j]]) > 1 && length(ind.f[[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])/diag(v.e1)^0.5
				}
				if(length(ind.r[[j]]) > 1 && length(ind.f[[j]]) == 0) {
					res1[[j]][i, s] <- ( - xmat[s, ind.r[[j]]] %*% 
						eb.beta[[j]][i, ind.r[[j]]] + ymat[i,
						s])/diag(v.e1)^0.5
				}
			}
		}
	}
	res1
#*********************************** End ****************************************#	
}