Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Automated predictor size detection needed for intercept-only models (Was: Error fitting homogeneous Poisson process with compressed cp likelihood) #237

Open
ASeatonSpatial opened this issue May 6, 2024 · 8 comments
Open

Automated predictor size detection needed for intercept-only models (Was: Error fitting homogeneous Poisson process with compressed cp likelihood) #237

ASeatonSpatial opened this issue May 6, 2024 · 8 comments
Assignees
@finnlindgren
Labels
enhancement

Comments

@ASeatonSpatial
Copy link
Collaborator

I am fitting a homogeneous Poisson process model using the cp likelihood. To do this I want to define the integration scheme as with a single integration location with a weight equal to the area of the observation window. However, this seems to cause an error for the compressed cp likelihood.

When I use the option bru_compress_cp = FALSE, the model will fit with a warning message about an array of length 1, but the results look sensible. When I use bru_compress_cp = TRUE I get an error and the model will not fit. It looks like in the second case inlabru tries to compute a linearised model configuration despite the formula being defined as geometry ~ ..

I thought initially there might be a problem with having just a single integration location. But I tried adding more and I still hit the same warnings and error messages.

Here is a reprex:

library(inlabru)
#> Loading required package: fmesher
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(ggplot2)
set.seed(183)

# region where points are observable defined as a
# square polygon with area 9 and stored as an sf polygon object.
loc.d <- 3 * cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0))
loc <- st_polygon(list(loc.d))
Omega <- st_sf(st_sfc(loc))

# set true intensity
true_lambda <- 10

# simulate the number of points in the pattern
N <- rpois(1, true_lambda * 9)

# simulate random locations using st_sample()
pp <- st_sf(
  geometry = st_sample(Omega, N)
)

# define integration scheme
# since homogeneous just 1 location with weight equal
# to area of the observation window
ips <- st_sf(
  geometry = st_sample(Omega, 1)
)
ips$weight <- 9

# components - just the log intensity parameter
cmp <- ~ 0 + loglambda(1)

# fit without compressed cp likelihood
# runs but with a warning message
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = FALSE
  )
)

fit <- bru(lik,
  components = cmp
)
#> Warning in pred_eps - pred0: Recycling array of length 1 in array-vector arithmetic is deprecated.
#>   Use c() or as.vector() instead.

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.504, Running = 0.194, Post = 0.036, Total = 0.734 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -609.29
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -182.79
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -340.04
#> Effective number of parameters .................: 7.80
#> 
#> Marginal log-Likelihood:  116.14 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

# this should be around 10
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# fit with compressed cp likelihood
# gives an error
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = TRUE
  )
)

fit <- bru(lik,
  components = cmp
)
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor 'loglambda' plus eps.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor derivatives for 'loglambda'; treated as 0.0.
#> 
#>  *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1
#> Warning in iinla(model = info[["model"]], lhoods = info[["lhoods"]], options = info[["options"]]): iinla: Problem in inla: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : 
#>   2 rows in 'control.predictor=list(A=A)'-matrix, contain only zeros. This is not allowed.
#> The inla program failed and the maximum number of tries has been reached.
#> iinla: Problem in inla: 1: base::tryCatch(base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("C:\\ [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "C:\\Users\\ANDYWO~1\\AppData\\Local\\Temp\\ [...]
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    }), error = function(e) {
#>        NULL
#>        if (FALSE) {
#>            base::try(base::stop(e))
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    }, interrupt = function(e) {
#>        NULL
#>        if (FALSE) {
#>            e
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    })
#> 2: tryCatchList(expr, classes, parentenv, handlers)
#> 3: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
#>        names[nh], parentenv, handlers[[nh]])
#> 4: doTryCatch(return(expr), name, parentenv, handler)
#> 5: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
#> 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
#> 7: doTryCatch(return(expr), name, parentenv, handler)
#> 8: base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("C:\\ [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "C:\\Users\\ANDYWO~1\\AppData\\Local\\Temp\\ [...]
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    })
#> 9: base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("C:\\User [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE), file = "C:\\Users\\ANDYWO~1\\AppData\\Local\\Temp\\Rtmp [...]
#>        compress = FALSE)
#> 10: base::do.call(base::do.call, base::c(base::readRDS("C:\\Users\\ANDYWO~1\\A [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE)
#> 11: (function (what, args, quote = FALSE, envir = parent.frame()) 
#>    {
#>        if (!is.list(args)) 
#>            stop("second argument must be a list")
#>        if (quote) 
#>            args <- lapply(args, enquote)
#>        .Internal(do.call(what, args, envir))
#>    })(base::quote(function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    }), base::quote(list(input = "grave-cobra_reprex.R")), envir = base::quote [...]
#>        quote = base::quote(TRUE))
#> 12: (function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    })(input = base::quote("grave-cobra_reprex.R"))
#> 13: rmarkdown::render(input, quiet = TRUE, envir = globalenv(), encoding = "UTF-8")
#> 14: knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
#> 15: process_file(text, output)
#> 16: xfun:::handle_error(withCallingHandlers(if (tangle) process_tangle(group)  [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#>        function(e, loc) {
#>            setwd(wd)
#>            write_utf8(res, output %n% stdout())
#>            paste0("\nQuitting from lines ", loc)
#>        }, if (labels[i] != "") sprintf(" [%s]", labels[i]), get_loc)
#> 17: withCallingHandlers(expr, error = function(e) {
#>        loc = if (is.function(fun)) 
#>            trimws(fun(label))
#>        else ""
#>        if (loc != "") 
#>            loc = sprintf(" at lines %s", loc)
#>        message(one_string(handler(e, loc)))
#>    })
#> 18: withCallingHandlers(if (tangle) process_tangle(group) else process_group(g [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#> 19: process_gr
#> iinla: Giving up and returning last successfully obtained result for diagnostic purposes.
#> Error in track[[1]]: subscript out of bounds

# try an integration scheme with 9 points
#########################################

ips <- st_sf(
  geometry = st_sample(Omega, 9)
)
ips$weight <- 1

# non-compressed cp likelihood
# still works
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = FALSE
  )
)

fit <- bru(lik,
           components = cmp
)
#> Warning in pred_eps - pred0: Recycling array of length 1 in array-vector arithmetic is deprecated.
#>   Use c() or as.vector() instead.

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.187, Running = 0.187, Post = 0.0313, Total = 0.405 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -241.72
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: 0.997
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -230.43
#> Effective number of parameters .................: 11.28
#> 
#> Marginal log-Likelihood:  116.14 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# compressed cp likelihood
# still gives an error
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = TRUE
  )
)

