R/PBS0.LR.r

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

#' Parametric bootstrap for the likelihood ratio statistic
#'
#' Parametric bootstrap for the likelihood ratio statistic
#' to implement the improved tests developed by Noma et al. (2017).
#'
#' 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 A vector of (unconditional) maximum likelihood estimate of the random effects meta-analysis, computed by ML.
#' @param mu0 The value of the first component of the grand mean vector (corresponding to null hypothesis).
#' @param B Number of resampling.
#' @return Resampled B samples of the efficient score statistics.
#' @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
#'
#' # C <- 0.95
#' # alpha <- 1 - C
#'
#' # 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)
#'
#' # mu13 <- ml1$Coefficients[1, 1]
#' # ci3 <- ml1$Coefficients[1, 3:4]
#'
#' # R <- qchisq(C, df = 1)
#'
#' # beta1 <- log(beta1e)
#' # mlike0 <- .5*ml1$Loglikelihood
#'
#' # bcb1 <- PBS0.LR(y, S, ml0 = a4, mu0 = beta1, B = 400)
#'
#' # mlike1 <- .5*RML(y, S, ml0 = a4, mu0 = beta1)
#' # LR0 <- -2*(mlike1 - mlike0)
#' # LR1 <- LR0/mean(bcb1)
#'
#' # LR0  # ordinary likelihood ratio statistic
#' # LR1  # Bartlett corrected likelihood ratio statistic
#' @export
PBS0.LR <- function(y, S, ml0, mu0, B = 400){

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

  rml1 <- RMLP(y, S, ml0, mu0)
  G1 <- gmat(rml1$RML[p + 1], rml1$RML[p + 2], p)

  mu1 <- rml1$RML[1:p]

  y.pb <- matrix(numeric(N*p), N)

  LR0.b <- numeric(B)

  for(b in 1:B){

    for(i in 1:N){

      yi <- as.vector(y[i, ])
      wi <- which(is.na(yi) == FALSE)
      mSi <- vmat(S[i, ], p) + G1

      mui <- mu1[wi]
      mSi <- pmat(mSi, wi)

      y.pb[i,] <- ivec2(mvrnorm(1, mui, mSi), wi, p)

    }

    mlike1.b <- .5*RML(y.pb, S, ml0 = ml0, mu0 = mu0, maxitr = 50)
    mlike0.b <- .5*ML(y.pb, S, maxitr = 50)$Loglikelihood

    LR0.b[b] <- -2*(mlike1.b - mlike0.b)

  }

  return(LR0.b)

}