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 <pauljoh...@gmail.com> wrote:
> On Tue, Jul 3, 2012 at 12:41 PM, jimmycloud <jimmycl...@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@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@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.
