eb.ran.beta.mix_
function(ymat, xmat, ind.r, beta, bv, ev, mispat)
{
	# class number specified in GGMM
	n.class <- length(bv)
	# sample size 
	n.obs <- nrow(ymat)
	# number of repeated measures
	repn <- ncol(ymat)
	# create a list of indices for fixed effect 
	len.beta <- length(beta[[1]])
	betaind <- c(1:length(beta[[1]]))
	ind.f <- ind.r
	for(j in 1:n.class) {
		ind.f[[j]] <- setdiff(betaind, ind.r[[j]])
	}
	# create lists of objects for placing class-specific empirical Bayes (eb) estimates at level 1 &
	# the associated covariance matrices 
	eb.beta <- as.list(c(1:n.class))
	for(j in 1:n.class) {
		eb.beta[[j]] <- matrix(0, n.obs, len.beta)/0
		bv[[j]] <- as.matrix(bv[[j]])
	}
	# compute class-specific eb estimates for standardized level 1 residuals 
	for(i in 1:n.obs) {
		# indices for nonmissing components for subject i
		s <- sort.list(mispat[i,  ])
		s <- s[1:(repn - sum(mispat[i,  ]))]
		l.s <- length(s)
		# note that the following is computed for subjects with sufficient number of nonmissing records 
		if(l.s > (len.beta - 1)) {
			for(j in 1:n.class) {
				u1 <- ginverse(t(xmat[s, ind.r[[j]]]) %*% xmat[s, ind.r[[
					j]]]) %*% t(xmat[s, ind.r[[j]]])
				newv <- ginverse(u1 %*% ev[[j]][s, s] %*% t(u1)) + 
					ginverse(bv[[j]][ind.r[[j]], ind.r[[j]]])
				inv.newv <- ginverse(newv)
				if(length(ind.r[[j]]) > 1 && length(ind.f[[j]]) > 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] - xmat[s, ind.f[[j]]] %*% beta[[
						j]][ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + t(
						beta[[j]][ind.r[[j]]]) %*% ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) %*% inv.newv
				}
				# 6 cases are considered to ensure compatability among the matrix operations 
				if(length(ind.r[[j]]) > 1 && length(ind.f[[j]]) < 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s])) %*% ginverse(u1 %*% ev[[j]][s,
						s] %*% t(u1)) %*% inv.newv + t(beta[[j]][
						ind.r[[j]]]) %*% ginverse(bv[[j]][ind.r[[
						j]], ind.r[[j]]]) %*% inv.newv
				}
				if(length(ind.r[[j]]) > 1 && length(ind.f[[j]]) == 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% c(ymat[
						i, s] - xmat[s, ind.f[[j]]] * beta[[j]][
						ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + t(
						beta[[j]][ind.r[[j]]]) %*% ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) %*% inv.newv
				}
				if(length(ind.r[[j]]) < 2 && length(ind.f[[j]]) <1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] )) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + (
						beta[[j]][ind.r[[j]]]) * ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) * inv.newv
				}
             if(length(ind.r[[j]]) < 2 && length(ind.f[[j]]) ==1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] - xmat[s, ind.f[[j]]] * beta[[j]][
						ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + (
						beta[[j]][ind.r[[j]]]) * ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) * inv.newv
				}
				if(length(ind.r[[j]]) < 2 && length(ind.f[[j]]) > 1) {
					eb.beta[[j]][i, ind.r[[j]]] <- t(u1 %*% (ymat[
						i, s] - xmat[s, ind.f[[j]]] %*% beta[[
						j]][ind.f[[j]]])) %*% ginverse(u1 %*% ev[[
						j]][s, s] %*% t(u1)) %*% inv.newv + (
						beta[[j]][ind.r[[j]]]) * ginverse(bv[[
						j]][ind.r[[j]], ind.r[[j]]]) * inv.newv
				}
				
			}
		}
	}
	eb.beta
}