R/LRT.R

In aejensen/TruncComp: Two-sample comparison of truncated continuous outcome by parametric or semi-parametric maximum likelihood

getLRstartingValues <- function(d) {
  logOdds <- function(p) log(p / (1-p))

  mu0Init <- base::mean(d[d$R == 0 & d$A == 1, "Y"])
  muDeltaInit <- base::mean(d[d$R == 1 & d$A == 1, "Y"])  - mu0Init
  sigmaInit <- stats::sd(d[d$A == 1, "Y"])
  alphaInit <- logOdds(base::mean(d[d$R == 0, "A"]))
  alphaDeltaInit <- logOdds(base::mean(d[d$R == 1, "A"]))- alphaInit

  list(mu = mu0Init, muDelta = muDeltaInit,
       alpha = alphaInit, alphaDelta = alphaDeltaInit,
       sigma = sigmaInit)
}


LRT <- function(data, init, conf.level = 0.95) {
  a <- data$A
  y <- data$Y
  z <- data$R

  #A = 1(y_i is observed)
  expit <- function(x) exp(x)/(exp(x) + 1)

  dQoL <- function(A, Y, pA, eta, sigma) {
    stats::dnorm(Y, eta, sigma)^A * pA^A * (1 - pA)^(1 - A)
  }

  logLikH0 <- function(alpha, mu, logSigma) {
    #The log-likelihood under H_0: no treatment effect
    -sum(log(dQoL(a, y, expit(alpha), mu, exp(logSigma))))
  }

  logLikHA <- function(alpha, alphaDelta, mu, muDelta, logSigma) {
    #The log-likelihood under H_A: effects of treatment
    piA <- expit(alpha + alphaDelta * z)
    eta <- mu + muDelta * z
    sigma <- exp(logSigma)
    -sum(log(dQoL(a, y, piA, eta, sigma)))
  }

  #First, fit the model under H_0
  initH0 <- list(alpha = init$alpha, mu = init$mu, logSigma = log(init$sigma))
  mH0 <-  tryCatch(bbmle::mle2(logLikH0, start = initH0),
                  error = function(e) { NULL })
  if(is.null(mH0)) {
    return(NULL) #Return immediately on error
  }

  #Then fit the model under H_A
  initHA <- list(alpha = init$alpha, alphaDelta = init$alphaDelta, mu = init$mu,
                 muDelta = init$muDelta, logSigma = log(init$sigma))
  mHA <- tryCatch(bbmle::mle2(logLikHA, start = initHA),
                  error = function(e) { NULL })
  if(is.null(mHA)) {
    return(NULL)  #Return immediately on error
  }

  #We have arrived here without errors.
  #Get difference contrasts
  coefHA <- bbmle::coef(mHA)
  muDelta <- as.numeric(coefHA["muDelta"])
  alphaDelta <- as.numeric(exp(coefHA["alphaDelta"]))

  #Get confidence intervals for differences
  confintHA <- bbmle::confint(mHA, level = conf.level, method = "quad", quietly = TRUE)
  muDeltaCI <- as.numeric(confintHA["muDelta", ])
  alphaDeltaCI <- as.numeric(exp(confintHA["alphaDelta",]))

  #Get Delta
  delta <- mean(data$Y[data$R == 1]) - mean(data$Y[data$R == 0])
  deltaCI <- c(NA, NA)


  #Joint likelihood ratio test
  LRT <- bbmle::anova(mH0, mHA)
  W <- LRT[2,"Chisq"]
  p <- LRT[2, "Pr(>Chisq)"]

  out <- list(muDelta = muDelta,
              muDeltaCI = muDeltaCI,
              alphaDelta = alphaDelta,
              alphaDeltaCI = alphaDeltaCI,
              Delta = delta,
              DeltaCI = deltaCI,
              W = W,
              p = p,
              method = "Parametric Likelihood Ratio Test",
              conf.level = conf.level,
              success = TRUE,
              error = "",
              init = init)
  class(out) <- append(class(out), "TruncComp")
  out
}