[R] EM algorithm to find MLE of coeff in mixed effects model
John Kane
jrkrideau at inbox.com
Thu Jul 5 15:28:26 CEST 2012
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!
More information about the R-help
mailing list