fit <- bru(lik,
           components = cmp
)
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor 'loglambda' plus eps.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor derivatives for 'loglambda'; treated as 0.0.
#> 
#>  *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1
#> Warning in iinla(model = info[["model"]], lhoods = info[["lhoods"]], options = info[["options"]]): iinla: Problem in inla: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : 
#>   10 rows in 'control.predictor=list(A=A)'-matrix, contain only zeros. This is not allowed.
#> The inla program failed and the maximum number of tries has been reached.
#> iinla: Problem in inla: 1: base::tryCatch(base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("C:\\ [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "C:\\Users\\ANDYWO~1\\AppData\\Local\\Temp\\ [...]
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    }), error = function(e) {
#>        NULL
#>        if (FALSE) {
#>            base::try(base::stop(e))
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    }, interrupt = function(e) {
#>        NULL
#>        if (FALSE) {
#>            e
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    })
#> 2: tryCatchList(expr, classes, parentenv, handlers)
#> 3: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
#>        names[nh], parentenv, handlers[[nh]])
#> 4: doTryCatch(return(expr), name, parentenv, handler)
#> 5: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
#> 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
#> 7: doTryCatch(return(expr), name, parentenv, handler)
#> 8: base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("C:\\ [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "C:\\Users\\ANDYWO~1\\AppData\\Local\\Temp\\ [...]
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("C:\ [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    })
#> 9: base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("C:\\User [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE), file = "C:\\Users\\ANDYWO~1\\AppData\\Local\\Temp\\Rtmp [...]
#>        compress = FALSE)
#> 10: base::do.call(base::do.call, base::c(base::readRDS("C:\\Users\\ANDYWO~1\\A [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE)
#> 11: (function (what, args, quote = FALSE, envir = parent.frame()) 
#>    {
#>        if (!is.list(args)) 
#>            stop("second argument must be a list")
#>        if (quote) 
#>            args <- lapply(args, enquote)
#>        .Internal(do.call(what, args, envir))
#>    })(base::quote(function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    }), base::quote(list(input = "grave-cobra_reprex.R")), envir = base::quote [...]
#>        quote = base::quote(TRUE))
#> 12: (function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    })(input = base::quote("grave-cobra_reprex.R"))
#> 13: rmarkdown::render(input, quiet = TRUE, envir = globalenv(), encoding = "UTF-8")
#> 14: knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
#> 15: process_file(text, output)
#> 16: xfun:::handle_error(withCallingHandlers(if (tangle) process_tangle(group)  [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#>        function(e, loc) {
#>            setwd(wd)
#>            write_utf8(res, output %n% stdout())
#>            paste0("\nQuitting from lines ", loc)
#>        }, if (labels[i] != "") sprintf(" [%s]", labels[i]), get_loc)
#> 17: withCallingHandlers(expr, error = function(e) {
#>        loc = if (is.function(fun)) 
#>            trimws(fun(label))
#>        else ""
#>        if (loc != "") 
#>            loc = sprintf(" at lines %s", loc)
#>        message(one_string(handler(e, loc)))
#>    })
#> 18: withCallingHandlers(if (tangle) process_tangle(group) else process_group(g [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#> 19: process_g
#> iinla: Giving up and returning last successfully obtained result for diagnostic purposes.
#> Error in track[[1]]: subscript out of bounds

# session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.3.3 (2024-02-29 ucrt)
#>  os       Windows 10 x64 (build 19045)
#>  system   x86_64, mingw32
#>  ui       RTerm
#>  language (EN)
#>  collate  English_United Kingdom.utf8
#>  ctype    English_United Kingdom.utf8
#>  tz       Europe/London
#>  date     2024-05-06
#>  pandoc   3.1.11 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package      * version     date (UTC) lib source
#>  cachem         1.0.8       2023-05-01 [1] CRAN (R 4.3.3)
#>  class          7.3-22      2023-05-03 [2] CRAN (R 4.3.3)
#>  classInt       0.4-10      2023-09-05 [1] CRAN (R 4.3.3)
#>  cli            3.6.2       2023-12-11 [1] CRAN (R 4.3.3)
#>  colorspace     2.1-0       2023-01-23 [1] CRAN (R 4.3.3)
#>  DBI            1.2.2       2024-02-16 [1] CRAN (R 4.3.3)
#>  devtools       2.4.5       2022-10-11 [1] CRAN (R 4.3.3)
#>  digest         0.6.35      2024-03-11 [1] CRAN (R 4.3.3)
#>  dplyr          1.1.4       2023-11-17 [1] CRAN (R 4.3.3)
#>  e1071          1.7-14      2023-12-06 [1] CRAN (R 4.3.3)
#>  ellipsis       0.3.2       2021-04-29 [1] CRAN (R 4.3.3)
#>  evaluate       0.23        2023-11-01 [1] CRAN (R 4.3.3)
#>  fansi          1.0.6       2023-12-08 [1] CRAN (R 4.3.3)
#>  fastmap        1.1.1       2023-02-24 [1] CRAN (R 4.3.3)
#>  fmesher      * 0.1.5.9004  2024-05-06 [1] Github (inlabru-org/fmesher@c4b05ad)
#>  fs             1.6.4       2024-04-25 [1] CRAN (R 4.3.3)
#>  generics       0.1.3       2022-07-05 [1] CRAN (R 4.3.3)
#>  ggplot2      * 3.5.1       2024-04-23 [1] CRAN (R 4.3.3)
#>  glue           1.7.0       2024-01-09 [1] CRAN (R 4.3.3)
#>  gtable         0.3.5       2024-04-22 [1] CRAN (R 4.3.3)
#>  htmltools      0.5.8.1     2024-04-04 [1] CRAN (R 4.3.3)
#>  htmlwidgets    1.6.4       2023-12-06 [1] CRAN (R 4.3.3)
#>  httpuv         1.6.15      2024-03-26 [1] CRAN (R 4.3.3)
#>  INLA           24.05.01-1  2024-05-01 [1] local
#>  inlabru      * 2.10.1.9005 2024-05-06 [1] Github (inlabru-org/inlabru@8d1d012)
#>  KernSmooth     2.23-22     2023-07-10 [2] CRAN (R 4.3.3)
#>  knitr          1.46        2024-04-06 [1] CRAN (R 4.3.3)
#>  later          1.3.2       2023-12-06 [1] CRAN (R 4.3.3)
#>  lattice        0.22-5      2023-10-24 [2] CRAN (R 4.3.3)
#>  lifecycle      1.0.4       2023-11-07 [1] CRAN (R 4.3.3)
#>  magrittr       2.0.3       2022-03-30 [1] CRAN (R 4.3.3)
#>  Matrix         1.6-5       2024-01-11 [2] CRAN (R 4.3.3)
#>  MatrixModels   0.5-3       2023-11-06 [1] CRAN (R 4.3.3)
#>  memoise        2.0.1       2021-11-26 [1] CRAN (R 4.3.3)
#>  mime           0.12        2021-09-28 [1] CRAN (R 4.3.1)
#>  miniUI         0.1.1.1     2018-05-18 [1] CRAN (R 4.3.3)
#>  munsell        0.5.1       2024-04-01 [1] CRAN (R 4.3.3)
#>  pillar         1.9.0       2023-03-22 [1] CRAN (R 4.3.3)
#>  pkgbuild       1.4.4       2024-03-17 [1] CRAN (R 4.3.3)
#>  pkgconfig      2.0.3       2019-09-22 [1] CRAN (R 4.3.3)
#>  pkgload        1.3.4       2024-01-16 [1] CRAN (R 4.3.3)
#>  plyr           1.8.9       2023-10-02 [1] CRAN (R 4.3.3)
#>  profvis        0.3.8       2023-05-02 [1] CRAN (R 4.3.3)
#>  promises       1.3.0       2024-04-05 [1] CRAN (R 4.3.3)
#>  proxy          0.4-27      2022-06-09 [1] CRAN (R 4.3.3)
#>  purrr          1.0.2       2023-08-10 [1] CRAN (R 4.3.3)
#>  R.cache        0.16.0      2022-07-21 [1] CRAN (R 4.3.3)
#>  R.methodsS3    1.8.2       2022-06-13 [1] CRAN (R 4.3.3)
#>  R.oo           1.26.0      2024-01-24 [1] CRAN (R 4.3.3)
#>  R.utils        2.12.3      2023-11-18 [1] CRAN (R 4.3.3)
#>  R6             2.5.1       2021-08-19 [1] CRAN (R 4.3.3)
#>  Rcpp           1.0.12      2024-01-09 [1] CRAN (R 4.3.3)
#>  remotes        2.5.0       2024-03-17 [1] CRAN (R 4.3.3)
#>  reprex         2.1.0       2024-01-11 [1] CRAN (R 4.3.3)
#>  rlang          1.1.3       2024-01-10 [1] CRAN (R 4.3.3)
#>  rmarkdown      2.26        2024-03-05 [1] CRAN (R 4.3.3)
#>  rstudioapi     0.16.0      2024-03-24 [1] CRAN (R 4.3.3)
#>  scales         1.3.0       2023-11-28 [1] CRAN (R 4.3.3)
#>  sessioninfo    1.2.2       2021-12-06 [1] CRAN (R 4.3.3)
#>  sf           * 1.0-16      2024-03-24 [1] CRAN (R 4.3.3)
#>  shiny          1.8.1.1     2024-04-02 [1] CRAN (R 4.3.3)
#>  sp             2.1-4       2024-04-30 [1] CRAN (R 4.3.3)
#>  stringi        1.8.3       2023-12-11 [1] CRAN (R 4.3.2)
#>  stringr        1.5.1       2023-11-14 [1] CRAN (R 4.3.3)
#>  styler         1.10.3      2024-04-07 [1] CRAN (R 4.3.3)
#>  tibble         3.2.1       2023-03-20 [1] CRAN (R 4.3.3)
#>  tidyselect     1.2.1       2024-03-11 [1] CRAN (R 4.3.3)
#>  units          0.8-5       2023-11-28 [1] CRAN (R 4.3.3)
#>  urlchecker     1.0.1       2021-11-30 [1] CRAN (R 4.3.3)
#>  usethis        2.2.3       2024-02-19 [1] CRAN (R 4.3.3)
#>  utf8           1.2.4       2023-10-22 [1] CRAN (R 4.3.3)
#>  vctrs          0.6.5       2023-12-01 [1] CRAN (R 4.3.3)
#>  withr          3.0.0       2024-01-16 [1] CRAN (R 4.3.3)
#>  xfun           0.43        2024-03-25 [1] CRAN (R 4.3.3)
#>  xtable         1.8-4       2019-04-21 [1] CRAN (R 4.3.3)
#>  yaml           2.3.8       2023-12-11 [1] CRAN (R 4.3.2)
#> 
#>  [1] C:/Users/Andy Work/AppData/Local/R/win-library/4.3
#>  [2] C:/Program Files/R/R-4.3.3/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Created on 2024-05-06 with reprex v2.1.0
@finnlindgren
Copy link
Collaborator

This looks like it’s related to the problem one runs into with models with just a single intercept component (and no other components), where it can’t know the size of the predictor, and/or whether the “1” needs to be expanded into a vector or not. Try loglambda(rep(1, nrow(.data.))) instead of loglambda(1). I haven’t spent much time trying to work out if this special case can be handled more smoothly or not, since models with only and intercept are so rare in practice and the workaround is simple.

@finnlindgren
Copy link
Collaborator

I ran it with cmp <- ~ 0 + loglambda(rep(1, nrow(.data.))) and that works, with both bru_compress_cp = FALSE and bru_compress_cp = TRUE, whereas TRUE failed for loglambda(1), so the workaround works. However, I'm not getting the same results for the original reprex. Investigating.

@finnlindgren
Copy link
Collaborator

My reprex output:

library(inlabru)
#> Loading required package: fmesher
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.7.1, PROJ 9.2.1; sf_use_s2() is TRUE
library(ggplot2)
set.seed(183)

# region where points are observable defined as a
# square polygon with area 9 and stored as an sf polygon object.
loc.d <- 3 * cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0))
loc <- st_polygon(list(loc.d))
Omega <- st_sf(st_sfc(loc))

