PGLMM reml

PGLMM.reml.R

######################################
## Created by Matthew Helmus
######################################


PGLMM.reml<-function(sp)
{
  #global tH tinvW tVV tX

  sp<-abs(Re(sp))
  Cd<-matrix(0,dim(tVV[[1]])[1],dim(tVV[[1]])[2])
  for(i in 1:length(sp))
  {
	 Cd=Cd + sp[i] * tVV[[i]]
  }
  V<- tinvW + Cd
  invV<- solve(V,diag(x=1,nrow=dim(V)[1],ncol=dim(V)[2]))

  # check to see if V is positive definite
  if(all( eigen(V)$values >0 ))
  {
    cholV<-chol(V)
    # ML
    # LL=.5*(2*sum(log(diag(chol(V))))+tH'*invV*tH)

    #REML
    LL=.5 * (2 * sum(log(diag(cholV))) + t(tH) %*% invV %*% tH + log(det(t(tX) %*% invV %*% tX)))
  } else {
    LL<-10^10
  }
  LL
}