R/Mlmmc.R

Mlmmc <-
  function(y,fit,cov){
    n <- length(y)
    X <- model.matrix(fit,fit$data)
    groups <- as.character(fit$groups[[1]])
    group <- factor(groups,ordered=FALSE)
    n.group <- length(levels(group))
    MLE <- c(diag(getVarCov(fit)),fit$sigma^2)
    p <- length(MLE)
    if (is.null(cov)){cov <- fit$apVar}
    covi <- vector("list", n.group)
    rbind(covi)
    Zi <- vector("list", n.group)
    Zi <- split(group,group)
    yi <- vector("list",n.group)
    yi <- split(y,group)
    Xi <- vector("list",n.group)
    Xi <- split(as.data.frame(X),group)
    ri <- vector("list",n.group)
    ri1 <- vector("list",n.group)
    ri2 <- vector("list",n.group)
    lik1 <- 0
    lik2 <- 0
    D <- 0
    H.lik <- function(yi, Xi, MLE) {
      for (j in 1:n.group) {
        covi[[j]] <- as.matrix(as.numeric(Zi[[j]])) %*% as.matrix(MLE[1:p-1]) %*% t(as.matrix(as.numeric(Zi[[j]])))+diag(MLE[p],length(yi[[j]]),length(yi[[j]]))
        D[j] <- log(det(covi[[j]]))
        ri1[[j]] <- t(as.matrix(Xi[[j]])) %*% solve(covi[[j]]) %*% as.matrix(Xi[[j]])
        ri2[[j]] <- t(as.matrix(Xi[[j]])) %*% solve(covi[[j]]) %*% yi[[j]]
      }
      for (j in 1:n.group) {
        ri[[j]] <- yi[[j]]-as.matrix(Xi[[j]]) %*% (solve(Reduce("+",ri1)) %*% Reduce("+",ri2))
        lik1[j] <- t(ri[[j]]) %*% solve(covi[[j]]) %*% ri[[j]]
        lik2[j] <- log(det(t(as.matrix(Xi[[j]])) %*% solve(covi[[j]]) %*% as.matrix(Xi[[j]])))
      }
      lik <- -0.5*(n-p)*log(2*pi)-0.5*sum(D)-0.5*sum(lik1)-0.5*sum(lik2)
      lik
    }
    max_loglik <- fit$logLik
    input <- list(y=yi,X=Xi,MLE=MLE,max_loglik=max_loglik)
    H <- function(epsilon,input,VZ,cstar){
      H.lik(y=yi,X=Xi,MLE=as.numeric(MLE+epsilon*VZ))-max_loglik+0.5*cstar
    }
    list(cov=cov,input=input,H=H,H.lik=H.lik)
  }
ppham27/setsim documentation built on May 25, 2019, 11:25 a.m.