# set true intensity
true_lambda <- 10

# simulate the number of points in the pattern
N <- rpois(1, true_lambda * 9)

# simulate random locations using st_sample()
pp <- st_sf(
  geometry = st_sample(Omega, N)
)

# define integration scheme
# since homogeneous just 1 location with weight equal
# to area of the observation window
ips <- st_sf(
  geometry = st_sample(Omega, 1)
)
ips$weight <- 9

# components - just the log intensity parameter
cmp <- ~ 0 + loglambda(1)

# fit without compressed cp likelihood
# runs but with a warning message
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = FALSE
  )
)

fit <- bru(lik,
           components = cmp
)
#> Warning in pred_eps - pred0: Recycling array of length 1 in array-vector arithmetic is deprecated.
#>   Use c() or as.vector() instead.

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.431, Running = 0.176, Post = 0.0192, Total = 0.627 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -609.29
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -182.79
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -340.04
#> Effective number of parameters .................: 7.80
#> 
#> Marginal log-Likelihood:  116.14 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

# this should be around 10
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# fit with compressed cp likelihood
# gives an error
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = TRUE
  )
)

fit <- bru(lik,
           components = cmp
)
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor 'loglambda' plus eps.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor derivatives for 'loglambda'; treated as 0.0.
#> 
#>  *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1
#> Warning in iinla(model = info[["model"]], lhoods = info[["lhoods"]], options = info[["options"]]): iinla: Problem in inla: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : 
#>   2 rows in 'control.predictor=list(A=A)'-matrix, contain only zeros. This is not allowed.
#> The inla program failed and the maximum number of tries has been reached.
#> iinla: Problem in inla: 1: base::tryCatch(base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "/tmp/RtmpWxJiB4/callr-res-29ba6f0677ca", 
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    }), error = function(e) {
#>        NULL
#>        if (FALSE) {
#>            base::try(base::stop(e))
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    }, interrupt = function(e) {
#>        NULL
#>        if (FALSE) {
#>            e
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    })
#> 2: tryCatchList(expr, classes, parentenv, handlers)
#> 3: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
#>        names[nh], parentenv, handlers[[nh]])
#> 4: doTryCatch(return(expr), name, parentenv, handler)
#> 5: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
#> 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
#> 7: doTryCatch(return(expr), name, parentenv, handler)
#> 8: base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "/tmp/RtmpWxJiB4/callr-res-29ba6f0677ca", 
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    })
#> 9: base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp/Rtm [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE), file = "/tmp/RtmpWxJiB4/callr-res-29ba6f0677ca", 
#>        compress = FALSE)
#> 10: base::do.call(base::do.call, base::c(base::readRDS("/tmp/RtmpWxJiB4/callr- [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE)
#> 11: (function (what, args, quote = FALSE, envir = parent.frame()) 
#>    {
#>        if (!is.list(args)) 
#>            stop("second argument must be a list")
#>        if (quote) 
#>            args <- lapply(args, enquote)
#>        .Internal(do.call(what, args, envir))
#>    })(base::quote(function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    }), base::quote(list(input = "used-aidi_reprex.R")), envir = base::quote(< [...]
#>        quote = base::quote(TRUE))
#> 12: (function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    })(input = base::quote("used-aidi_reprex.R"))
#> 13: rmarkdown::render(input, quiet = TRUE, envir = globalenv(), encoding = "UTF-8")
#> 14: knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
#> 15: process_file(text, output)
#> 16: xfun:::handle_error(withCallingHandlers(if (tangle) process_tangle(group)  [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#>        function(e, loc) {
#>            setwd(wd)
#>            write_utf8(res, output %n% stdout())
#>            paste0("\nQuitting from lines ", loc)
#>        }, if (labels[i] != "") sprintf(" [%s]", labels[i]), get_loc)
#> 17: withCallingHandlers(expr, error = function(e) {
#>        loc = if (is.function(fun)) 
#>            trimws(fun(label))
#>        else ""
#>        if (loc != "") 
#>            loc = sprintf(" at lines %s", loc)
#>        message(one_string(handler(e, loc)))
#>    })
#> 18: withCallingHandlers(if (tangle) process_tangle(group) else process_group(g [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#> 19: process_group(group)
#> 20: call_block(x)
#> 2
#> iinla: Giving up and returning last successfully obtained result for diagnostic purposes.
#> Error in track[[1]]: subscript out of bounds


