search

hrue / r-inla

This is the public repository for the r-inla project
GNU General Public License v2.0
87 stars 24 forks source link

Hyperparameter sample name problem for experimental inla.mode #32

Closed finnlindgren closed 3 years ago

finnlindgren commented 3 years ago

With the experimental inla.mode="experimental" setting, the hyperparameter names from inla.posterior.sample are incorrect, with -- unknown added to the internal name, instead of just the plain non-internal name.

inla.posterior.sample, where the "rfake" object that is created doesn't contain internal.summary.hyperpar and summary.hyperpar, but those are the places the function INLA::inla.transform.names() looks for the parameter names to use. The result of those missing objects is that " -- unknown" gets added to the names, breaking the generate and predict() calls for inlabru calls that reference hyperparameters (and presumably also for any inla.posterior.sample.eval that does that). The problematic line is

rfake <- list(mlik = result$mlik, misc = list(from.theta = result$misc$from.theta, 
      to.theta = result$misc$to.theta, configs = result$misc$configs))

that should probably pass on the hyperpar summaries as well; not just mlik and misc. A possible fix:

rfake <- list(mlik = result$mlik, misc = list(from.theta = result$misc$from.theta, 
      to.theta = result$misc$to.theta, configs = result$misc$configs),
    summary.hyperpar = result$summary.hyperpar,
    internal.summary.hyperpar = result$internal.summary.hyperpar
    )

Reproducible example:

library(INLA)
#> Loading required package: Matrix
#> Loading required package: foreach
#> Loading required package: parallel
#> Loading required package: sp
#> This is INLA_21.07.10-1 built 2021-07-10 12:04:41 UTC.
#>  - See www.r-inla.org/contact-us for how to get help.
#>  - To enable PARDISO sparse library; see inla.pardiso()

# Fit a model with random effect z
# and repeated measurements (for identifiability) with
# additive Gaussian noise

df <- data.frame(z = rep(c(1, 10), each = 5))
df <- within(df,
                   y <- 5 +
                     rnorm(10, mean = 0, sd = 1)[z] +
                     rnorm(10, mean = 0, sd = 0.1))

# Estimation is OK for all three modes:
fit <- list()
for (inla.mode in c("classic", "twostage", "experimental")) {
  fit[[inla.mode]] <-
    inla(y ~ f(z, model = "iid"),
         family = "gaussian",
         data = df,
         control.compute = list(config = TRUE),
         inla.mode = inla.mode)
}

# inla.posterior.sample names are wrong for experimental, non-internal scale:
cat("inla.posterior.sample:")
#> inla.posterior.sample:
for (intern in c(FALSE, TRUE)) {
  for (inla.mode in c("classic", "twostage", "experimental")) {
    samples <-
      inla.posterior.sample(n = 1, result = fit[[inla.mode]], intern = intern)

    cat("Names for inla.mode='", inla.mode, "', intern = ", intern, ":\n  ",
        paste0(names(samples[[1]]$hyperpar), collapse = "\n  "),
        "\n", sep = "")
  }
}
#> Names for inla.mode='classic', intern = FALSE:
#>   Precision for the Gaussian observations
#>   Precision for z
#> Names for inla.mode='twostage', intern = FALSE:
#>   Precision for the Gaussian observations
#>   Precision for z
#> Names for inla.mode='experimental', intern = FALSE:
#>   Log precision for the Gaussian observations -- unknown
#>   Log precision for z -- unknown
#> Names for inla.mode='classic', intern = TRUE:
#>   Log precision for the Gaussian observations
#>   Log precision for z
#> Names for inla.mode='twostage', intern = TRUE:
#>   Log precision for the Gaussian observations
#>   Log precision for z
#> Names for inla.mode='experimental', intern = TRUE:
#>   Log precision for the Gaussian observations
#>   Log precision for z

# inla.hyperpar.sample names are correct:
cat("inla.hyperpar.sample:")
#> inla.hyperpar.sample:
for (intern in c(FALSE, TRUE)) {
  for (inla.mode in c("classic", "twostage", "experimental")) {
    samples <-
      inla.hyperpar.sample(n = 1, result = fit[[inla.mode]], intern = intern)

    cat("Names for inla.mode='", inla.mode, "', intern = ", intern, ":\n  ",
        paste0(colnames(samples), collapse = "\n  "),
        "\n", sep = "")
  }
}
#> Names for inla.mode='classic', intern = FALSE:
#>   Precision for the Gaussian observations
#>   Precision for z
#> Names for inla.mode='twostage', intern = FALSE:
#>   Precision for the Gaussian observations
#>   Precision for z
#> Names for inla.mode='experimental', intern = FALSE:
#>   Precision for the Gaussian observations
#>   Precision for z
#> Names for inla.mode='classic', intern = TRUE:
#>   Log precision for the Gaussian observations
#>   Log precision for z
#> Names for inla.mode='twostage', intern = TRUE:
#>   Log precision for the Gaussian observations
#>   Log precision for z
#> Names for inla.mode='experimental', intern = TRUE:
#>   Log precision for the Gaussian observations
#>   Log precision for z

Created on 2021-07-14 by the reprex package (v2.0.0)

