test_that("elicitDispersion and elicitGLM and fitGLM for unknown dispersion", { skip_on_cran() ## unknown dispersion: ## generate fake elicitations based on sample mean and covariance of data(longley) ## for each link function # true location and scale will be based on longley data (all positive values) dat <- longley # possible link functions links <- c("logit", "log", "cloglog", "probit", "cauchit", "log.sqrt", "identity", "log.inverse", "log.inverse.square") # test objects delta.ls <- S.ls <- delta0.ls <- S0.ls <- vector("list", length(links)) # location and scale names(delta.ls) <- names(S.ls) <- links r.vec <- s.vec <- numeric(length(links)) # dispersion names(r.vec) <- names(s.vec) <- links cond.probs <- matrix(NA, nrow = ncol(dat) - 1, # tail probabilities for conditioning values eta ncol = length(links), dimnames = list( level = 0:(ncol(dat) - 2), link = links )) # arbritrary true values specified for shape and rate dispersion prior parameters s.true <- 15 r.true <- 75 # true location and scale given elicited parameters r and s for dispersion prior # based on first and second moments of generalised multivariate t # on scale of linear predictor # scale longley dat by max value dat.true <- dat/max(dat) m.true <- colMeans(dat.true) V.true <- cov(dat.true)*(s.true - 2)/r.true for (j in links) { # link function link.true <- eglm:::linkFunction(j) # choose an observation model with an appropriate mean domain if (j%in%c("logit", "probit", "cauchit", "cloglog")) { obs.true <- "simplex" } else if (j%in%c("log.inverse")) { obs.true <- "inverseGaussian" } else if (j%in%c("log", "log.sqrt", "log.inverse.square")) { obs.true <- "gamma" } else if (j == "identity") { obs.true <- "normal" } ## specify elicitation scenarios matrix and initiate eglm object X <- matrix(1:ncol(dat), nrow = ncol(dat)) # scenarios Z <- eglm(scenarios = X, link = j, target = "target", expertID = "Expert", facilitator = "facilitator", known.dispersion = FALSE, obs.name = obs.true) ## true values for dispersion prediction intervals muhat <- mean(m.true) # linear predictor scale w <- 100 intervals <- eglm:::qgt(c((1 - Z$target$CI.prob)/2, 0.05), mu = muhat, sigma = sqrt(eglm:::var.fun(muhat, obs.name = obs.true)/w), r = r.true, s = s.true ) ######################## # dispersion elicitation Z <- elicitDispersion(Z, mu.hat = muhat, lower.bound1 = intervals[1], lower.bound2 = intervals[2], w = w) r.vec[j] <- Z$obs.model$r s.vec[j] <- Z$obs.model$s ####################### # marginal elicitations for (i in 1:ncol(dat)) { # interval on linear predictor scale then inverse link transformed interval <- link.true$linkinv(eglm:::qgt(c((1 - Z$target$CI.prob)/2, (1 + Z$target$CI.prob)/2), mu = m.true[i], sigma = sqrt(V.true[i,i]), r = Z$obs.model$r, s = Z$obs.model$s)) # inverted link functions must swap interval bounds if (j%in%c("log.inverse", "log.inverse.square")) interval <- rev(interval) Z <- elicitGLM(Z, scenario = i, CI.bounds = interval) } ################################# # conditional median elicitations for (l in 1:(ncol(V.true) - 1)) { Z <- elicitGLM(Z, scenario = l + 1, level = l) if (l == 1) { m.l.L <- Z$linpred$m[l] V.l.L <- Z$linpred$V[l, l] r.L <- Z$obs.model$r s.L <- Z$obs.model$s } else { inv.Lm1 <- chol2inv(chol(Z$linpred$V[1:(l - 1), 1:(l - 1)])) m.l.L <- Z$linpred$m[l] + + Z$linpred$V[l, 1:(l - 1)]%*%inv.Lm1%*%(Z$linpred$eta.cond[1:(l - 1)] - Z$linpred$m[1:(l - 1)]) V.l.L <- Z$linpred$V[l, l] - Z$linpred$V[l, 1:(l - 1)]%*%inv.Lm1%*%Z$linpred$V[1:(l - 1), l] z.l.L <- t(Z$linpred$eta.cond[1:(l - 1)] - Z$linpred$m[1:(l - 1)])%*%inv.Lm1%*%(Z$linpred$eta.cond[1:(l - 1)] - Z$linpred$m[1:(l - 1)]) r.L <- Z$obs.model$r + z.l.L s.L <- Z$obs.model$s + l - 1 } ############################################################################ # test: tail probabilities for conditioning values on linear predictor scale cond.probs[l, j] <- eglm:::pgt(Z$linpred$eta.cond[l], m.l.L, sqrt(V.l.L), r.L, s.L) if (Z$linpred$sgn.cond[l]) { true.tail.prob <- (1 - Z$target$CI.prob)/2 } else { true.tail.prob <- (1 + Z$target$CI.prob)/2 } expect_equal(cond.probs[!!l, !!j], true.tail.prob) inv <- chol2inv(chol(V.true[1:l, 1:l])) for (k in (l + 1):ncol(V.true)) { m.condtrue <- m.true[k] + V.true[k, 1:l]%*%inv%*%(Z$linpred$eta.cond[1:l] - m.true[1:l]) Z <- elicitGLM(Z, scenario = k, level = l, median = link.true$linkinv(m.condtrue)) } } ## elicited parametrisation of saturated model # truncated to level 0 (marginal elicitation only) eglm.t0 <- fitGLM(Z, X = diag(nrow(V.true)), trunc.level = 0) delta0.ls[[j]] <- eglm.t0$delta S0.ls[[j]] <- eglm.t0$S # full dependency eglm.t6 <- fitGLM(Z, X = diag(nrow(V.true))) delta.ls[[j]] <- eglm.t6$delta S.ls[[j]] <- eglm.t6$S ################################################################## # tests ## dispersion expect_equal(s.vec[!!j], s.true, tolerance = 1e-4, ignore_attr = TRUE) expect_equal(r.vec[!!j], r.true, tolerance = 1e-4, ignore_attr = TRUE) ## truncated eglm to level 0 (an Ieglm) matches true marginal scale expect_equal(delta0.ls[[!!j]], m.true, ignore_attr = TRUE) expect_equal(S0.ls[[!!j]], diag(diag(V.true)), tolerance = 1e-4, ignore_attr = TRUE) ## untruncated eglm matches true scale expect_equal(delta.ls[[!!j]], m.true, ignore_attr = TRUE) expect_equal(S.ls[[!!j]], V.true, tolerance = 1e-4, ignore_attr = TRUE) } # alternative variance function test: elicited gamma, fit inverse Gaussian fit.eglm.t6 <- fitGLM(Z, X = diag(nrow(V.true)), v.fun = function(mu) mu^3) # refitted rate r.refit <- r.true*eglm:::var.fun(muhat, obs.name = obs.true)/muhat^3 expect_equal(r.refit, fit.eglm.t6$r, tolerance = 1e-4, ignore_attr = TRUE) # refitted Sigma == V given that X = I expect_equal(Z$linpred$V*r.true/r.refit, fit.eglm.t6$Sigma, tolerance = 1e-4, ignore_attr = TRUE) # refitted delta should be unchanged expect_equal(eglm.t6$delta, fit.eglm.t6$delta, tolerance = 1e-4, ignore_attr = TRUE) })