# try an integration scheme with 9 points
#########################################

ips <- st_sf(
  geometry = st_sample(Omega, 9)
)
ips$weight <- 1

# non-compressed cp likelihood
# still works
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = FALSE
  )
)

fit <- bru(lik,
           components = cmp
)
#> Warning in pred_eps - pred0: Recycling array of length 1 in array-vector arithmetic is deprecated.
#>   Use c() or as.vector() instead.

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(1), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.145, Running = 0.176, Post = 0.00956, Total = 0.331 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -241.72
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: 0.997
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -230.43
#> Effective number of parameters .................: 11.28
#> 
#> Marginal log-Likelihood:  116.14 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# compressed cp likelihood
# still gives an error
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = TRUE
  )
)

fit <- bru(lik,
           components = cmp
)
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor 'loglambda' plus eps.
#> Warning in bru_compute_linearisation.component(model[["effects"]][[label]], : Non-finite (-Inf/Inf/NaN) entries detected in predictor derivatives for 'loglambda'; treated as 0.0.
#> 
#>  *** inla.core.safe:  The inla program failed, but will rerun in case better initial values may help. try=1/1
#> Warning in iinla(model = info[["model"]], lhoods = info[["lhoods"]], options = info[["options"]]): iinla: Problem in inla: Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts,  : 
#>   10 rows in 'control.predictor=list(A=A)'-matrix, contain only zeros. This is not allowed.
#> The inla program failed and the maximum number of tries has been reached.
#> iinla: Problem in inla: 1: base::tryCatch(base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "/tmp/RtmpWxJiB4/callr-res-29ba6f0677ca", 
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    }), error = function(e) {
#>        NULL
#>        if (FALSE) {
#>            base::try(base::stop(e))
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    }, interrupt = function(e) {
#>        NULL
#>        if (FALSE) {
#>            e
#>        }
#>        else {
#>            base::invisible()
#>        }
#>    })
#> 2: tryCatchList(expr, classes, parentenv, handlers)
#> 3: tryCatchOne(tryCatchList(expr, names[-nh], parentenv, handlers[-nh]), 
#>        names[nh], parentenv, handlers[[nh]])
#> 4: doTryCatch(return(expr), name, parentenv, handler)
#> 5: tryCatchList(expr, names[-nh], parentenv, handlers[-nh])
#> 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
#> 7: doTryCatch(return(expr), name, parentenv, handler)
#> 8: base::withCallingHandlers({
#>        NULL
#>        base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp [...]
#>            base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>            quote = TRUE), file = "/tmp/RtmpWxJiB4/callr-res-29ba6f0677ca", 
#>            compress = FALSE)
#>        base::flush(base::stdout())
#>        base::flush(base::stderr())
#>        NULL
#>        base::invisible()
#>    }, error = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, interrupt = function(e) {
#>        {
#>            callr_data <- base::as.environment("tools:callr")$`__callr_data__`
#>            err <- callr_data$err
#>            if (FALSE) {
#>                base::assign(".Traceback", base::.traceback(4), envir = callr_data)
#>                utils::dump.frames("__callr_dump__")
#>                base::assign(".Last.dump", .GlobalEnv$`__callr_dump__`, 
#>                    envir = callr_data)
#>                base::rm("__callr_dump__", envir = .GlobalEnv)
#>            }
#>            e <- err$process_call(e)
#>            e2 <- err$new_error("error in callr subprocess")
#>            class <- base::class
#>            class(e2) <- base::c("callr_remote_error", class(e2))
#>            e2 <- err$add_trace_back(e2)
#>            cut <- base::which(e2$trace$scope == "global")[1]
#>            if (!base::is.na(cut)) {
#>                e2$trace <- e2$trace[-(1:cut), ]
#>            }
#>            base::saveRDS(base::list("error", e2, e), file = base::paste0("/tm [...]
#>                ".error"))
#>        }
#>    }, callr_message = function(e) {
#>        base::try(base::signalCondition(e))
#>    })
#> 9: base::saveRDS(base::do.call(base::do.call, base::c(base::readRDS("/tmp/Rtm [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE), file = "/tmp/RtmpWxJiB4/callr-res-29ba6f0677ca", 
#>        compress = FALSE)
#> 10: base::do.call(base::do.call, base::c(base::readRDS("/tmp/RtmpWxJiB4/callr- [...]
#>        base::list(envir = .GlobalEnv, quote = TRUE)), envir = .GlobalEnv, 
#>        quote = TRUE)
#> 11: (function (what, args, quote = FALSE, envir = parent.frame()) 
#>    {
#>        if (!is.list(args)) 
#>            stop("second argument must be a list")
#>        if (quote) 
#>            args <- lapply(args, enquote)
#>        .Internal(do.call(what, args, envir))
#>    })(base::quote(function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    }), base::quote(list(input = "used-aidi_reprex.R")), envir = base::quote(< [...]
#>        quote = base::quote(TRUE))
#> 12: (function (input) 
#>    {
#>        rmarkdown::render(input, quiet = TRUE, envir = globalenv(), 
#>            encoding = "UTF-8")
#>    })(input = base::quote("used-aidi_reprex.R"))
#> 13: rmarkdown::render(input, quiet = TRUE, envir = globalenv(), encoding = "UTF-8")
#> 14: knitr::knit(knit_input, knit_output, envir = envir, quiet = quiet)
#> 15: process_file(text, output)
#> 16: xfun:::handle_error(withCallingHandlers(if (tangle) process_tangle(group)  [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#>        function(e, loc) {
#>            setwd(wd)
#>            write_utf8(res, output %n% stdout())
#>            paste0("\nQuitting from lines ", loc)
#>        }, if (labels[i] != "") sprintf(" [%s]", labels[i]), get_loc)
#> 17: withCallingHandlers(expr, error = function(e) {
#>        loc = if (is.function(fun)) 
#>            trimws(fun(label))
#>        else ""
#>        if (loc != "") 
#>            loc = sprintf(" at lines %s", loc)
#>        message(one_string(handler(e, loc)))
#>    })
#> 18: withCallingHandlers(if (tangle) process_tangle(group) else process_group(g [...]
#>        error = function(e) if (xfun::pkg_available("rlang", "1.0.0")) rlang:: [...]
#> 19: process_group(group)
#> 20: call_block(x)
#> iinla: Giving up and returning last successfully obtained result for diagnostic purposes.
#> Error in track[[1]]: subscript out of bounds

# session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.4.0 (2024-04-24)
#>  os       Ubuntu 23.10
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language en_GB:en
#>  collate  en_GB.UTF-8
#>  ctype    en_GB.UTF-8
#>  tz       Europe/London
#>  date     2024-05-07
#>  pandoc   3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package      * version     date (UTC) lib source
#>  cachem         1.0.8       2023-05-01 [1] CRAN (R 4.4.0)
#>  class          7.3-22      2023-05-03 [2] CRAN (R 4.4.0)
#>  classInt       0.4-10      2023-09-05 [1] CRAN (R 4.4.0)
#>  cli            3.6.2       2023-12-11 [1] CRAN (R 4.4.0)
#>  colorspace     2.1-0       2023-01-23 [1] CRAN (R 4.4.0)
#>  DBI            1.2.2       2024-02-16 [1] CRAN (R 4.4.0)
#>  devtools       2.4.5       2022-10-11 [1] CRAN (R 4.4.0)
#>  digest         0.6.35      2024-03-11 [1] CRAN (R 4.4.0)
#>  dplyr          1.1.4       2023-11-17 [1] CRAN (R 4.4.0)
#>  e1071          1.7-14      2023-12-06 [1] CRAN (R 4.4.0)
#>  ellipsis       0.3.2       2021-04-29 [1] CRAN (R 4.4.0)
#>  evaluate       0.23        2023-11-01 [1] CRAN (R 4.4.0)
#>  fansi          1.0.6       2023-12-08 [1] CRAN (R 4.4.0)
#>  fastmap        1.1.1       2023-02-24 [1] CRAN (R 4.4.0)
#>  fmesher      * 0.1.5.9004  2024-04-29 [1] local (/home/flindgre/git/fmesher)
#>  fs             1.6.4       2024-04-25 [1] CRAN (R 4.4.0)
#>  generics       0.1.3       2022-07-05 [1] CRAN (R 4.4.0)
#>  ggplot2      * 3.5.1       2024-04-23 [1] CRAN (R 4.4.0)
#>  glue           1.7.0       2024-01-09 [1] CRAN (R 4.4.0)
#>  gtable         0.3.5       2024-04-22 [1] CRAN (R 4.4.0)
#>  htmltools      0.5.8.1     2024-04-04 [1] CRAN (R 4.4.0)
#>  htmlwidgets    1.6.4       2023-12-06 [1] CRAN (R 4.4.0)
#>  httpuv         1.6.15      2024-03-26 [1] CRAN (R 4.4.0)
#>  INLA           24.05.01-1  2024-05-01 [1] local
#>  inlabru      * 2.10.1.9005 2024-04-30 [1] local
#>  KernSmooth     2.23-22     2023-07-10 [2] CRAN (R 4.4.0)
#>  knitr          1.46        2024-04-06 [1] CRAN (R 4.4.0)
#>  later          1.3.2       2023-12-06 [1] CRAN (R 4.4.0)
#>  lattice        0.22-6      2024-03-20 [2] CRAN (R 4.4.0)
#>  lifecycle      1.0.4       2023-11-07 [1] CRAN (R 4.4.0)
#>  magrittr       2.0.3       2022-03-30 [1] CRAN (R 4.4.0)
#>  Matrix         1.7-0       2024-03-22 [2] CRAN (R 4.4.0)
#>  MatrixModels   0.5-3       2023-11-06 [1] CRAN (R 4.4.0)
#>  memoise        2.0.1       2021-11-26 [1] CRAN (R 4.4.0)
#>  mime           0.12        2021-09-28 [1] CRAN (R 4.4.0)
#>  miniUI         0.1.1.1     2018-05-18 [1] CRAN (R 4.4.0)
#>  munsell        0.5.1       2024-04-01 [1] CRAN (R 4.4.0)
#>  pillar         1.9.0       2023-03-22 [1] CRAN (R 4.4.0)
#>  pkgbuild       1.4.4       2024-03-17 [1] CRAN (R 4.4.0)
#>  pkgconfig      2.0.3       2019-09-22 [1] CRAN (R 4.4.0)
#>  pkgload        1.3.4       2024-01-16 [1] CRAN (R 4.4.0)
#>  plyr           1.8.9       2023-10-02 [1] CRAN (R 4.4.0)
#>  profvis        0.3.8       2023-05-02 [1] CRAN (R 4.4.0)
#>  promises       1.3.0       2024-04-05 [1] CRAN (R 4.4.0)
#>  proxy          0.4-27      2022-06-09 [1] CRAN (R 4.4.0)
#>  purrr          1.0.2       2023-08-10 [1] CRAN (R 4.4.0)
#>  R.cache        0.16.0      2022-07-21 [1] CRAN (R 4.4.0)
#>  R.methodsS3    1.8.2       2022-06-13 [1] CRAN (R 4.4.0)
#>  R.oo           1.26.0      2024-01-24 [1] CRAN (R 4.4.0)
#>  R.utils        2.12.3      2023-11-18 [1] CRAN (R 4.4.0)
#>  R6             2.5.1       2021-08-19 [1] CRAN (R 4.4.0)
#>  Rcpp           1.0.12      2024-01-09 [1] CRAN (R 4.4.0)
#>  remotes        2.5.0       2024-03-17 [1] CRAN (R 4.4.0)
#>  reprex         2.1.0       2024-01-11 [1] CRAN (R 4.4.0)
#>  rlang          1.1.3       2024-01-10 [1] CRAN (R 4.4.0)
#>  rmarkdown      2.26        2024-03-05 [1] CRAN (R 4.4.0)
#>  rstudioapi     0.16.0      2024-03-24 [1] CRAN (R 4.4.0)
#>  scales         1.3.0       2023-11-28 [1] CRAN (R 4.4.0)
#>  sessioninfo    1.2.2       2021-12-06 [1] CRAN (R 4.4.0)
#>  sf           * 1.0-16      2024-03-24 [1] CRAN (R 4.4.0)
#>  shiny          1.8.1.1     2024-04-02 [1] CRAN (R 4.4.0)
#>  sp             2.1-3       2024-01-30 [1] CRAN (R 4.4.0)
#>  stringi        1.8.3       2023-12-11 [1] CRAN (R 4.4.0)
#>  stringr        1.5.1       2023-11-14 [1] CRAN (R 4.4.0)
#>  styler         1.10.3      2024-04-07 [1] CRAN (R 4.4.0)
#>  tibble         3.2.1       2023-03-20 [1] CRAN (R 4.4.0)
#>  tidyselect     1.2.1       2024-03-11 [1] CRAN (R 4.4.0)
#>  units          0.8-5       2023-11-28 [1] CRAN (R 4.4.0)
#>  urlchecker     1.0.1       2021-11-30 [1] CRAN (R 4.4.0)
#>  usethis        2.2.3       2024-02-19 [1] CRAN (R 4.4.0)
#>  utf8           1.2.4       2023-10-22 [1] CRAN (R 4.4.0)
#>  vctrs          0.6.5       2023-12-01 [1] CRAN (R 4.4.0)
#>  withr          3.0.0       2024-01-16 [1] CRAN (R 4.4.0)
#>  xfun           0.43        2024-03-25 [1] CRAN (R 4.4.0)
#>  xtable         1.8-4       2019-04-21 [1] CRAN (R 4.4.0)
#>  yaml           2.3.8       2023-12-11 [1] CRAN (R 4.4.0)
#> 
#>  [1] /home/flindgre/R/x86_64-pc-linux-gnu-library/4.4
#>  [2] /home/flindgre/local/R-4.4.0/lib/R/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Created on 2024-05-07 with reprex v2.1.0

@finnlindgren
Copy link
Collaborator

The reprex::reprex() run does give the same warnings as for you. But when I run manually, the first warning,

fit <- bru(lik,
           components = cmp
)
#> Warning in pred_eps - pred0: Recycling array of length 1 in array-vector arithmetic is deprecated.
#>   Use c() or as.vector() instead.

is missing for me. I suspect this is related to the tricky issue of warnings/error capture when running inla in a tryCatch block; I'll need to look into that further, to make sure all warnings are displayed at all times.

@finnlindgren
Copy link
Collaborator

finnlindgren commented May 7, 2024 •
edited
Loading

For reference, here's the reprex with the loglambda(rep(1, nrow(.data.))) workaround, and added result output for all bru() runs:

library(inlabru)
#> Loading required package: fmesher
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.7.1, PROJ 9.2.1; sf_use_s2() is TRUE
library(ggplot2)
set.seed(183)

# region where points are observable defined as a
# square polygon with area 9 and stored as an sf polygon object.
loc.d <- 3 * cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0))
loc <- st_polygon(list(loc.d))
Omega <- st_sf(st_sfc(loc))

