R/ES.r

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

#' Calculating the efficient score statistic
#'
#' Calculating the efficient score statistic
#'
#' 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 ml0 Initial value of the grand mean vector except for the first component.
#' @param mu0 The value of the first component of the grand mean vector.
#' @param maxitr The maximum iteration number of the Newton-Raphson algorithm.
#' @return The value of the efficient score statistic.
#' @references
#' 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
#'
#' # beta1e <- 0.80
#'
#' # ml1 <- ML(y, S)
#'
#' # a1 <- ml1$Coefficients[, 1]
#' # a2 <- (ml1$`Between-studies_SD`)^2
#' # a3 <- a2*(ml1$`Between-studies_COR`)
#' # a4 <- c(a1, a2, a3)
#'
#' # beta1 <- log(beta1e)
#'
#' # ES0 <- ES(y, S, ml0 = a4 , mu0 = beta1)
#' @export
ES <- function(y, S, ml0, mu0, maxitr = 200){

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

  muc <- ml0[2:p]		# initial values
  mu <- c(mu0, muc)
  g1 <- ml0[p+1]
  g2 <- ml0[p+2]

  Qc0 <- c(muc, 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(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)   # return "minus loglikelihood"

  }

  for(itr in 1:maxitr){

    A1 <- A2 <- A3 <- numeric(p - 1)
    A4 <- matrix(numeric((p - 1)*(p - 1)),(p - 1))

    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)

      Wi <- imat(Wi, wi, p)
      yi <- ivec(yi, wi, p)

      Wi21 <- Wi[2:p,1]
      Wi22 <- Wi[2:p,2:p]

      A1 <- A1 + Wi21 * yi[1]
      A2 <- A2 + Wi22 %*% yi[2:p]
      A3 <- A3 + Wi21 * mu0
      A4 <- A4 + Wi22

    }

    muc <- as.vector(A1 + A2 - A3) %*% ginv(A4)
    mu <- c(mu0, muc)

    g1 <- optimize(LL1, lower = 0, upper = 5)$minimum
    g2 <- 0.5*g1

    Qc <- c(muc,g1,g2)

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

    Qc0 <- Qc

  }

  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)

    Wi <- imat(Wi, wi, p)
    yi <- ivec(yi, wi, p)

    A1 <- A1 + as.numeric(Wi %*% (yi - mu))
    A2 <- A2 + Wi

  }

  R1 <- as.numeric(t(A1) %*% ginv(A2) %*% A1)

  return(R1)

}