test_that("elicitGLM for known dispersion matches densities", {

  skip_on_cran()

  # scenarios matrix: two scenarios
  X <- matrix(c(1, 1, 0, 1), nrow = 2)
  rownames(X) <- c("scenario1", "scenario2")
  colnames(X) <- c("covariate1", "covariate2")

  # logit link
  # central credible intervals with tertiles
  Z <- eglm(scenarios = X, link = "logit")

  # 1st scenario
  # no elicited fractiles
  Z <- elicitGLM(Z, scenario = 1)
  # elicited median
  Z <- elicitGLM(Z, scenario = 1, CI.bounds = c(0.2, 0.4))

  # second scenario
  # elicit upper and lower tertiles
  Z <- elicitGLM(Z, scenario = 2, CI.bounds = c(0.1, 0.3))

  # conditional median at level 1, scenario 2
  # initially assumed conditionally independent
  # grey area indicates coherent range for median
  # conditional on responses at level 0
  Z <- elicitGLM(Z, scenario = 2, level = 1)

  # specify conditional median
  Z <- elicitGLM(Z, scenario = 2, level = 1, median = 0.25)

  # the current elicitation can be removed by setting equal to NA
  Z <- elicitGLM(Z, scenario = 2, level = 1, median = NA)

  # Change conditional median to 0.2 and calculate
  # cumulative probability that mu is less than 0.5
  Z <- elicitGLM(Z, scenario = 2, level = 1, median = 0.2, evalCDF = 0.5)

  #################
  # check densities
  tZ <- fitGLM(Z, X = diag(2))

  # coordinates
  xx <- seq(from = 0, to = 1, length.out = 101)

  # scenario 1, level 0
  # graphical density check
  elicitGLM(Z, scenario = 1)
  curve(eglm:::dLogitNorm(x, tZ$delta[1], sqrt(tZ$Sigma[1, 1])), add = TRUE, col = 'red')
  # probs
  capture <- capture.output(result <- elicitGLM(Z, scenario = 1))
  out <- strsplit(capture[13], " ")[[1]]
  nums <- as.numeric(out)
  fracs <- nums[!is.na(nums)]
  expect_equal(c((1 - Z$target$CI.prob)/2, 1/2, (1 + Z$target$CI.prob)/2),
               pnorm(log(fracs/(1 - fracs)), tZ$delta[1], sqrt(tZ$Sigma[1, 1])),
               tolerance = 1e-4)


  # scenario 2, level 0
  elicitGLM(Z, scenario = 2)
  curve(eglm:::dLogitNorm(x, tZ$delta[2], sqrt(tZ$Sigma[2, 2])), add = TRUE, col = 'red')
  # probs
  capture <- capture.output(result <- elicitGLM(Z, scenario = 2))
  out <- strsplit(capture[13], " ")[[1]]
  nums <- as.numeric(out)
  fracs <- nums[!is.na(nums)]
  expect_equal(c((1 - Z$target$CI.prob)/2, 1/2, (1 + Z$target$CI.prob)/2),
               pnorm(log(fracs/(1 - fracs)), tZ$delta[2], sqrt(tZ$Sigma[2, 2])),
               tolerance = 1e-4)

  # scenario 2, level = 1
  elicitGLM(Z, scenario = 2, level = 1)
  curve(eglm:::dLogitNorm(x,
    tZ$delta[2] + tZ$Sigma[1, 2]*(Z$linpred$eta.cond[1] - tZ$delta[1])/tZ$Sigma[1, 1],
     sqrt(tZ$Sigma[2, 2] - tZ$Sigma[1, 2]^2/tZ$Sigma[1, 1])),
     add = TRUE, col = 'red')
  # probs
  capture <- capture.output(result <- elicitGLM(Z, scenario = 2, level = 1))
  out <- strsplit(capture[13], " ")[[1]]
  nums <- as.numeric(out)
  fracs <- nums[!is.na(nums)]
  expect_equal(c((1 - Z$target$CI.prob)/2, 1/2, (1 + Z$target$CI.prob)/2),
               pnorm(log(fracs/(1 - fracs)),
                     tZ$delta[2] + tZ$Sigma[1, 2]*(Z$linpred$eta.cond[1] - tZ$delta[1])/tZ$Sigma[1, 1],
                     sqrt(tZ$Sigma[2, 2] - tZ$Sigma[1, 2]^2/tZ$Sigma[1, 1])),
               tolerance = 1e-4)
})