R help - Jul 2012 - EM algorithm to find MLE of coeff in mixed effects model

I have a general question about coefficients estimation of the mixed model.

I simulated a very basic model: Y|b=X*\beta+Z*b +\sigma^2* diag(ni);
                                                             b follows
N(0,\psi)  #i.e. bivariate normal
where b is the latent variable, Z and X are ni*2 design matrices, sigma is 
the error variance,
Y are longitudinal data, i.e. there are ni measurements for object i.
Parameters are \beta, \sigma, \psi; call them \theta.

I wrote a regular EM, the M step is to maximize the log(f(Y,b;\theta))  by
taking first order derivatives, setting to 0 and solving the equation.
The E step involves the evaluation of E step, using Gauss Hermite 
approximation. (10 points)

All are simulated data. X and Z are naive like cbind(rep(1,m),1:m)
After 200 iterations, the estimated \beta converges to the true value while 
\sigma and \psi do not. Even after one iteration, the \sigma goes up from
about 10^0 to about 10^1.
I am confused since the \hat{\beta} requires \sigma and \psi from previous 
iteration. If something wrong then all estimations should be incorrect...

Another question is that I calculated the logf(Y;\theta) to see if it 
increases after updating \theta.
Seems decreasing.....

I thought the X and Z are linearly dependent would cause some issue but I 
also changed the X and Z to some random numbers from normal distribution.

I also tried ECM, which converges fast but sigma and psi are not close to
the desired values.
Is this due to the limitation of some methods that I used?

Can any one give some help? I am stuck for a week. I can send the code to
you. 
First time to raise a question here. Not sure if it is proper to post all
code.
##########################################################################
# the main R script
n=100
beta=c(-0.5,1)
vvar=2   #sigma^2=2
psi=matrix(c(1,0.2,0.2,1),2,2)
b.m.true=mvrnorm(n=n,mu=c(0,0),Sigma=psi)  #100*2 matrix, each row is the
b_i
Xi=cbind(rnorm(7,mean=2,sd=0.5),log(2:8)) #Xi=cbind(rep(1,7),1:7)
y.m=matrix(NA,nrow=n,ncol=nrow(Xi))   #100*7, each row is a y_i
Zi=Xi

b.list=as.list(data.frame(t(b.m.true)))
psi.old=matrix(c(0.5,0.4,0.4,0.5),2,2)      #starting psi, beta and var, not
exactly the same as the true value
beta.old=c(-0.3,0.7)
var.old=1.7

gausspar=gauss.quad(10,kind="hermite",alpha=0,beta=0)

data.list.wob=list()
for (i in 1:n)
{
data.list.wob[[i]]=list(Xi=Xi,yi=y.m[i,],Zi=Zi)
}

#compute true loglikelihood and initial loglikelihood
truelog=0
for (i in 1:length(data.list.wob))
{
truelog=truelog+loglike(data.list.wob[[i]],vvar,beta,psi)
}

truelog

loglikeseq=c()
loglikeseq[1]=sum(sapply(data.list.wob,loglike))

ECM=F


for (m in 1:300)
{

Sig.hat=Zi%*%psi.old%*%t(Zi)+var.old*diag(nrow(Zi))
W.hat=psi.old-psi.old%*%t(Zi)%*%solve(Sig.hat)%*%Zi%*%psi.old

Sigb=psi.old-psi.old%*%t(Zi)%*%solve(Sig.hat)%*%Zi%*%psi.old
det(Sigb)^(-0.5)

Y.minus.X.beta=t(t(y.m)-c(Xi%*%beta.old))
miu.m=t(apply(Y.minus.X.beta,MARGIN=1,function(s,B=psi.old%*%t(Zi)%*%solve(Sig.hat))
						{
						B%*%s
						}
		))  ### each row is the miu_i


tmp1=permutations(length(gausspar$nodes),2,repeats.allowed=T)
tmp2=c(tmp1)
a.mat=matrix(gausspar$nodes[tmp2],nrow=nrow(tmp1)) #a1,a1
								   #a1,a2
								   #...
								   #a10,a9
								   #a10,a10
a.mat.list=as.list(data.frame(t(a.mat)))
tmp1=permutations(length(gausspar$weights),2,repeats.allowed=T)
tmp2=c(tmp1)
weights.mat=matrix(gausspar$weights[tmp2],nrow=nrow(tmp1)) #w1,w1
								   #w1,w2
								   #...
								   #w10,w9
								   #w10,w10
L=chol(solve(W.hat))
LL=sqrt(2)*solve(L)
halfb=t(LL%*%t(a.mat))

# each page of b.array is all values of bi_k and bi_j for b_i
b.list=list()
for (i in 1:n)
{
b.list[[i]]=t(t(halfb)+miu.m[i,])
}

#generate a list, each page contains Xi,yi,Zi,
data.list=list()
for (i in 1:n)
{
data.list[[i]]=list(Xi=Xi,yi=y.m[i,],Zi=Zi,b=b.list[[i]])
}

#update sigma^2
t1=proc.time()
tempaa=c()
tempbb=c()
for (j in 1:n)
{
#tempaa[j]=Eh4new(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
tempbb[j]=Eh4newv2(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)

}
var.new=mean(tempbb)
if (ECM==T){var.old=var.new}

sumXiXi=matrix(rowSums(sapply(data.list,function(s){t(s$Xi)%*%(s$Xi)})),ncol=ncol(Xi))
tempbb=c()
for (j in 1:n)
{
tempbb=cbind(tempbb,Eh2new(data.list[[j]],weights.m=weights.mat))
}
beta.new=solve(sumXiXi)%*%rowSums(tempbb)

if (ECM==T){beta.old=beta.new}

#update psi
tempcc=array(NA,c(2,2,n))
sumtempcc=matrix(0,2,2)
for (j in 1:n)
{
tempcc[,,j]=Eh1new(data.list[[j]],weights.m=weights.mat)
sumtempcc=sumtempcc+tempcc[,,j]
}
psi.new=sumtempcc/n

#stop
if(sum(abs(beta.old-beta.new))+sum(abs(psi.old-psi.new))+sum(abs(var.old-var.new))<0.01)
{print("converge, stop");break;}

#update 
var.old=var.new
psi.old=psi.new
beta.old=beta.new
loglikeseq[m+1]=sum(sapply(data.list,loglike))
} # end of M loop

########################################################################
#function to calculate loglikelihood
loglike=function(datai=data.list[[i]],vvar=var.old,beta=beta.old,psi=psi.old)
{
temp1=-0.5*log(det(vvar*diag(nrow(datai$Zi))+datai$Zi%*%psi%*%t(datai$Zi)))
temp2=-0.5*t(datai$yi-datai$Xi%*%beta)%*%solve(vvar*diag(nrow(datai$Zi))+datai$Zi%*%psi%*%t(datai$Zi))%*%(datai$yi-datai$Xi%*%beta)
return(temp1+temp2)
}

#######################################################################
#functions to evaluate the conditional expection, saved as Efun v2.R
#Eh1new=E(bibi')
#Eh2new=E(X'(y-Zbi))
#Eh3new=E(bi'Z'Zbi)
#Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
#Eh5new=E(Xibeta-yi)'Zibi
#Eh6new=E(bi)