# set true intensity
true_lambda <- 10

# simulate the number of points in the pattern
N <- rpois(1, true_lambda * 9)

# simulate random locations using st_sample()
pp <- st_sf(
  geometry = st_sample(Omega, N)
)

# define integration scheme
# since homogeneous just 1 location with weight equal
# to area of the observation window
ips <- st_sf(
  geometry = st_sample(Omega, 1)
)
ips$weight <- 9

# components - just the log intensity parameter
cmp <- ~ 0 + loglambda(rep(1, nrow(.data.)))

# fit without compressed cp likelihood
# runs but with a warning message
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = FALSE
  )
)

fit <- bru(lik,
           components = cmp
)

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(rep(1, nrow(.data.))), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.43, Running = 0.186, Post = 0.0179, Total = 0.634 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -609.29
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -182.79
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -340.04
#> Effective number of parameters .................: 7.80
#> 
#> Marginal log-Likelihood:  116.14 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

# this should be around 10
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# fit with compressed cp likelihood
# gives an error
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = TRUE
  )
)

fit <- bru(lik,
           components = cmp
)

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(rep(1, nrow(.data.))), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.2, Running = 0.187, Post = 0.0232, Total = 0.411 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -410.23
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -410.44
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 157.77
#> Effective number of parameters .................: 6.80
#> 
#> Marginal log-Likelihood:  -211.04 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
exp(fit$summary.fixed$mean)
#> [1] 10.16656



