R/jEM.R

In syspremed/exploBATCH: Explore batch effect (exploBATCH) is a package for discovering and correcting for batch effect using an approach in Nyamundanda et al (2017).

jEM <-
function(n,Y,idx,nj,p,Covars,Alphaopt,Wopt,Sigopt,Iq,repV,eps,Np,Nvp,C2p,Ip,picon,V){
  Covarsj=Covars[,-idx[n]]
  J=Y[-idx[n],]
  muhat=colMeans(J)
  tJc=t(J)-muhat
  Jc=t(tJc)
  sumdiSj=mdiag(fM(tJc,Jc)/nj,p)
  Wj=Wopt
  Sigj=Sigopt
  iM=chol.inv(crossprod(Wj)+Sigj*Iq)
  iMtW=tcrossprod(iM,Wj)
  Alphaj=Alphaopt
  delj=Alphaj%*%Covarsj
  Xproj=crossprod(Covarsj,chol.inv(tcrossprod(Covarsj)))
  ll=lla=repV
  tol=eps+1; v=0
  while(tol>eps){
    v=v+1
    u=fM(iMtW,tJc)+iM%*%(Sigj*delj)
    tu=t(u)
    sumEuu=nj*Sigj*iM+u%*%tu
    Alphaj=u%*%Xproj
    delj=Alphaj%*%Covarsj
    Wj=(tJc%*%tu)%*%chol.inv(sumEuu)
    ww=crossprod(Wj)
    MLESigj=(nj*sumdiSj+sum(ww*sumEuu)-2*sum(fM(Jc,Wj)*tu))/(nj*p)
    Sigj=(Np*MLESigj+C2p)/Nvp
    SIGj=Sigj*Iq
    M=ww+SIGj
    iM=chol.inv(M)
    iMtW=tcrossprod(iM,Wj)
    InvS=(Ip-Wj%*%iMtW)/Sigj
    tJcc=tJc-Wj%*%delj
    MhalDist=colSums(tJcc*fM(InvS,tJcc))
    logdetD=p*log(Sigj)+log(det(M/Sigj))
    ll[v]=sum(-picon-.5*logdetD-.5*MhalDist)
    converge=Aitkenf(ll,lla,v,eps)
    tol=converge[[1]]
    lla[v]=converge[[2]]
    if(v==V){tol=eps-1}}
  Alphaj}