Saving warning messages
Fernanda Vallejos
Hello everyone! I'm conducting a simulation study adjusting BYM models to different simulated datasets. In some scenarios, the model fits, but warning messages like these are generated:
I would like to save these messages to identify in which simulation they occurred, but I'm having trouble doing so. Please help me!
besag.iid.fun <- function(result) {
ady <- inla.read.graph(filename = "map.adj")
train_df <- dplyr::bind_cols(
data.frame(result$codigo_comuna),
x1masE = result$x1masE,
x2masE = result$x2masE,
E = result$E,
y = result$y,
real.RR = result$real.RR,
re_u = result$re_u,
re_v = result$re_v
)
# Lista para almacenar warnings
warnings_list <- list()
# Función para manejar advertencias
warning_handler <- function(w) {
warnings_list[[length(warnings_list) + 1]] <<- conditionMessage(w)
invokeRestart("muffleWarning")
}
# Ejecutar inla con tryCatch para capturar advertencias
resultadobym <- tryCatch({
withCallingHandlers(
inla(
formula = y ~ x1masE + x2masE +
f(re_u, model = "besag", graph = ady, scale.model = TRUE, constr = TRUE) +
f(re_v, model = "iid"),
family = "poisson",
data = train_df,
E = E,
control.compute = list(waic = TRUE, cpo = TRUE, dic = TRUE, return.marginals.predictor = TRUE),
control.predictor = list(compute = TRUE)
),
warning = warning_handler
)
}, error = function(e) e)
# Verificar si hubo un error
if (inherits(resultadobym, "error")) {
return(list(
error = "Se produjo un error al ejecutar inla."
))
}
# Extracción de resultados
efectos.fijos <- resultadobym$summary.fixed
waic <- resultadobym$waic$waic
dic <- resultadobym$dic$dic
cpo <- resultadobym$cpo$cpo
precision_u <- resultadobym$summary.hyperpar$mean[1]
precision_v <- resultadobym$summary.hyperpar$mean[2]
RR.ajustada.post <- resultadobym$summary.fitted.values[,"mean"]
LL.ajustada.post <- resultadobym$summary.fitted.values[, "0.025quant"]
UL.ajustada.post <- resultadobym$summary.fitted.values[, "0.975quant"]
sd_precision_u <- resultadobym$summary.hyperpar[1,2]
q1.precision_u <- resultadobym$summary.hyperpar[1,3]
q0.5precision_u <- resultadobym$summary.hyperpar[1,4]
q3precision_u <- resultadobym$summary.hyperpar[1,5]
sd_precision_v <- resultadobym$summary.hyperpar[2,2]
q1.precision_v <- resultadobym$summary.hyperpar[2,3]
q0.5precision_v <- resultadobym$summary.hyperpar[2,4]
q3precision_v <- resultadobym$summary.hyperpar[2,5]
x1masE <- result$x1masE
x2masE <- result$x2masE
# Determinar el mensaje de warning
warning_message <- if (length(warnings_list) > 0) {
paste(warnings_list, collapse = "; ")
} else {
"OK"
}
return(list(
efectos.fijos = efectos.fijos,
waic = waic,
dic = dic,
cpo = cpo,
x1masE = x1masE,
x2masE = x2masE,
precision_u = precision_u,
precision_v = precision_v,
real.RR = result$real.RR,
RR.ajustada.post = RR.ajustada.post,
UL.ajustada.post = UL.ajustada.post,
LL.ajustada.post = LL.ajustada.post,
sd_precision_u = sd_precision_u,
q1.precision_u = q1.precision_u,
q0.5precision_u = q0.5precision_u,
q3precision_u = q3precision_u,
sd_precision_v = sd_precision_v,
q1.precision_v = q1.precision_v,
q0.5precision_v = q0.5precision_v,
q3precision_v = q3precision_v,
warnings = warning_message # Agregar el mensaje de warning
))
}
Håvard Rue
x1masE + x2masE +
are somewhat scaled and centered.
the 'warning' that is show, says that the correction for the mean computed, is
to large to be trusted, and this usually happens when the model is to complex
for the data, or the other way around
On Sun, 2024-06-23 at 07:18 -0700, Fernanda Vallejos wrote:
> x1masE + x2masE +
--
Håvard Rue
hr...@r-inla.org
Fernanda Vallejos
Helpdesk (Haavard Rue)
try to rescale covariates and put informative priors for the fixed effects and
also intercept. you want to reach a model without the warning
Helpdesk (Haavard Rue)
formula = y ~ x1masE + x2masE +
f(re_u, model = "besag", graph = ady, scale.model = TRUE,
hyper=list(prec=list(prior="pc.prec", param=c(0.5, 0.01)))) +
f(re_v, model = "iid",
hyper=list(prec=list(prior="pc.prec", param=c(0.5, 0.01)))),
family = "poisson",
control.fixed = list(prec.intercept=1, prec=1)
Fernanda Vallejos
Helpdesk (Haavard Rue)
y <- rpois(n, exp(2))
idx <- 1:n
for (emerg in c(0.001, 25)) {
cat("emergency = ", emerg, "\n")
r <- inla(y ~ 1 + f(idx, vb.correct = T, initial = log(1), fixed = T),
data = data.frame(y, idx),
control.inla = list(control.vb = list(emergency = emerg),
int.strategy = "eb"),
family = "poisson")
if (length(grep("max_correction = ", r$logfile)) > 0) {
cat("VB ABORTED\n")
} else {
cat("VB OK\n")