# try an integration scheme with 9 points
#########################################

ips <- st_sf(
  geometry = st_sample(Omega, 9)
)
ips$weight <- 1

# non-compressed cp likelihood
# still works
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = FALSE
  )
)

fit <- bru(lik,
           components = cmp
)

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(rep(1, nrow(.data.))), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.152, Running = 0.187, Post = 0.0109, Total = 0.35 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -241.72
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: 0.997
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: -230.43
#> Effective number of parameters .................: 11.28
#> 
#> Marginal log-Likelihood:  116.14 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# compressed cp likelihood
# still gives an error
lik <- like(
  data = pp,
  family = "cp",
  formula = geometry ~ .,
  ips = ips,
  options = list(
    bru_compress_cp = TRUE
  )
)

fit <- bru(lik,
           components = cmp
)

summary(fit)
#> inlabru version: 2.10.1.9005
#> INLA version: 24.05.01-1
#> Components:
#> loglambda: main = linear(rep(1, nrow(.data.))), group = exchangeable(1L), replicate = iid(1L)
#> Likelihoods:
#>   Family: 'cp'
#>     Data class: 'sf', 'data.frame'
#>     Predictor: geometry ~ .
#> Time used:
#>     Pre = 0.147, Running = 0.176, Post = 0.0132, Total = 0.335 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> loglambda 2.319 0.104      2.115    2.319      2.523 2.319   0
#> 
#> Deviance Information Criterion (DIC) ...............: -42.66
#> Deviance Information Criterion (DIC, saturated) ....: NA
#> Effective number of parameters .....................: -226.66
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 267.37
#> Effective number of parameters .................: 10.28
#> 
#> Marginal log-Likelihood:  -211.04 
#>  is computed 
#> Posterior summaries for the linear predictor and the fitted values are computed
#> (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
exp(fit$summary.fixed$mean)
#> [1] 10.16655