Eh1new=function(datai=data.list[[i]],
		 weights.m=weights.mat)
{
#one way
#temp=matrix(0,2,2)
#for (i in 1:nrow(bi))
#{
#temp=temp+bi[i,]%*%t(bi[i,])*weights.m[i,1]*weights.m[i,2]
#}
#print(temp)

#the other way
bi=datai$b
tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b, so need
sqrt
#deno=sum(weights.m[,1]*weights.m[,2])

return(t(tempb)%*%tempb/pi)
} # end of Eh1


#Eh2new=E(X'(y-Zbi))
Eh2new=function(datai=data.list[[i]],
		 weights.m=weights.mat)
{
#one way
#temp=rep(0,2)
#for (j in 1:nrow(bi))
#{
#temp=temp+c(t(datai$Xi)%*%(datai$yi-datai$Zi%*%bi[j,])*weights.m[j,1]*weights.m[j,2])
#}
#deno=sum(weights.m[,1]*weights.m[,2])
#print(temp/deno)

#another way
bi=datai$b
tempb=bi*rep(weights.m[,1]*weights.m[,2],2)

tt=t(datai$Xi)%*%datai$yi-t(datai$Xi)%*%datai$Zi%*%colSums(tempb)/pi
return(c(tt))
} # end of Eh2


#Eh3new=E(b'Z'Zbi)
Eh3new=function(datai=data.list[[i]],
		 weights.m=weights.mat)
{
#one way
#deno=sum(weights.m[,1]*weights.m[,2])
#tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b, so need
sqrt
#sum(apply(datai$Zi%*%t(tempb),MARGIN=2,function(s){sum(s^2)}))/deno

#another way
bi=datai$b
tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b, so need
sqrt
return(sum(sapply(as.list(data.frame(datai$Zi%*%t(tempb))),function(s){sum(s^2)}))/pi)
}  # end of Eh3

#Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
Eh4new=function(datai=data.list[[i]],
		 weights.m=weights.mat,beta=beta.old)
{
#one way
#temp=0
#bi=datai$b
#tt=c()
#for (j in 1:nrow(bi))
#{
#tt[j]=sum((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[j,])^2)*weights.m[j,1]*weights.m[j,2]
#temp=temp+sum((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[j,])^2)*weights.m[j,1]*weights.m[j,2]
#}
#temp/deno

# another way
bi=datai$b

tt=sapply(as.list(ata.frame(c(datai$yi-datai$Xi%*%beta)-datai$Zi%*%t(bi))),
	    function(s){sum(s^2)})*(weights.m[,1]*weights.m[,2])
return(sum(tt)/pi)
} # end of Eh4


Eh4newv2=function(datai=data.list[[i]],
		 weights.m=weights.mat,beta=beta.old)
{
bi=datai$b
v=weights.m[,1]*weights.m[,2]
temp=c()
for (i in 1:length(v))
{
temp[i]=sum(((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[i,])*sqrt(v[i]))^2)
}
return(sum(temp)/pi)
} # end of Eh4


#Eh5new=E(Xibeta-yi)'Zib
Eh5new=function(datai=data.list[[i]],
		 weights.m=weights.mat,beta=beta.old)
{
bi=datai$b
tempb=bi*rep(weights.m[,1]*weights.m[,2],2)
t(datai$Xi%*%beta-datai$yi)%*%(datai$Zi%*%t(tempb))
return(sum(t(datai$Xi%*%beta-datai$yi)%*%(datai$Zi%*%t(tempb)))/pi)
}



Eh6new=function(datai=data.list[[i]],
		 weights.m=weights.mat)
{
bi=datai$b
tempw=weights.m[,1]*weights.m[,2]
for (i in 1:nrow(bi))
{
bi[i,]=bi[i,]*tempw[i]
}
return(colMeans(bi)/pi)
} # end of Eh1








--
View this message in context:
http://r.789695.n4.nabble.com/EM-algorithm-to-find-MLE-of-coeff-in-mixed-effects-model-tp4635315.html
Sent from the R help mailing list archive at Nabble.com.

On Tue, Jul 3, 2012 at 12:41 PM, jimmycloud <jimmycloud at gmail.com>
wrote:
> I have a general question about coefficients estimation of the mixed model.
>
I have 2 ideas for you.

1. Fit with lme4 package, using the lmer function. That's what it is for.

2. If you really want to write your own EM algorithm, I don't feel
sure that very many R and EM experts are going to want to read through
the code you have because you don't follow some of the minimal
readability guidelines.  I accept the fact that there is no officially
mandated R style guide, except for "indent with 4 spaces, not tabs."
But there are some almost universally accepted basics like

1. Use <-, not =, for assignment
2. put a space before and  after symbols like <-, = , + , * , and so forth.
3. I'd suggest you get used to putting squiggly braces in the K&R style.

I have found the formatR package's tool tidy.source can do this
nicely. From tidy.source, here's what I get with your code

http://pj.freefaculty.org/R/em2.R

Much more readable!
(I inserted library(MASS) for you at the top, otherwise this doesn't
even survive the syntax check.)

When I copy from email to Emacs, some line-breaking monkey-business
occurs, but I expect you get the idea.

Now, looking at your cleaned up code, I can see some things to tidy up
to improve the chances that some expert will look.

First, R functions don't require "return" at the end, many experts
consider it distracting. (Run "lm" or such and review how they write.
No return at the end of functions).

Second, about that big for loop at the top, the one that goes from m 1:300

I don't know what that loop is doing, but there's some code in it that
I'm suspicious of. For example, this part:

 W.hat <- psi.old - psi.old %*% t(Zi) %*% solve(Sig.hat) %*%  Zi %*% psi.old

    Sigb <- psi.old - psi.old %*% t(Zi) %*% solve(Sig.hat) %*%  Zi %*%
psi.old


First, you've caused the possibly slow calculation of solve
 (Sig.hat) to run two times.  If you really need it, run it once and
save the result.