hrue commented 3 years ago

Thx

— Håvard Rue Professor of Statistics Statistics Program, CEMSE Division King Abdullah University of Science and Technology Thuwal 23955-6900, Saudi Arabia

@. Office: +966 (0)12 808 0640 Mobile: +966 (0)54 470 0421 Research group: bayescomp.kaust.edu.sa R-INLA project: www.r-inla.org Zoom: kaust.zoom.us/my/haavard.rue On 14 Jul 2021, 17:14 +0200, Finn Lindgren @.>, wrote:

With the experimental inla.mode="experimental" setting, the hyperparameter names from inla.posterior.sample are incorrect, with -- unknown added to the internal name, instead of just the plain non-internal name. inla.posterior.sample, where the "rfake" object that is created doesn't contain internal.summary.hyperpar and summary.hyperpar, but those are the places the function INLA::inla.transform.names() looks for the parameter names to use. The result of those missing objects is that " -- unknown" gets added to the names, breaking the generate and predict() calls for inlabru calls that reference hyperparameters (and presumably also for any inla.posterior.sample.eval that does that). The problematic line is rfake <- list(mlik = result$mlik, misc = list(from.theta = result$misc$from.theta, to.theta = result$misc$to.theta, configs = result$misc$configs)) that should probably pass on the hyperpar summaries as well; not just mlik and misc. A possible fix: rfake <- list(mlik = result$mlik, misc = list(from.theta = result$misc$from.theta, to.theta = result$misc$to.theta, configs = result$misc$configs), summary.hyperpar = result$summary.hyperpar, internal.summary.hyperpar = result$internal.summary.hyperpar ) Reproducible example: library(INLA)

> Loading required package: Matrix

> Loading required package: foreach

> Loading required package: parallel

> Loading required package: sp

> This is INLA_21.07.10-1 built 2021-07-10 12:04:41 UTC.

> - See www.r-inla.org/contact-us for how to get help.

> - To enable PARDISO sparse library; see inla.pardiso()

Fit a model with random effect z

and repeated measurements (for identifiability) with

additive Gaussian noise

df <- data.frame(z = rep(c(1, 10), each = 5)) df <- within(df, y <- 5 + rnorm(10, mean = 0, sd = 1)[z] + rnorm(10, mean = 0, sd = 0.1))

Estimation is OK for all three modes:

fit <- list() for (inla.mode in c("classic", "twostage", "experimental")) { fit[[inla.mode]] <- inla(y ~ f(z, model = "iid"), family = "gaussian", data = df, control.compute = list(config = TRUE), inla.mode = inla.mode) }

inla.posterior.sample names are wrong for experimental, non-internal scale:

cat("inla.posterior.sample:")

> inla.posterior.sample:

for (intern in c(FALSE, TRUE)) { for (inla.mode in c("classic", "twostage", "experimental")) { samples <- inla.posterior.sample(n = 1, result = fit[[inla.mode]], intern = intern)

cat("Names for inla.mode='", inla.mode, "', intern = ", intern, ":\n ", paste0(names(samples[[1]]$hyperpar), collapse = "\n "), "\n", sep = "") } }

> Names for inla.mode='classic', intern = FALSE:

> Precision for the Gaussian observations

> Precision for z

> Names for inla.mode='twostage', intern = FALSE:

> Precision for the Gaussian observations

> Precision for z

> Names for inla.mode='experimental', intern = FALSE:

> Log precision for the Gaussian observations -- unknown

> Log precision for z -- unknown

> Names for inla.mode='classic', intern = TRUE:

> Log precision for the Gaussian observations

> Log precision for z

> Names for inla.mode='twostage', intern = TRUE:

> Log precision for the Gaussian observations

> Log precision for z

> Names for inla.mode='experimental', intern = TRUE:

> Log precision for the Gaussian observations

> Log precision for z

inla.hyperpar.sample names are correct:

cat("inla.hyperpar.sample:")

> inla.hyperpar.sample:

for (intern in c(FALSE, TRUE)) { for (inla.mode in c("classic", "twostage", "experimental")) { samples <- inla.hyperpar.sample(n = 1, result = fit[[inla.mode]], intern = intern)

cat("Names for inla.mode='", inla.mode, "', intern = ", intern, ":\n ", paste0(colnames(samples), collapse = "\n "), "\n", sep = "") } }

> Names for inla.mode='classic', intern = FALSE:

> Precision for the Gaussian observations

> Precision for z

> Names for inla.mode='twostage', intern = FALSE:

> Precision for the Gaussian observations

> Precision for z

> Names for inla.mode='experimental', intern = FALSE:

> Precision for the Gaussian observations

> Precision for z

> Names for inla.mode='classic', intern = TRUE:

> Log precision for the Gaussian observations

> Log precision for z

> Names for inla.mode='twostage', intern = TRUE:

> Log precision for the Gaussian observations

> Log precision for z

> Names for inla.mode='experimental', intern = TRUE:

> Log precision for the Gaussian observations

> Log precision for z

Created on 2021-07-14 by the reprex package (v2.0.0) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.

hrue commented 3 years ago

applied you fix. seems correct now

hrue commented 3 years ago

close