# session info
devtools::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.4.0 (2024-04-24)
#>  os       Ubuntu 23.10
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language en_GB:en
#>  collate  en_GB.UTF-8
#>  ctype    en_GB.UTF-8
#>  tz       Europe/London
#>  date     2024-05-07
#>  pandoc   3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package      * version     date (UTC) lib source
#>  cachem         1.0.8       2023-05-01 [1] CRAN (R 4.4.0)
#>  class          7.3-22      2023-05-03 [2] CRAN (R 4.4.0)
#>  classInt       0.4-10      2023-09-05 [1] CRAN (R 4.4.0)
#>  cli            3.6.2       2023-12-11 [1] CRAN (R 4.4.0)
#>  colorspace     2.1-0       2023-01-23 [1] CRAN (R 4.4.0)
#>  DBI            1.2.2       2024-02-16 [1] CRAN (R 4.4.0)
#>  devtools       2.4.5       2022-10-11 [1] CRAN (R 4.4.0)
#>  digest         0.6.35      2024-03-11 [1] CRAN (R 4.4.0)
#>  dplyr          1.1.4       2023-11-17 [1] CRAN (R 4.4.0)
#>  e1071          1.7-14      2023-12-06 [1] CRAN (R 4.4.0)
#>  ellipsis       0.3.2       2021-04-29 [1] CRAN (R 4.4.0)
#>  evaluate       0.23        2023-11-01 [1] CRAN (R 4.4.0)
#>  fansi          1.0.6       2023-12-08 [1] CRAN (R 4.4.0)
#>  fastmap        1.1.1       2023-02-24 [1] CRAN (R 4.4.0)
#>  fmesher      * 0.1.5.9004  2024-04-29 [1] local (/home/flindgre/git/fmesher)
#>  fs             1.6.4       2024-04-25 [1] CRAN (R 4.4.0)
#>  generics       0.1.3       2022-07-05 [1] CRAN (R 4.4.0)
#>  ggplot2      * 3.5.1       2024-04-23 [1] CRAN (R 4.4.0)
#>  glue           1.7.0       2024-01-09 [1] CRAN (R 4.4.0)
#>  gtable         0.3.5       2024-04-22 [1] CRAN (R 4.4.0)
#>  htmltools      0.5.8.1     2024-04-04 [1] CRAN (R 4.4.0)
#>  htmlwidgets    1.6.4       2023-12-06 [1] CRAN (R 4.4.0)
#>  httpuv         1.6.15      2024-03-26 [1] CRAN (R 4.4.0)
#>  INLA           24.05.01-1  2024-05-01 [1] local
#>  inlabru      * 2.10.1.9005 2024-04-30 [1] local
#>  KernSmooth     2.23-22     2023-07-10 [2] CRAN (R 4.4.0)
#>  knitr          1.46        2024-04-06 [1] CRAN (R 4.4.0)
#>  later          1.3.2       2023-12-06 [1] CRAN (R 4.4.0)
#>  lattice        0.22-6      2024-03-20 [2] CRAN (R 4.4.0)
#>  lifecycle      1.0.4       2023-11-07 [1] CRAN (R 4.4.0)
#>  magrittr       2.0.3       2022-03-30 [1] CRAN (R 4.4.0)
#>  Matrix         1.7-0       2024-03-22 [2] CRAN (R 4.4.0)
#>  MatrixModels   0.5-3       2023-11-06 [1] CRAN (R 4.4.0)
#>  memoise        2.0.1       2021-11-26 [1] CRAN (R 4.4.0)
#>  mime           0.12        2021-09-28 [1] CRAN (R 4.4.0)
#>  miniUI         0.1.1.1     2018-05-18 [1] CRAN (R 4.4.0)
#>  munsell        0.5.1       2024-04-01 [1] CRAN (R 4.4.0)
#>  pillar         1.9.0       2023-03-22 [1] CRAN (R 4.4.0)
#>  pkgbuild       1.4.4       2024-03-17 [1] CRAN (R 4.4.0)
#>  pkgconfig      2.0.3       2019-09-22 [1] CRAN (R 4.4.0)
#>  pkgload        1.3.4       2024-01-16 [1] CRAN (R 4.4.0)
#>  plyr           1.8.9       2023-10-02 [1] CRAN (R 4.4.0)
#>  profvis        0.3.8       2023-05-02 [1] CRAN (R 4.4.0)
#>  promises       1.3.0       2024-04-05 [1] CRAN (R 4.4.0)
#>  proxy          0.4-27      2022-06-09 [1] CRAN (R 4.4.0)
#>  purrr          1.0.2       2023-08-10 [1] CRAN (R 4.4.0)
#>  R.cache        0.16.0      2022-07-21 [1] CRAN (R 4.4.0)
#>  R.methodsS3    1.8.2       2022-06-13 [1] CRAN (R 4.4.0)
#>  R.oo           1.26.0      2024-01-24 [1] CRAN (R 4.4.0)
#>  R.utils        2.12.3      2023-11-18 [1] CRAN (R 4.4.0)
#>  R6             2.5.1       2021-08-19 [1] CRAN (R 4.4.0)
#>  Rcpp           1.0.12      2024-01-09 [1] CRAN (R 4.4.0)
#>  remotes        2.5.0       2024-03-17 [1] CRAN (R 4.4.0)
#>  reprex         2.1.0       2024-01-11 [1] CRAN (R 4.4.0)
#>  rlang          1.1.3       2024-01-10 [1] CRAN (R 4.4.0)
#>  rmarkdown      2.26        2024-03-05 [1] CRAN (R 4.4.0)
#>  rstudioapi     0.16.0      2024-03-24 [1] CRAN (R 4.4.0)
#>  scales         1.3.0       2023-11-28 [1] CRAN (R 4.4.0)
#>  sessioninfo    1.2.2       2021-12-06 [1] CRAN (R 4.4.0)
#>  sf           * 1.0-16      2024-03-24 [1] CRAN (R 4.4.0)
#>  shiny          1.8.1.1     2024-04-02 [1] CRAN (R 4.4.0)
#>  sp             2.1-3       2024-01-30 [1] CRAN (R 4.4.0)
#>  stringi        1.8.3       2023-12-11 [1] CRAN (R 4.4.0)
#>  stringr        1.5.1       2023-11-14 [1] CRAN (R 4.4.0)
#>  styler         1.10.3      2024-04-07 [1] CRAN (R 4.4.0)
#>  tibble         3.2.1       2023-03-20 [1] CRAN (R 4.4.0)
#>  tidyselect     1.2.1       2024-03-11 [1] CRAN (R 4.4.0)
#>  units          0.8-5       2023-11-28 [1] CRAN (R 4.4.0)
#>  urlchecker     1.0.1       2021-11-30 [1] CRAN (R 4.4.0)
#>  usethis        2.2.3       2024-02-19 [1] CRAN (R 4.4.0)
#>  utf8           1.2.4       2023-10-22 [1] CRAN (R 4.4.0)
#>  vctrs          0.6.5       2023-12-01 [1] CRAN (R 4.4.0)
#>  withr          3.0.0       2024-01-16 [1] CRAN (R 4.4.0)
#>  xfun           0.43        2024-03-25 [1] CRAN (R 4.4.0)
#>  xtable         1.8-4       2019-04-21 [1] CRAN (R 4.4.0)
#>  yaml           2.3.8       2023-12-11 [1] CRAN (R 4.4.0)
#> 
#>  [1] /home/flindgre/R/x86_64-pc-linux-gnu-library/4.4
#>  [2] /home/flindgre/local/R-4.4.0/lib/R/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────

