##############################并############################ Classical EM algorithm for Linear Mixed Model ##################并##########################并############并 em.mixed <-function(y,x,z,beta,var0,var1,maxiter=2000,tolerance =1e-0010) { time <-proc.time() n <-nrow(y) q1 <-nrow(z) conv <-1 LO <-loglike(y,x,z,beta,var0,var1) i<-0 cat("Iter. sigma0 sigmal Likelihood",fill=T) repeat if(i>maxiter){conv<-0 break} V <-c(var1)*z %*t(z)+c(var0)diag(n) Vinv <-solve(V) xb <-x %*beta resid <-(y-xb) temp1 <-Vinv *resid s0 <-c(var0)2 t(temp1)%*%temp1 c(var0)n -c(var0)2 tr(Vinv) s1 <-c(var1)2 t(temp1)%*%z%*'%t(z)%*%temp1+c(var1)*q1- c(var1)2 *tr(t(z)%*%Vinv%*%z) w <-xb c(var0)temp1 var0 <-s0/n var1 <-s1/q1 beta <-ginverse(t(x)%*%x)%*%t(x)%*%w L1 <-loglike(y,x,z,beta,var0,var1) if(L1 LO){print("log-likelihood must increase,llikel <llike0,break.") conv <-0 break i<-i+1 cat(""i,""varo,"",var1,",L1,fill=T) if(abs(L1 -LO)<tolerance){break}#check for convergence L0<-L1 list(beta=beta,var0=var0,var1=var1,Loglikelihood=LO) 21########################################################### # Classical EM algorithm for Linear Mixed Model # ########################################################### em.mixed <- function(y, x, z, beta, var0, var1,maxiter=2000,tolerance = 1e-0010) { time <-proc.time() n <- nrow(y) q1 <- nrow(z) conv <- 1 L0 <- loglike(y, x, z, beta, var0, var1) i<-0 cat(" Iter. sigma0 sigma1 Likelihood",fill=T) repeat { if(i>maxiter) {conv<-0 break} V <- c(var1) * z %*% t(z) + c(var0) * diag(n) Vinv <- solve(V) xb <- x %*% beta resid <- (y-xb) temp1 <- Vinv %*% resid s0 <- c(var0)^2 * t(temp1)%*%temp1 + c(var0) * n - c(var0)^2 * tr(Vinv) s1 <- c(var1)^2 * t(temp1)%*%z%*%t(z)%*%temp1+ c(var1)*q1 - c(var1)^2 *tr(t(z)%*%Vinv%*%z) w <- xb + c(var0) * temp1 var0 <- s0/n var1 <- s1/q1 beta <- ginverse( t(x) %*% x) %*% t(x)%*% w L1 <- loglike(y, x, z, beta, var0, var1) if(L1 < L0) { print("log-likelihood must increase, llikel <llikeO, break.") conv <- 0 break } i <- i + 1 cat(" ", i," ",var0," ",var1," ",L1,fill=T) if(abs(L1 - L0) < tolerance) {break} #check for convergence L0 <- L1 } list(beta=beta, var0=var0,var1=var1,Loglikelihood=L0) } 21