R/ML.r

In nshi-stat/netiim3: Improved Inference Methods for Network Meta-Analysis

#' Maximum likelihood estimation for a random effects multivariate meta-analysis
#'
#' Maximum likelihood estimation for a random effects model for multivariate meta-analysis.
#'
#' The correlation matrix of the between-studies covariance matrix is set to compound symmetry
#'with the correlation coefficient 0.50. It can be changed by modifying the source code.
#' Please see Noma et al. (2017) for details.
#'
#' @param y N x p matrix of outcome variables.
#' @param S Series of within-study covariance matrices of the outcome variables.
#' A matrix or data frame with N rows and p(p+1)/2 columns.
#' @param maxitr The maximum iteration number of the Newton-Raphson algorithm.
#' @return
#' \item{Coefficients}{The maximum likeliood estimate of the grand mean vector and its standard error (SE) with the 95\% confidence interval.}
#' \item{Between-studies_SD}{The maximum likeliood estimate of the between-studies standard deviance (SD).}
#' \item{Between-studies_COR}{The correlation coefficient of the between-studies correlation matrix.}
#' \item{Loglikelihood}{The loglikelihood at the converged point.}
#' @references
#' Jackson, D., Riley, R., White, I. R. (2011).
#' Multivariate meta-analysis: Potential and promise.
#' \emph{Statistics in Medicine}. \strong{30}: 2481-2498.
#'
#' Noma, H., Nagashima, K., Maruo, K., Gosho, M., Furukawa, T. A. (2017).
#' Bartlett-type corrections and bootstrap adjustments of likelihood-based inference methods for network meta-analysis.
#' \emph{ISM Research Memorandum} 1205.
#' @examples
#' # dae <- data.aug.edit(smoking)
#' # y <- dae$y
#' # S <- dae$S
#'
#' # ML(y, S)
#' @family ML
#' @export
ML <- function(y, S, maxitr = 200){

  # ML algorithm for the consistent model (Sigma is exchangeble)

  N <- dim(y)[1]
  p <- dim(y)[2]

  mu <- rnorm(p)	# initial values
  g1 <- 0.2
  g2 <- 0.1

  Qc0 <- c(mu, g1, g2)

  LL1 <- function(g){

    #G <- gmat(g, g2, p)
    G <- gmat(g, (g/2), p)

    ll1 <- 0

    for(i in 1:N){

      yi <- as.vector(y[i,])
      wi <- which(is.na(yi) == FALSE)
      pl <- length(wi)

      Si <- vmat(S[i,], p)

      yi <- yi[wi]
      Si <- pmat(Si, wi)
      mui <- mu[wi]
      Gi <- pmat(G, wi)

      B1 <- (yi - mui)
      B2 <- ginv(Gi + Si)

      A1 <- log(det(Gi + Si))
      A2 <- t(B1) %*% B2 %*% B1
      A3 <- pl * log(2*pi)

      ll1 <- ll1 + A1 + A2 + A3

    }

    return(ll1)

  }

  LL2 <- function(g){

    G <- gmat(g1, g, p)

    ll1 <- 0

    for(i in 1:N){

      yi <- as.vector(y[i,])
      wi <- which(is.na(yi) == FALSE)
      pl <- length(wi)

      Si <- vmat(S[i,], p)

      yi <- yi[wi]
      Si <- pmat(Si, wi)
      mui <- mu[wi]
      Gi <- pmat(G, wi)

      B1 <- (yi - mui)
      B2 <- ginv(Gi + Si)

      A1 <- log(det(Gi + Si))
      A2 <- t(B1) %*% B2 %*% B1
      A3 <- pl * log(2*pi)

      ll1 <- ll1 + A1 + A2 + A3

    }

    return(ll1)   # return "minus loglikelihood"

  }

  for(itr in 1:maxitr){

    A1 <- numeric(p)
    A2 <- matrix(numeric(p*p), p)

    G <- gmat(g1, g2, p)

    for(i in 1:N){

      yi <- as.vector(y[i,])
      wi <- which(is.na(yi) == FALSE)
      pl <- length(wi)

      Si <- vmat(S[i,], p)

      yi <- yi[wi]
      Si <- pmat(Si, wi)
      Gi <- pmat(G, wi)

      Wi <- ginv(Gi + Si)

      A1 <- A1 + ivec(yi %*% Wi, wi, p)
      A2 <- A2 + imat(Wi, wi, p)

    }

    mu <- A1 %*% ginv(A2)
    g1 <- optimize(LL1, lower = 0, upper = 5)$minimum
    g2 <- 0.5*g1

    V.mu <- ginv(A2)

    Qc <- c(mu, g1, g2)

    rb <- abs(Qc - Qc0)/abs(Qc0); rb[is.nan(rb)] <- 0
    if(max(rb) < 10^-4) break


    Qc0 <- Qc

  }



  SE <- sqrt(diag(V.mu))

  R1 <- as.vector(mu)
  R2 <- as.vector(SE)
  R3 <- as.vector(mu - qnorm(.975)*SE)
  R4 <- as.vector(mu + qnorm(.975)*SE)

  R5 <- cbind(R1,R2,R3,R4); colnames(R5) <- c("Coef.", "SE", "95%CL", "95%CU")

  R6 <- sqrt(g1)
  R7 <- g2/g1

  R8 <- -LL1(g1)

  R9 <- list("Coefficients" = R5, "Between-studies_SD" = R6,
             "Between-studies_COR" = R7, "Loglikelihood" = R8)

  return(R9)

}