Created on 2024-05-07 with reprex v2.1.0

@ASeatonSpatial
Copy link
Collaborator Author

Using loglambda(rep(1, nrow(.data.))) also works for me. I am writing an example for the textbook we are writing, this will be a good opportunity to point out that the (1) shorthand is not infallible and how to do it safely.

I also don't see the warning when I run the code in Rstudio but they do appear in the reprex.

@finnlindgren
Copy link
Collaborator

I may find a way to make the workaround unnecessary. The main issue is that it needs to extract information about how long the predictor vector is supposed to be; for normal models, that's nrow(.data.), but not for models with aggregation. For those, it would internally be nrow(response_data), but only if that's a data.frame. So the main obstacle to getting this to work is to sort out how to safely and accurately extract the needed information.

@finnlindgren finnlindgren changed the title Error fitting homogeneous Poisson process with compressed cp likelihood Automated predictor size detection needed for intercept-only models (Was: Error fitting homogeneous Poisson process with compressed cp likelihood) May 17, 2024
@finnlindgren finnlindgren self-assigned this May 17, 2024
@finnlindgren
Copy link
Collaborator

This has been partially solved by aab9e90, that introduced bru_response_size() and further logic for expanding scalar effects to full vectors in unambiguous cases, that may include inspecting the response_data contents to determine the required predictor size. However, there are further cases where it may be impossible to implement such logic, such as in predict() and generate(), where response data is inherently unavailable, and nrow(.data.) is the only available information, but this may not be the predictor size. Such cases include aggregation predictors.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees

@finnlindgren finnlindgren

Labels
enhancement
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@finnlindgren @ASeatonSpatial