Second, a for loop is not necessarily slow, but it may be easier to
read if you re-write this:

 for (j in 1:n) {
tempaa[j]=Eh4new(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
        tempbb[j] <- Eh4newv2(datai = data.list[[j]], weights.m weights.mat,
beta = beta.old)
 }

like this:

tempaa <- lapply(data.list, Eh4new, weights.m, beta.old)
tempbb <- lapply(data.list, Eh4newv2, weights.m, beta.old)

Third, here is a no-no

tempbb <- c()
for (j in 1:n) {
        tempbb <- cbind(tempbb, Eh2new(data.list[[j]], weights.m =
weights.mat))
    }

That will call cbind over and over, causing a relatively slow memory
re-allocation.  See
(http://pj.freefaculty.org/R/WorkingExamples/stackListItems.R)

Instead, do this to apply Eh2new to each item in data.list

tempbbtemp <- lapply(data.list, Eh2new, weights.mat)

and then smash the results together in one go

tempbb <- do.call("cbind", tempbbtemp)


Fourth, I'm not sure on the matrix algebra. Are you sure you need
solve to get the full inverse of Sig.hat?  Once you start checking
into how estimates are calculated in R, you find that the
paper-and-pencil matrix algebra style of formula is generally frowned
upon. OLS (in lm.fit) doesn't do solve(X'X), it uses a QR
decomposition of matrices.  OR look in MASS package ridge regression
code, where the SVD is used to get the inverse.

I wish I knew enough about the EM algorithm to gaze at your code and
say "aha, error in line 332", but I'm not.  But if you clean up
the
presentation and tighten up the most obvious things, you improve your
chances that somebody who is an expert will exert him/herself to do
it.

pj

>                                                              b follows
> N(0,\psi)  #i.e. bivariate normal
> where b is the latent variable, Z and X are ni*2 design matrices, sigma is
> the error variance,
> Y are longitudinal data, i.e. there are ni measurements for object i.
> Parameters are \beta, \sigma, \psi; call them \theta.
>
> I wrote a regular EM, the M step is to maximize the log(f(Y,b;\theta))  by
> taking first order derivatives, setting to 0 and solving the equation.
> The E step involves the evaluation of E step, using Gauss Hermite
> approximation. (10 points)
>
> All are simulated data. X and Z are naive like cbind(rep(1,m),1:m)
> After 200 iterations, the estimated \beta converges to the true value while
> \sigma and \psi do not. Even after one iteration, the \sigma goes up from
> about 10^0 to about 10^1.
> I am confused since the \hat{\beta} requires \sigma and \psi from previous
> iteration. If something wrong then all estimations should be incorrect...
>
> Another question is that I calculated the logf(Y;\theta) to see if it
> increases after updating \theta.
> Seems decreasing.....
>
> I thought the X and Z are linearly dependent would cause some issue but I
> also changed the X and Z to some random numbers from normal distribution.
>
> I also tried ECM, which converges fast but sigma and psi are not close to
> the desired values.
> Is this due to the limitation of some methods that I used?
>
> Can any one give some help? I am stuck for a week. I can send the code to
> you.
> First time to raise a question here. Not sure if it is proper to post all
> code.
> ##########################################################################
> # the main R script
> n=100
> beta=c(-0.5,1)
> vvar=2   #sigma^2=2
> psi=matrix(c(1,0.2,0.2,1),2,2)
> b.m.true=mvrnorm(n=n,mu=c(0,0),Sigma=psi)  #100*2 matrix, each row is the
> b_i
> Xi=cbind(rnorm(7,mean=2,sd=0.5),log(2:8)) #Xi=cbind(rep(1,7),1:7)
> y.m=matrix(NA,nrow=n,ncol=nrow(Xi))   #100*7, each row is a y_i
> Zi=Xi
>
> b.list=as.list(data.frame(t(b.m.true)))
> psi.old=matrix(c(0.5,0.4,0.4,0.5),2,2)      #starting psi, beta and var,
not
> exactly the same as the true value
> beta.old=c(-0.3,0.7)
> var.old=1.7
>
> gausspar=gauss.quad(10,kind="hermite",alpha=0,beta=0)
>
> data.list.wob=list()
> for (i in 1:n)
> {
> data.list.wob[[i]]=list(Xi=Xi,yi=y.m[i,],Zi=Zi)
> }
>
> #compute true loglikelihood and initial loglikelihood
> truelog=0
> for (i in 1:length(data.list.wob))
> {
> truelog=truelog+loglike(data.list.wob[[i]],vvar,beta,psi)
> }
>
> truelog
>
> loglikeseq=c()
> loglikeseq[1]=sum(sapply(data.list.wob,loglike))
>
> ECM=F
>
>
> for (m in 1:300)
> {
>
> Sig.hat=Zi%*%psi.old%*%t(Zi)+var.old*diag(nrow(Zi))
> W.hat=psi.old-psi.old%*%t(Zi)%*%solve(Sig.hat)%*%Zi%*%psi.old
>
> Sigb=psi.old-psi.old%*%t(Zi)%*%solve(Sig.hat)%*%Zi%*%psi.old
> det(Sigb)^(-0.5)
>
> Y.minus.X.beta=t(t(y.m)-c(Xi%*%beta.old))
>
miu.m=t(apply(Y.minus.X.beta,MARGIN=1,function(s,B=psi.old%*%t(Zi)%*%solve(Sig.hat))
>                                                 {
>                                                 B%*%s
>                                                 }
>                 ))  ### each row is the miu_i
>
>
> tmp1=permutations(length(gausspar$nodes),2,repeats.allowed=T)
> tmp2=c(tmp1)
> a.mat=matrix(gausspar$nodes[tmp2],nrow=nrow(tmp1)) #a1,a1
>                                                                    #a1,a2
>                                                                    #...
>                                                                    #a10,a9
>                                                                    #a10,a10
> a.mat.list=as.list(data.frame(t(a.mat)))
> tmp1=permutations(length(gausspar$weights),2,repeats.allowed=T)
> tmp2=c(tmp1)
> weights.mat=matrix(gausspar$weights[tmp2],nrow=nrow(tmp1)) #w1,w1
>                                                                    #w1,w2
>                                                                    #...
>                                                                    #w10,w9
>                                                                    #w10,w10
> L=chol(solve(W.hat))
> LL=sqrt(2)*solve(L)
> halfb=t(LL%*%t(a.mat))
>
> # each page of b.array is all values of bi_k and bi_j for b_i
> b.list=list()
> for (i in 1:n)
> {
> b.list[[i]]=t(t(halfb)+miu.m[i,])
> }
>
> #generate a list, each page contains Xi,yi,Zi,
> data.list=list()
> for (i in 1:n)
> {
> data.list[[i]]=list(Xi=Xi,yi=y.m[i,],Zi=Zi,b=b.list[[i]])
> }
>
> #update sigma^2
> t1=proc.time()
> tempaa=c()
> tempbb=c()
> for (j in 1:n)
> {
> #tempaa[j]=Eh4new(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
>
tempbb[j]=Eh4newv2(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
>
> }
> var.new=mean(tempbb)
> if (ECM==T){var.old=var.new}
>
>
sumXiXi=matrix(rowSums(sapply(data.list,function(s){t(s$Xi)%*%(s$Xi)})),ncol=ncol(Xi))
> tempbb=c()
> for (j in 1:n)
> {
> tempbb=cbind(tempbb,Eh2new(data.list[[j]],weights.m=weights.mat))
> }
> beta.new=solve(sumXiXi)%*%rowSums(tempbb)
>
> if (ECM==T){beta.old=beta.new}
>
> #update psi
> tempcc=array(NA,c(2,2,n))
> sumtempcc=matrix(0,2,2)
> for (j in 1:n)
> {
> tempcc[,,j]=Eh1new(data.list[[j]],weights.m=weights.mat)
> sumtempcc=sumtempcc+tempcc[,,j]
> }
> psi.new=sumtempcc/n
>
> #stop
>
if(sum(abs(beta.old-beta.new))+sum(abs(psi.old-psi.new))+sum(abs(var.old-var.new))<0.01)
> {print("converge, stop");break;}
>
> #update
> var.old=var.new
> psi.old=psi.new
> beta.old=beta.new
> loglikeseq[m+1]=sum(sapply(data.list,loglike))
> } # end of M loop
>
> ########################################################################
> #function to calculate loglikelihood
>
loglike=function(datai=data.list[[i]],vvar=var.old,beta=beta.old,psi=psi.old)
> {
> temp1=-0.5*log(det(vvar*diag(nrow(datai$Zi))+datai$Zi%*%psi%*%t(datai$Zi)))
>
temp2=-0.5*t(datai$yi-datai$Xi%*%beta)%*%solve(vvar*diag(nrow(datai$Zi))+datai$Zi%*%psi%*%t(datai$Zi))%*%(datai$yi-datai$Xi%*%beta)
> return(temp1+temp2)
> }
>
> #######################################################################
> #functions to evaluate the conditional expection, saved as Efun v2.R
> #Eh1new=E(bibi')
> #Eh2new=E(X'(y-Zbi))
> #Eh3new=E(bi'Z'Zbi)
> #Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
> #Eh5new=E(Xibeta-yi)'Zibi
> #Eh6new=E(bi)
>
> Eh1new=function(datai=data.list[[i]],
>                  weights.m=weights.mat)
> {
> #one way
> #temp=matrix(0,2,2)
> #for (i in 1:nrow(bi))
> #{
> #temp=temp+bi[i,]%*%t(bi[i,])*weights.m[i,1]*weights.m[i,2]
> #}
> #print(temp)
>
> #the other way
> bi=datai$b
> tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b, so need
> sqrt
> #deno=sum(weights.m[,1]*weights.m[,2])
>
> return(t(tempb)%*%tempb/pi)
> } # end of Eh1
>
>
> #Eh2new=E(X'(y-Zbi))
> Eh2new=function(datai=data.list[[i]],
>                  weights.m=weights.mat)
> {
> #one way
> #temp=rep(0,2)
> #for (j in 1:nrow(bi))
> #{
>
#temp=temp+c(t(datai$Xi)%*%(datai$yi-datai$Zi%*%bi[j,])*weights.m[j,1]*weights.m[j,2])
> #}
> #deno=sum(weights.m[,1]*weights.m[,2])
> #print(temp/deno)
>
> #another way
> bi=datai$b
> tempb=bi*rep(weights.m[,1]*weights.m[,2],2)
>
> tt=t(datai$Xi)%*%datai$yi-t(datai$Xi)%*%datai$Zi%*%colSums(tempb)/pi
> return(c(tt))
> } # end of Eh2
>
>
> #Eh3new=E(b'Z'Zbi)
> Eh3new=function(datai=data.list[[i]],
>                  weights.m=weights.mat)
> {
> #one way
> #deno=sum(weights.m[,1]*weights.m[,2])
> #tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b, so need
> sqrt
> #sum(apply(datai$Zi%*%t(tempb),MARGIN=2,function(s){sum(s^2)}))/deno
>
> #another way
> bi=datai$b
> tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b, so need
> sqrt
>
return(sum(sapply(as.list(data.frame(datai$Zi%*%t(tempb))),function(s){sum(s^2)}))/pi)
> }  # end of Eh3
>
> #Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
> Eh4new=function(datai=data.list[[i]],
>                  weights.m=weights.mat,beta=beta.old)
> {
> #one way
> #temp=0
> #bi=datai$b
> #tt=c()
> #for (j in 1:nrow(bi))
> #{
>
#tt[j]=sum((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[j,])^2)*weights.m[j,1]*weights.m[j,2]
>
#temp=temp+sum((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[j,])^2)*weights.m[j,1]*weights.m[j,2]
> #}
> #temp/deno
>
> # another way
> bi=datai$b
>
> tt=sapply(as.list(ata.frame(c(datai$yi-datai$Xi%*%beta)-datai$Zi%*%t(bi))),
>             function(s){sum(s^2)})*(weights.m[,1]*weights.m[,2])
> return(sum(tt)/pi)
> } # end of Eh4
>
>
> Eh4newv2=function(datai=data.list[[i]],
>                  weights.m=weights.mat,beta=beta.old)
> {
> bi=datai$b
> v=weights.m[,1]*weights.m[,2]
> temp=c()
> for (i in 1:length(v))
> {
> temp[i]=sum(((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[i,])*sqrt(v[i]))^2)
> }
> return(sum(temp)/pi)
> } # end of Eh4
>
>
> #Eh5new=E(Xibeta-yi)'Zib
> Eh5new=function(datai=data.list[[i]],
>                  weights.m=weights.mat,beta=beta.old)
> {
> bi=datai$b
> tempb=bi*rep(weights.m[,1]*weights.m[,2],2)
> t(datai$Xi%*%beta-datai$yi)%*%(datai$Zi%*%t(tempb))
> return(sum(t(datai$Xi%*%beta-datai$yi)%*%(datai$Zi%*%t(tempb)))/pi)
> }
>
>
>
> Eh6new=function(datai=data.list[[i]],
>                  weights.m=weights.mat)
> {
> bi=datai$b
> tempw=weights.m[,1]*weights.m[,2]
> for (i in 1:nrow(bi))
> {
> bi[i,]=bi[i,]*tempw[i]
> }
> return(colMeans(bi)/pi)
> } # end of Eh1
>
>
>
>
>
>
>
>
> --
> View this message in context:
http://r.789695.n4.nabble.com/EM-algorithm-to-find-MLE-of-coeff-in-mixed-effects-model-tp4635315.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu

I don't believe that R-help permits pdf files.  A useful workaround is to
post it to a file hosting site like MediaFire and post the link here.
John Kane
Kingston ON Canada

> -----Original Message-----
> From: jimmycloud at gmail.com
> Sent: Wed, 4 Jul 2012 22:56:02 -0400
> To: pauljohn32 at gmail.com
> Subject: Re: [R] EM algorithm to find MLE of coeff in mixed effects model
> 
> Dear Paul,
> 
> Thank you for your suggestion. I was moved by the fact that people are so
> nice to help learners and ask for nothing.
> With your help, I made some changes to my code and add more comments, try
> to make things more clear.
> If this R email list allow me to upload a pdf file to illustrate the
> formula, it will be great.
> The reason I do not use lme4 is that later I plan to use this for a more
> complicated model (Y,T), Y is the same response of mixed model and T is a
> drop out process.
> (Frankly I did it for the more complicated model first and got NaN
> hessian
> after one iteration, so I turned to the simple model.)
> In the code, the m loop is the EM iterations. About efficiency, thank you
> again for pointing it out. I compared sapply and for loop, they are tie
> and
> for loop is even better sometimes.
> In a paper by Bates, he mentioned that the asymptotic properties of beta
> is
> kind of independent of the other two parameters. But it seems that paper
> was submitted but I did not find it on google scholar.
> Not sure if it is the reason for my results. That is why I hope some
> expert
> who have done some simulation similar may be willing to give some
> suggestion.
> I may turn to other methods to calculate the conditional expection of the
> latent variable as the same time.
> 
> Modified code is as below:
> 
> # library(PerformanceAnalytics)
> # install.packages("gregmisc")
> # install.packages("statmod")
> library(gregmisc)
> library(MASS)
> library(statmod)
> 
> ## function to calculate loglikelihood
> loglike <- function(datai = data.list[[i]], vvar = var.old,
>     beta = beta.old, psi = psi.old) {
>     temp1 <- -0.5 * log(det(vvar * diag(nrow(datai$Zi)) + datai$Zi %*%
>         psi %*% t(datai$Zi)))
>     temp2 <- -0.5 * t(datai$yi - datai$Xi %*% beta) %*% solve(vvar *
>         diag(nrow(datai$Zi)) + datai$Zi %*% psi %*% t(datai$Zi)) %*%
>         (datai$yi - datai$Xi %*% beta)
>     temp1 + temp2
> }
> 
> ## functions to evaluate the conditional expection, saved as Efun v2.R
> ## Eh1new=E(bibi')
> ## Eh2new=E(X'(y-Zbi))
> ## Eh3new=E(bi'Z'Zbi)
> ## Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
> ## Eh5new=E(Xibeta-yi)'Zibi
> ## Eh6new=E(bi)
> 
> Eh1new <- function(datai = data.list[[i]], weights.m = weights.mat) {
>     bi <- datai$b
>     tempb <- bi * rep(sqrt(weights.m[, 1] * weights.m[, 2]),
>         2)  #quadratic b, so need    sqrt
>     t(tempb) %*% tempb/pi
> }  # end of Eh1
> 
> 
> # Eh2new=E(X'(y-Zbi))
> Eh2new <- function(datai = data.list[[i]], weights.m = weights.mat) {
>     bi <- datai$b
>     tempb <- bi * rep(weights.m[, 1] * weights.m[, 2], 2)
>     tt <- t(datai$Xi) %*% datai$yi - t(datai$Xi) %*% datai$Zi %*%
>         colSums(tempb)/pi
>     c(tt)
> }  # end of Eh2
> 
> 
> # Eh3new=E(b'Z'Zbi)
> Eh3new <- function(datai = data.list[[i]], weights.m = weights.mat) {
>     bi <- datai$b
>     tempb <- bi * rep(sqrt(weights.m[, 1] * weights.m[, 2]),
>         2)  #quadratic b, so need sqrt
>     sum(sapply(as.list(data.frame(datai$Zi %*% t(tempb))),
>         function(s) {
>             sum(s^2)
>         }))/pi
> }  # end of Eh3
> 
> # Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
> Eh4new <- function(datai = data.list[[i]], weights.m = weights.mat,
>     beta = beta.old) {
>     bi <- datai$b
>     tt <- sapply(as.list(ata.frame(c(datai$yi - datai$Xi %*%
>         beta) - datai$Zi %*% t(bi))), function(s) {
>         sum(s^2)
>     }) * (weights.m[, 1] * weights.m[, 2])
>     sum(tt)/pi
> }  # end of Eh4
> 
> Eh4newv2 <- function(datai = data.list[[i]], weights.m = weights.mat,
>     beta = beta.old) {
>     bi <- datai$b
>     v <- weights.m[, 1] * weights.m[, 2]
>     temp <- c()
>     for (i in 1:length(v)) {
>         temp[i] <- sum(((datai$yi - datai$Xi %*% beta - datai$Zi %*%
>             bi[i, ]) * sqrt(v[i]))^2)
>     }
>     sum(temp)/pi
> }  # end of Eh4
> 
> # Eh5new=E(Xibeta-yi)'Zib
> Eh5new <- function(datai = data.list[[i]], weights.m = weights.mat,
>     beta = beta.old) {
>     bi <- datai$b
>     tempb <- bi * rep(weights.m[, 1] * weights.m[, 2], 2)
>     t(datai$Xi %*% beta - datai$yi) %*% (datai$Zi %*% t(tempb))
>     sum(t(datai$Xi %*% beta - datai$yi) %*% (datai$Zi %*%
>         t(tempb)))/pi
> }
> 
> Eh6new <- function(datai = data.list[[i]], weights.m = weights.mat) {
>     bi <- datai$b
>     tempw <- weights.m[, 1] * weights.m[, 2]
>     for (i in 1:nrow(bi)) {
>         bi[i, ] <- bi[i, ] * tempw[i]
>     }
>     colMeans(bi)/pi
> }  # end of Eh6
> 
> ## the main R script
> ################### initial the values and generate the data
> ##################
> n <- 100                                               #number of
> observations
> beta <- c(-0.5, 1)                               #true coefficient of
> fixed
> effects
> vvar <- 2  #sigma^2=2                                  #true error
> variance
> of epsilon
> psi <- matrix(c(1, 0.2, 0.2, 1), 2, 2)                 #covariance
matrix
> of random effects b
> b.m.true <- mvrnorm(n = n, mu = c(0, 0), Sigma = psi)  #100*2 matrix,
> each
> row is the b_i
> Xi <- cbind(rnorm(7, mean = 2, sd = 0.5), log(2:8))
> Zi <- Xi
> y.m <- matrix(NA, nrow = n, ncol = nrow(Xi))            #100*7, each row
> is
> a y_i Zi=Xi
> 
> for( i in 1:n)
> {
> y.m[i,]=mvrnorm(1,mu=(Xi%*%beta+Zi%*%b.m.true[i,]),vvar*diag(nrow(Xi)))
> }
> b.list <- as.list(data.frame(t(b.m.true)))
> psi.old <- matrix(c(0.5, 0.4, 0.4, 0.5), 2, 2)          #starting psi,
> beta
> and var,not the true values
> beta.old <- c(-0.3, 0.7)
> var.old <- 1.7
> 
> gausspar <- gauss.quad(10, kind = "hermite", alpha = 0, beta =
0)
> 
> # data.list.wob contains X,Y and Z, but no b's
> data.list.wob <- list()
> for (i in 1:n) {
>     data.list.wob[[i]] <- list(Xi = Xi, yi = y.m[i, ], Zi = Zi)
> }
> 
> # compute true loglikelihood and initial loglikelihood
> truelog <- 0
> for (i in 1:length(data.list.wob)) {
>     truelog <- truelog + loglike(data.list.wob[[i]], vvar, beta, psi)
> }
> truelog
> 
> loglikeseq <- c()
> loglikeseq[1] <- sum(sapply(data.list.wob, loglike))
> 
> ECM <- F
> 
> ############################# E-M steps
> ################################################
> # m loop is the EM iteration
> for (m in 1:300) {
>     Sig.hat <- Zi %*% psi.old %*% t(Zi) + var.old * diag(nrow(Zi))
> #estimated covariance matrix of y
>     Sig.hat.inv <- solve(Sig.hat)
>     W.hat <- psi.old - psi.old %*% t(Zi) %*% Sig.hat.inv %*%
>         Zi %*% psi.old
>     Sigb <- psi.old - psi.old %*% t(Zi) %*% Sig.hat.inv %*% Zi %*%
> psi.old
>  #estimated covariance matrix of b
> 
>     Y.minus.X.beta <- t(t(y.m) - c(Xi %*% beta.old))
>     miu.m <- t(apply(Y.minus.X.beta, MARGIN = 1, function(s,
>         B = psi.old %*% t(Zi) %*% solve(Sig.hat)) {
>         B %*% s
>     }))  ### each row is the miu_i
> 
>     tmp1 <- permutations(length(gausspar$nodes), 2, repeats.allowed = T)
>     tmp2 <- c(tmp1)
>     a.mat <- matrix(gausspar$nodes[tmp2], nrow = nrow(tmp1))  # a1,a1
>      # a1,a2
>   # ...
>   # a10,a9
>   # a10,a10
>   # a.mat are all ordered pairs of gauss hermite nodes
>     a.mat.list <- as.list(data.frame(t(a.mat)))
>     tmp1 <- permutations(length(gausspar$weights), 2, repeats.allowed
> T)
>     tmp2 <- c(tmp1)
>     weights.mat <- matrix(gausspar$weights[tmp2], nrow = nrow(tmp1))  #
> w1,w1
>        # w1,w2
>       # ...
>       # w10,w9
>       # w10,w10
>       # weights.mat are all ordered pairs of gauss hermite weights
>     L <- chol(solve(W.hat))
>     LL <- sqrt(2) * solve(L)
>     halfb <- t(LL %*% t(a.mat))
>     # each page of b.array is all values of bi_k and bi_j for b_i
>     # need to calculate those b's as linear functions of nodes since
the
> original integral needs transformation to use gauss approximation
>     b.list <- list()
>     for (i in 1:n) {
>         b.list[[i]] <- t(t(halfb) + miu.m[i, ])
>     }
> 
>     # generate a list, each page contains Xi,yi,Zi, and b's
>     data.list <- list()
>     for (i in 1:n) {
>         data.list[[i]] <- list(Xi = Xi, yi = y.m[i, ], Zi = Zi,
>             b = b.list[[i]])
>     }
> 
>     # update sigma^2
>     t1 <- proc.time()
>     tempaa <- c()
>     tempbb <- c()
>     for (j in 1:n) {
>         tempbb[j] <- Eh4newv2(datai = data.list[[j]], weights.m >
weights.mat, beta = beta.old)
>     }
>     var.new <- mean(tempbb)
>     if (ECM == T) {
>         var.old <- var.new
>     }
> 
>     sumXiXi <- matrix(rowSums(sapply(data.list, function(s) {
>         t(s$Xi) %*% (s$Xi)
>     })), ncol = ncol(Xi))
>     tempbb=sapply(data.list,Eh2new)
>     beta.new <- solve(sumXiXi) %*% rowSums(tempbb)
> 
>     if (ECM == T) {
>         beta.old <- beta.new
>     }
> 
>     # update psi
>     tempcc <- array(NA, c(2, 2, n))
>     sumtempcc <- matrix(0, 2, 2)
>     for (j in 1:n) {
>         tempcc[, , j] <- Eh1new(data.list[[j]], weights.m = weights.mat)
>         sumtempcc <- sumtempcc + tempcc[, , j]
>     }
>     psi.new <- sumtempcc/n
> 
>     # stop condition
>     if (sum(abs(beta.old - beta.new)) + sum(abs(psi.old - psi.new)) +
>         sum(abs(var.old - var.new)) < 0.01) {
>         print("converge, stop")
>         break
>     }
> 
>     # update parameters
>     var.old <- var.new
>     psi.old <- psi.new
>     beta.old <- beta.new
>     loglikeseq[m + 1] <- sum(sapply(data.list, loglike))
> }  # end of M loop
> 
> 
> On Tue, Jul 3, 2012 at 5:06 PM, Paul Johnson <pauljohn32 at
gmail.com>
> wrote:
> 
>> On Tue, Jul 3, 2012 at 12:41 PM, jimmycloud <jimmycloud at
gmail.com>
>> wrote:
>>> I have a general question about coefficients estimation of the
mixed
>> model.
>>> 
>> 
>> I have 2 ideas for you.
>> 
>> 1. Fit with lme4 package, using the lmer function. That's what it
is
>> for.
>> 
>> 2. If you really want to write your own EM algorithm, I don't feel
>> sure that very many R and EM experts are going to want to read through
>> the code you have because you don't follow some of the minimal
>> readability guidelines.  I accept the fact that there is no officially
>> mandated R style guide, except for "indent with 4 spaces, not
tabs."
>> But there are some almost universally accepted basics like
>> 
>> 1. Use <-, not =, for assignment
>> 2. put a space before and  after symbols like <-, = , + , * , and so
>> forth.
>> 3. I'd suggest you get used to putting squiggly braces in the
K&R style.
>> 
>> I have found the formatR package's tool tidy.source can do this
>> nicely. From tidy.source, here's what I get with your code
>> 
>> http://pj.freefaculty.org/R/em2.R
>> 
>> Much more readable!
>> (I inserted library(MASS) for you at the top, otherwise this
doesn't
>> even survive the syntax check.)
>> 
>> When I copy from email to Emacs, some line-breaking monkey-business
>> occurs, but I expect you get the idea.
>> 
>> Now, looking at your cleaned up code, I can see some things to tidy up
>> to improve the chances that some expert will look.
>> 
>> First, R functions don't require "return" at the end,
many experts
>> consider it distracting. (Run "lm" or such and review how
they write.
>> No return at the end of functions).
>> 
>> Second, about that big for loop at the top, the one that goes from m
>> 1:300
>> 
>> I don't know what that loop is doing, but there's some code in
it that
>> I'm suspicious of. For example, this part:
>> 
>>  W.hat <- psi.old - psi.old %*% t(Zi) %*% solve(Sig.hat) %*%  Zi %*%
>> psi.old
>> 
>>     Sigb <- psi.old - psi.old %*% t(Zi) %*% solve(Sig.hat) %*%  Zi
%*%
>> psi.old
>> 
>> 
>> First, you've caused the possibly slow calculation of solve
>>  (Sig.hat) to run two times.  If you really need it, run it once and
>> save the result.
>> 
>> 
>> Second, a for loop is not necessarily slow, but it may be easier to
>> read if you re-write this:
>> 
>>  for (j in 1:n) {
>>
tempaa[j]=Eh4new(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
>>         tempbb[j] <- Eh4newv2(datai = data.list[[j]], weights.m
>> weights.mat, beta = beta.old)
>>  }
>> 
>> like this:
>> 
>> tempaa <- lapply(data.list, Eh4new, weights.m, beta.old)
>> tempbb <- lapply(data.list, Eh4newv2, weights.m, beta.old)
>> 
>> Third, here is a no-no
>> 
>> tempbb <- c()
>> for (j in 1:n) {
>>         tempbb <- cbind(tempbb, Eh2new(data.list[[j]], weights.m
>> weights.mat))
>>     }
>> 
>> That will call cbind over and over, causing a relatively slow memory
>> re-allocation.  See
>> (http://pj.freefaculty.org/R/WorkingExamples/stackListItems.R)
>> 
>> Instead, do this to apply Eh2new to each item in data.list
>> 
>> tempbbtemp <- lapply(data.list, Eh2new, weights.mat)
>> 
>> and then smash the results together in one go
>> 
>> tempbb <- do.call("cbind", tempbbtemp)
>> 
>> 
>> Fourth, I'm not sure on the matrix algebra. Are you sure you need
>> solve to get the full inverse of Sig.hat?  Once you start checking
>> into how estimates are calculated in R, you find that the
>> paper-and-pencil matrix algebra style of formula is generally frowned
>> upon. OLS (in lm.fit) doesn't do solve(X'X), it uses a QR
>> decomposition of matrices.  OR look in MASS package ridge regression
>> code, where the SVD is used to get the inverse.
>> 
>> I wish I knew enough about the EM algorithm to gaze at your code and
>> say "aha, error in line 332", but I'm not.  But if you
clean up the
>> presentation and tighten up the most obvious things, you improve your
>> chances that somebody who is an expert will exert him/herself to do
>> it.
>> 
>> pj
>> 
>> 
>>>                                                              b
follows
>>> N(0,\psi)  #i.e. bivariate normal
>>> where b is the latent variable, Z and X are ni*2 design matrices,
sigma
>> is
>>> the error variance,
>>> Y are longitudinal data, i.e. there are ni measurements for object
i.
>>> Parameters are \beta, \sigma, \psi; call them \theta.
>>> 
>>> I wrote a regular EM, the M step is to maximize the
log(f(Y,b;\theta))
>>  by
>>> taking first order derivatives, setting to 0 and solving the
equation.
>>> The E step involves the evaluation of E step, using Gauss Hermite
>>> approximation. (10 points)
>>> 
>>> All are simulated data. X and Z are naive like cbind(rep(1,m),1:m)
>>> After 200 iterations, the estimated \beta converges to the true
value
>> while
>>> \sigma and \psi do not. Even after one iteration, the \sigma goes
up
>>> from
>>> about 10^0 to about 10^1.
>>> I am confused since the \hat{\beta} requires \sigma and \psi from
>> previous
>>> iteration. If something wrong then all estimations should be
>>> incorrect...
>>> 
>>> Another question is that I calculated the logf(Y;\theta) to see if
it
>>> increases after updating \theta.
>>> Seems decreasing.....
>>> 
>>> I thought the X and Z are linearly dependent would cause some issue
but
>>> I
>>> also changed the X and Z to some random numbers from normal
>>> distribution.
>>> 
>>> I also tried ECM, which converges fast but sigma and psi are not
close
>>> to
>>> the desired values.
>>> Is this due to the limitation of some methods that I used?
>>> 
>>> Can any one give some help? I am stuck for a week. I can send the
code
>>> to
>>> you.
>>> First time to raise a question here. Not sure if it is proper to
post
>>> all
>>> code.
>>> 
>>
##########################################################################
>>> # the main R script
>>> n=100
>>> beta=c(-0.5,1)
>>> vvar=2   #sigma^2=2
>>> psi=matrix(c(1,0.2,0.2,1),2,2)
>>> b.m.true=mvrnorm(n=n,mu=c(0,0),Sigma=psi)  #100*2 matrix, each row
is
>>> the
>>> b_i
>>> Xi=cbind(rnorm(7,mean=2,sd=0.5),log(2:8)) #Xi=cbind(rep(1,7),1:7)
>>> y.m=matrix(NA,nrow=n,ncol=nrow(Xi))   #100*7, each row is a y_i
>>> Zi=Xi
>>> 
>>> b.list=as.list(data.frame(t(b.m.true)))
>>> psi.old=matrix(c(0.5,0.4,0.4,0.5),2,2)      #starting psi, beta and
>>> var,
>> not
>>> exactly the same as the true value
>>> beta.old=c(-0.3,0.7)
>>> var.old=1.7
>>> 
>>> gausspar=gauss.quad(10,kind="hermite",alpha=0,beta=0)
>>> 
>>> data.list.wob=list()
>>> for (i in 1:n)
>>> {
>>> data.list.wob[[i]]=list(Xi=Xi,yi=y.m[i,],Zi=Zi)
>>> }
>>> 
>>> #compute true loglikelihood and initial loglikelihood
>>> truelog=0
>>> for (i in 1:length(data.list.wob))
>>> {
>>> truelog=truelog+loglike(data.list.wob[[i]],vvar,beta,psi)
>>> }
>>> 
>>> truelog
>>> 
>>> loglikeseq=c()
>>> loglikeseq[1]=sum(sapply(data.list.wob,loglike))
>>> 
>>> ECM=F
>>> 
>>> 
>>> for (m in 1:300)
>>> {
>>> 
>>> Sig.hat=Zi%*%psi.old%*%t(Zi)+var.old*diag(nrow(Zi))
>>> W.hat=psi.old-psi.old%*%t(Zi)%*%solve(Sig.hat)%*%Zi%*%psi.old
>>> 
>>> Sigb=psi.old-psi.old%*%t(Zi)%*%solve(Sig.hat)%*%Zi%*%psi.old
>>> det(Sigb)^(-0.5)
>>> 
>>> Y.minus.X.beta=t(t(y.m)-c(Xi%*%beta.old))
>>> 
>>
miu.m=t(apply(Y.minus.X.beta,MARGIN=1,function(s,B=psi.old%*%t(Zi)%*%solve(Sig.hat))
>>>                                                 {
>>>                                                 B%*%s
>>>                                                 }
>>>                 ))  ### each row is the miu_i
>>> 
>>> 
>>> tmp1=permutations(length(gausspar$nodes),2,repeats.allowed=T)
>>> tmp2=c(tmp1)
>>> a.mat=matrix(gausspar$nodes[tmp2],nrow=nrow(tmp1)) #a1,a1
>>> 
>>> #a1,a2
>>>                                                                   
#...
>>> 
>>  #a10,a9
>>> 
>>  #a10,a10
>>> a.mat.list=as.list(data.frame(t(a.mat)))
>>> tmp1=permutations(length(gausspar$weights),2,repeats.allowed=T)
>>> tmp2=c(tmp1)
>>> weights.mat=matrix(gausspar$weights[tmp2],nrow=nrow(tmp1)) #w1,w1
>>> 
>>> #w1,w2
>>>                                                                   
#...
>>> 
>>  #w10,w9
>>> 
>>  #w10,w10
>>> L=chol(solve(W.hat))
>>> LL=sqrt(2)*solve(L)
>>> halfb=t(LL%*%t(a.mat))
>>> 
>>> # each page of b.array is all values of bi_k and bi_j for b_i
>>> b.list=list()
>>> for (i in 1:n)
>>> {
>>> b.list[[i]]=t(t(halfb)+miu.m[i,])
>>> }
>>> 
>>> #generate a list, each page contains Xi,yi,Zi,
>>> data.list=list()
>>> for (i in 1:n)
>>> {
>>> data.list[[i]]=list(Xi=Xi,yi=y.m[i,],Zi=Zi,b=b.list[[i]])
>>> }
>>> 
>>> #update sigma^2
>>> t1=proc.time()
>>> tempaa=c()
>>> tempbb=c()
>>> for (j in 1:n)
>>> {
>>> 
>>
#tempaa[j]=Eh4new(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
>>> 
>>
tempbb[j]=Eh4newv2(datai=data.list[[j]],weights.m=weights.mat,beta=beta.old)
>>> 
>>> }
>>> var.new=mean(tempbb)
>>> if (ECM==T){var.old=var.new}
>>> 
>>> 
>>
sumXiXi=matrix(rowSums(sapply(data.list,function(s){t(s$Xi)%*%(s$Xi)})),ncol=ncol(Xi))
>>> tempbb=c()
>>> for (j in 1:n)
>>> {
>>> tempbb=cbind(tempbb,Eh2new(data.list[[j]],weights.m=weights.mat))
>>> }
>>> beta.new=solve(sumXiXi)%*%rowSums(tempbb)
>>> 
>>> if (ECM==T){beta.old=beta.new}
>>> 
>>> #update psi
>>> tempcc=array(NA,c(2,2,n))
>>> sumtempcc=matrix(0,2,2)
>>> for (j in 1:n)
>>> {
>>> tempcc[,,j]=Eh1new(data.list[[j]],weights.m=weights.mat)
>>> sumtempcc=sumtempcc+tempcc[,,j]
>>> }
>>> psi.new=sumtempcc/n
>>> 
>>> #stop
>>> 
>>
if(sum(abs(beta.old-beta.new))+sum(abs(psi.old-psi.new))+sum(abs(var.old-var.new))<0.01)
>>> {print("converge, stop");break;}
>>> 
>>> #update
>>> var.old=var.new
>>> psi.old=psi.new
>>> beta.old=beta.new
>>> loglikeseq[m+1]=sum(sapply(data.list,loglike))
>>> } # end of M loop
>>> 
>>>
########################################################################
>>> #function to calculate loglikelihood
>>> 
>>
loglike=function(datai=data.list[[i]],vvar=var.old,beta=beta.old,psi=psi.old)
>>> {
>>> 
>>
temp1=-0.5*log(det(vvar*diag(nrow(datai$Zi))+datai$Zi%*%psi%*%t(datai$Zi)))
>>> 
>>
temp2=-0.5*t(datai$yi-datai$Xi%*%beta)%*%solve(vvar*diag(nrow(datai$Zi))+datai$Zi%*%psi%*%t(datai$Zi))%*%(datai$yi-datai$Xi%*%beta)
>>> return(temp1+temp2)
>>> }
>>> 
>>>
#######################################################################
>>> #functions to evaluate the conditional expection, saved as Efun
v2.R
>>> #Eh1new=E(bibi')
>>> #Eh2new=E(X'(y-Zbi))
>>> #Eh3new=E(bi'Z'Zbi)
>>> #Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
>>> #Eh5new=E(Xibeta-yi)'Zibi
>>> #Eh6new=E(bi)
>>> 
>>> Eh1new=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat)
>>> {
>>> #one way
>>> #temp=matrix(0,2,2)
>>> #for (i in 1:nrow(bi))
>>> #{
>>> #temp=temp+bi[i,]%*%t(bi[i,])*weights.m[i,1]*weights.m[i,2]
>>> #}
>>> #print(temp)
>>> 
>>> #the other way
>>> bi=datai$b
>>> tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b,
so
>>> need
>>> sqrt
>>> #deno=sum(weights.m[,1]*weights.m[,2])
>>> 
>>> return(t(tempb)%*%tempb/pi)
>>> } # end of Eh1
>>> 
>>> 
>>> #Eh2new=E(X'(y-Zbi))
>>> Eh2new=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat)
>>> {
>>> #one way
>>> #temp=rep(0,2)
>>> #for (j in 1:nrow(bi))
>>> #{
>>> 
>>
#temp=temp+c(t(datai$Xi)%*%(datai$yi-datai$Zi%*%bi[j,])*weights.m[j,1]*weights.m[j,2])
>>> #}
>>> #deno=sum(weights.m[,1]*weights.m[,2])
>>> #print(temp/deno)
>>> 
>>> #another way
>>> bi=datai$b
>>> tempb=bi*rep(weights.m[,1]*weights.m[,2],2)
>>> 
>>>
tt=t(datai$Xi)%*%datai$yi-t(datai$Xi)%*%datai$Zi%*%colSums(tempb)/pi
>>> return(c(tt))
>>> } # end of Eh2
>>> 
>>> 
>>> #Eh3new=E(b'Z'Zbi)
>>> Eh3new=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat)
>>> {
>>> #one way
>>> #deno=sum(weights.m[,1]*weights.m[,2])
>>> #tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b,
so
>> need
>>> sqrt
>>>
#sum(apply(datai$Zi%*%t(tempb),MARGIN=2,function(s){sum(s^2)}))/deno
>>> 
>>> #another way
>>> bi=datai$b
>>> tempb=bi*rep(sqrt(weights.m[,1]*weights.m[,2]),2)   #quadratic b,
so
>>> need
>>> sqrt
>>> 
>>
return(sum(sapply(as.list(data.frame(datai$Zi%*%t(tempb))),function(s){sum(s^2)}))/pi)
>>> }  # end of Eh3
>>> 
>>> #Eh4new=E(Y-Xibeta-Zibi)'(Y-Xibeta-Zibi)
>>> Eh4new=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat,beta=beta.old)
>>> {
>>> #one way
>>> #temp=0
>>> #bi=datai$b
>>> #tt=c()
>>> #for (j in 1:nrow(bi))
>>> #{
>>> 
>>
#tt[j]=sum((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[j,])^2)*weights.m[j,1]*weights.m[j,2]
>>> 
>>
#temp=temp+sum((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[j,])^2)*weights.m[j,1]*weights.m[j,2]
>>> #}
>>> #temp/deno
>>> 
>>> # another way
>>> bi=datai$b
>>> 
>>> 
>>
tt=sapply(as.list(ata.frame(c(datai$yi-datai$Xi%*%beta)-datai$Zi%*%t(bi))),
>>>             function(s){sum(s^2)})*(weights.m[,1]*weights.m[,2])
>>> return(sum(tt)/pi)
>>> } # end of Eh4
>>> 
>>> 
>>> Eh4newv2=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat,beta=beta.old)
>>> {
>>> bi=datai$b
>>> v=weights.m[,1]*weights.m[,2]
>>> temp=c()
>>> for (i in 1:length(v))
>>> {
>>>
temp[i]=sum(((datai$yi-datai$Xi%*%beta-datai$Zi%*%bi[i,])*sqrt(v[i]))^2)
>>> }
>>> return(sum(temp)/pi)
>>> } # end of Eh4
>>> 
>>> 
>>> #Eh5new=E(Xibeta-yi)'Zib
>>> Eh5new=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat,beta=beta.old)
>>> {
>>> bi=datai$b
>>> tempb=bi*rep(weights.m[,1]*weights.m[,2],2)
>>> t(datai$Xi%*%beta-datai$yi)%*%(datai$Zi%*%t(tempb))
>>> return(sum(t(datai$Xi%*%beta-datai$yi)%*%(datai$Zi%*%t(tempb)))/pi)
>>> }
>>> 
>>> 
>>> 
>>> Eh6new=function(datai=data.list[[i]],
>>>                  weights.m=weights.mat)
>>> {
>>> bi=datai$b
>>> tempw=weights.m[,1]*weights.m[,2]
>>> for (i in 1:nrow(bi))
>>> {
>>> bi[i,]=bi[i,]*tempw[i]
>>> }
>>> return(colMeans(bi)/pi)
>>> } # end of Eh1
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> --
>>> View this message in context:
>>
http://r.789695.n4.nabble.com/EM-algorithm-to-find-MLE-of-coeff-in-mixed-effects-model-tp4635315.html
>>> Sent from the R help mailing list archive at Nabble.com.
>>> 
>>> ______________________________________________
>>> R-help at r-project.org mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>> 
>> 
>> 
>> --
>> Paul E. Johnson
>> Professor, Political Science    Assoc. Director
>> 1541 Lilac Lane, Room 504     Center for Research Methods
>> University of Kansas               University of Kansas
>> http://pj.freefaculty.org            http://quant.ku.edu
>> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
____________________________________________________________
FREE ONLINE PHOTOSHARING - Share your photos online with your friends and
family!
Visit http://www.inbox.com/photosharing to find out more!