INLA crash when fitting the Mathematical model
19 views
Skip to first unread message
Barry Thierno Souleymane
Jul 24, 2025, 6:24:55 AMJul 24
to R-inla discussion group, hr...@r-inla.org, Finn Lindgren
Hi INLA Team,
I am running the codes below.
Everything looks good, except for the part about INLA .
After fitting the model using INLA formula, the following results were obtained
Please, I need your help.
ERRORS
Error in inla.core(formula = formula, family = family, contrasts = contrasts, :
In f(region_car): 'covariate' must match 'values', and both must either be 'numeric', or 'factor'/'character'.
*** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=1/2
Error in inla.core(formula = formula, family = family, contrasts = contrasts, :
In f(region_car): 'covariate' must match 'values', and both must either be 'numeric', or 'factor'/'character'.
*** inla.core.safe: inla.program has crashed: rerun to get better initial values. try=2/2
Error in inla.core(formula = formula, family = family, contrasts = contrasts, :
In f(region_car): 'covariate' must match 'values', and both must either be 'numeric', or 'factor'/'character'.
Error in inla.core.safe(formula = formula, family = family, contrasts = contrasts, :
*** Fail to get good enough initial values. Maybe it is due to something else.The codes I have used are below:
###MY CODES
###1. Étape : To Install packages
install.packages(c("deSolve", "BayesianTools", "ggplot2"))
#####2. Étape : Simuler les données avec un modèle SIR
install.packages("ggplot2", repos = "https://cloud.r-project.org/")
install.packages("lme4", type = "source")
install.packages("spdep")
install.packages("ggplot2")
library(deSolve)
library(BayesianTools)
library(ggplot2)
library(MASS) # pour mvrnorm
library(spdep) # pour créer la matrice W
rm(list = ls())
# Modèle SIR
slir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dS <- omega-beta * S * I-(enn+zeta)*S
dL <- beta * S * I-(enn+zeta+teta )* L
dI <- teta * L-(lamda+enn+zeta+teta+eta )* I
dR <- lamda-(enn+zeta)* R
list(c(dS, dL, dI, dR))
})
}
### multi-régions data Simulation with random effects
set.seed(123)
N <- 8 # nombre de régions
times <- seq(0, 30, by = 1)
# Random Effets
# Hiérarchy Parameters
mu_omega <- log(0.1); sigma_v_omega <- 0.2
mu_beta <- log(0.3); sigma_v_beta <- 0.3
mu_enn <- log(0.2); sigma_v_enn <- 0.1
mu_zeta <- log(0.2); sigma_v_zeta <- 0.1
mu_teta <- log(0.1); sigma_v_teta <- 0.2
mu_eta <- log(0.2); sigma_v_eta <- 0.2
mu_lamda <- log(0.3); sigma_v_lamda <- 0.3
# Non spatial random Effects
v_omega <- rnorm(N, 0, sigma_v_omega)
v_beta <- rnorm(N, 0, sigma_v_beta)
v_enn <- rnorm(N, 0, sigma_v_enn)
v_zeta <- rnorm(N, 0, sigma_v_zeta)
v_teta <- rnorm(N, 0, sigma_v_teta)
v_eta <- rnorm(N, 0, sigma_v_eta)
v_lamda <- rnorm(N, 0, sigma_v_lamda)
### spatial random effects
#install.packages(c("sp", "sf", "rgeos", "rgdal", "spdep"))
#install.packages("spdep")
#library(spdep)
#install.packages("sp")
coords <- cbind(1:N, rep(1, N)) # disposition en ligne
neighbors <- dnearneigh(coords, 0, 1.5) # voisin immédiat
W <- nb2mat(neighbors, style = "B", zero.policy = TRUE)
D <- diag(rowSums(W))
rho <- 0.9
Q <- D - rho * W
# Generation of the spatial effect u ~ N(0, Q^-1)
tau <- 1
Sigma_u <- ginv(Q) / tau
u_omega <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_beta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_enn <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_zeta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_teta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_eta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_lamda <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
# regional final Parameters
log_omega <- mu_omega + u_omega + v_omega
log_beta <- mu_beta + u_beta + v_beta
log_enn <- mu_enn + u_enn + v_enn
log_zeta <- mu_zeta + u_zeta + v_zeta
log_teta <- mu_teta + u_teta + v_teta
log_eta <- mu_eta + u_eta + v_eta
log_lamda <- mu_lamda + u_lamda + v_lamda
omega_i <- exp(log_omega)
beta_i <- exp(log_beta)
enn_i <- exp(log_enn)
zeta_i <- exp(log_zeta)
teta_i <- exp(log_teta)
eta_i <- exp(log_eta)
lamda_i <- exp(log_lamda)
##A measurement variance is also set for each region
sigma_i <- rep(0.01, N) # or variable by region
init_state <- c(S = 0.97, L=0.02, I = 0.01, R = 0)
### Simulation for regions
all_data <- list()
for (i in 1:N) {
params <- c(omega= omega_i[i], eta= eta_i[i], zeta= zeta_i[i], teta= teta_i[i], beta = beta_i[i], enn= enn_i[i], lamda = lamda_i[i])
sim <- ode(y = init_state, times = times, func = slir_model, parms = params)
sim_df <- as.data.frame(sim)
# Add observation noise
sim_df$I_obs <- sim_df$I + rnorm(length(sim_df$I), mean = 0, sd = sigma_i[i])
sim_df$region <- i
sim_df$omega <- omega_i[i]
sim_df$beta <- beta_i[i]
sim_df$enn<- enn_i[i]
sim_df$zeta <- zeta_i[i]
sim_df$teta <- teta_i[i]
sim_df$eta <- eta_i[i]
sim_df$lamda <- lamda_i[i]
all_data[[i]] <- sim_df
}
# Merge data
sim_data <- do.call(rbind, all_data)
#head(sim_data)
### 6. Visualisation
ggplot(sim_data, aes(x = time, y = I_obs, color = factor(region))) +
geom_point() +
geom_line(aes(y = I), linetype = "dashed") +
labs(title = "slir multi-régions avec effets aléatoires", y = "Infectés", x = "Temps") +
theme_minimal()
####Step: Bayesian modeling using INLA
#install.packages("INLA", repos = c(INLA = "https://inla.r-inla-download.org/R/stable"), type = "binary")
#library(INLA)
# Simulated data
data <- sim_data
data$log_I_obs <-log(pmax(data$I_obs, 1e-6)) # pour éviter log(0)
#head(data)
#head(data$log_I_obs)
data$region_iid <- data$region
data$region_car <- data$region
#head(data$region_car)
data$region_car <- as.character(data$region_car)
#class(data$region_car)
# Replace with the path of your .shp file
shp <- st_read("C:/Users/BARRY/Downloads/gin_admbnda_ocha_fis/gin_admbnda_adm1_ocha.shp")
#plot(shp)
#names(shp)
#names(data)
#data$region
#data$region_car
#data$region_iid
#table(region_name)
#unique(sim_data$region)
#unique(shp$region_car)
#Create the lookup table:
region_lut <- data.frame(
region = 1:8,
region_name = c("Boke", "Conakry", "Faranah", "Kankan",
"Kindia", "Labe", "Mamou", "Nzerekore")
)
names(sim_data)
names(region_lut)
#Join this table to your data:
data$region
region_lut$region
sim_data <- merge(data, region_lut, by = "region")
#names(sim_data)
#Merge with the shapefile (without geometry):
#names(shp)
#data$region_car
#data$region_iid
#data$admin1Name
#data$shp$admin1Pcod
shp$region_car <- shp$admin1Name # names of regions (character)
shp$region_iid <- shp$admin1Name # even if these are the names
shp_data <- st_drop_geometry(shp)[, c("region_car", "region_iid")]
#shp_data$region_car
#shp_data$region_iid
#sim_data$region_name
#sim_data$region_car
dataF <- merge(sim_data, shp_data, by.x = "region_name", by.y = "region_car")
#names(dataF)
#dataF$region_car
#dataF$region_iid.x
#dataF$region_iid.y
#Add the required columns for INLA:
#dataF$region_name
#dataF$region_iid.y
dataF$region_car <- as.character(dataF$region_name) # Pour l'effet spatial
dataF$region_iid <- as.character(dataF$region_iid.y) # Si pas déjà en caractère
#names(dataF)
#Option: With ggplot2 for more aesthetic
#library(ggplot2)
#ggplot(data = shp) +
# geom_sf(fill = "lightblue", color = "black") +
# theme_minimal()
# Convert to Spatial object if necessary
#shp_sp <- as(shp, "Spatial")
###############################################################
# Create neighbors from polygons
nb <- poly2nb(shp, row.names = shp$region_car)
nb_inla <- nb2INLA("adjacency.graph", nb)
# Model with spatial effect (CAR) + random effect (iid)
#data$region_iid
#data$region_car <- as.character(data$region) # Si data$region va de 1 à 8
#data$region_iid <- as.character(data$region) # Si data$region va de 1 à 8
formula <- log_I_obs ~ 1 +
f(region_car, model = "besag", graph = "adjacency.graph") +
f(region_iid, model = "iid")
# Model adjustment
#names(sim_data)
#names(data)
#data$log_I_obs
head(dataF)
result <- inla(
formula,
family = "gaussian",
data =dataF,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE)
)
# Sommary
summary(result)
summary(result)$random$region_car
car_summary <- result$summary.random$region_car
ggplot(car_summary, aes(x = ID, y = mean)) +
geom_point(color = "blue") +
geom_errorbar(aes(ymin = `0.025quant`, ymax = `0.975quant`), width = 0.2) +
labs(title = "Effet spatial estimé par région", y = "Effet CAR", x = "Région") +
theme_minimal()
install.packages(c("deSolve", "BayesianTools", "ggplot2"))
#####2. Étape : Simuler les données avec un modèle SIR
install.packages("ggplot2", repos = "https://cloud.r-project.org/")
install.packages("lme4", type = "source")
install.packages("spdep")
install.packages("ggplot2")
library(deSolve)
library(BayesianTools)
library(ggplot2)
library(MASS) # pour mvrnorm
library(spdep) # pour créer la matrice W
rm(list = ls())
# Modèle SIR
slir_model <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dS <- omega-beta * S * I-(enn+zeta)*S
dL <- beta * S * I-(enn+zeta+teta )* L
dI <- teta * L-(lamda+enn+zeta+teta+eta )* I
dR <- lamda-(enn+zeta)* R
list(c(dS, dL, dI, dR))
})
}
### multi-régions data Simulation with random effects
set.seed(123)
N <- 8 # nombre de régions
times <- seq(0, 30, by = 1)
# Random Effets
# Hiérarchy Parameters
mu_omega <- log(0.1); sigma_v_omega <- 0.2
mu_beta <- log(0.3); sigma_v_beta <- 0.3
mu_enn <- log(0.2); sigma_v_enn <- 0.1
mu_zeta <- log(0.2); sigma_v_zeta <- 0.1
mu_teta <- log(0.1); sigma_v_teta <- 0.2
mu_eta <- log(0.2); sigma_v_eta <- 0.2
mu_lamda <- log(0.3); sigma_v_lamda <- 0.3
# Non spatial random Effects
v_omega <- rnorm(N, 0, sigma_v_omega)
v_beta <- rnorm(N, 0, sigma_v_beta)
v_enn <- rnorm(N, 0, sigma_v_enn)
v_zeta <- rnorm(N, 0, sigma_v_zeta)
v_teta <- rnorm(N, 0, sigma_v_teta)
v_eta <- rnorm(N, 0, sigma_v_eta)
v_lamda <- rnorm(N, 0, sigma_v_lamda)
### spatial random effects
#install.packages(c("sp", "sf", "rgeos", "rgdal", "spdep"))
#install.packages("spdep")
#library(spdep)
#install.packages("sp")
coords <- cbind(1:N, rep(1, N)) # disposition en ligne
neighbors <- dnearneigh(coords, 0, 1.5) # voisin immédiat
W <- nb2mat(neighbors, style = "B", zero.policy = TRUE)
D <- diag(rowSums(W))
rho <- 0.9
Q <- D - rho * W
# Generation of the spatial effect u ~ N(0, Q^-1)
tau <- 1
Sigma_u <- ginv(Q) / tau
u_omega <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_beta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_enn <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_zeta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_teta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_eta <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
u_lamda <- mvrnorm(1, mu = rep(0, N), Sigma = Sigma_u)
# regional final Parameters
log_omega <- mu_omega + u_omega + v_omega
log_beta <- mu_beta + u_beta + v_beta
log_enn <- mu_enn + u_enn + v_enn
log_zeta <- mu_zeta + u_zeta + v_zeta
log_teta <- mu_teta + u_teta + v_teta
log_eta <- mu_eta + u_eta + v_eta
log_lamda <- mu_lamda + u_lamda + v_lamda
omega_i <- exp(log_omega)
beta_i <- exp(log_beta)
enn_i <- exp(log_enn)
zeta_i <- exp(log_zeta)
teta_i <- exp(log_teta)
eta_i <- exp(log_eta)
lamda_i <- exp(log_lamda)
##A measurement variance is also set for each region
sigma_i <- rep(0.01, N) # or variable by region
init_state <- c(S = 0.97, L=0.02, I = 0.01, R = 0)
### Simulation for regions
all_data <- list()
for (i in 1:N) {
params <- c(omega= omega_i[i], eta= eta_i[i], zeta= zeta_i[i], teta= teta_i[i], beta = beta_i[i], enn= enn_i[i], lamda = lamda_i[i])
sim <- ode(y = init_state, times = times, func = slir_model, parms = params)
sim_df <- as.data.frame(sim)
# Add observation noise
sim_df$I_obs <- sim_df$I + rnorm(length(sim_df$I), mean = 0, sd = sigma_i[i])
sim_df$region <- i
sim_df$omega <- omega_i[i]
sim_df$beta <- beta_i[i]
sim_df$enn<- enn_i[i]
sim_df$zeta <- zeta_i[i]
sim_df$teta <- teta_i[i]
sim_df$eta <- eta_i[i]
sim_df$lamda <- lamda_i[i]
all_data[[i]] <- sim_df
}
# Merge data
sim_data <- do.call(rbind, all_data)
#head(sim_data)
### 6. Visualisation
ggplot(sim_data, aes(x = time, y = I_obs, color = factor(region))) +
geom_point() +
geom_line(aes(y = I), linetype = "dashed") +
labs(title = "slir multi-régions avec effets aléatoires", y = "Infectés", x = "Temps") +
theme_minimal()
####Step: Bayesian modeling using INLA
#install.packages("INLA", repos = c(INLA = "https://inla.r-inla-download.org/R/stable"), type = "binary")
#library(INLA)
# Simulated data
data <- sim_data
data$log_I_obs <-log(pmax(data$I_obs, 1e-6)) # pour éviter log(0)
#head(data)
#head(data$log_I_obs)
data$region_iid <- data$region
data$region_car <- data$region
#head(data$region_car)
data$region_car <- as.character(data$region_car)
#class(data$region_car)
# Replace with the path of your .shp file
shp <- st_read("C:/Users/BARRY/Downloads/gin_admbnda_ocha_fis/gin_admbnda_adm1_ocha.shp")
#plot(shp)
#names(shp)
#names(data)
#data$region
#data$region_car
#data$region_iid
#table(region_name)
#unique(sim_data$region)
#unique(shp$region_car)
#Create the lookup table:
region_lut <- data.frame(
region = 1:8,
region_name = c("Boke", "Conakry", "Faranah", "Kankan",
"Kindia", "Labe", "Mamou", "Nzerekore")
)
names(sim_data)
names(region_lut)
#Join this table to your data:
data$region
region_lut$region
sim_data <- merge(data, region_lut, by = "region")
#names(sim_data)
#Merge with the shapefile (without geometry):
#names(shp)
#data$region_car
#data$region_iid
#data$admin1Name
#data$shp$admin1Pcod
shp$region_car <- shp$admin1Name # names of regions (character)
shp$region_iid <- shp$admin1Name # even if these are the names
shp_data <- st_drop_geometry(shp)[, c("region_car", "region_iid")]
#shp_data$region_car
#shp_data$region_iid
#sim_data$region_name
#sim_data$region_car
dataF <- merge(sim_data, shp_data, by.x = "region_name", by.y = "region_car")
#names(dataF)
#dataF$region_car
#dataF$region_iid.x
#dataF$region_iid.y
#Add the required columns for INLA:
#dataF$region_name
#dataF$region_iid.y
dataF$region_car <- as.character(dataF$region_name) # Pour l'effet spatial
dataF$region_iid <- as.character(dataF$region_iid.y) # Si pas déjà en caractère
#names(dataF)
#Option: With ggplot2 for more aesthetic
#library(ggplot2)
#ggplot(data = shp) +
# geom_sf(fill = "lightblue", color = "black") +
# theme_minimal()
# Convert to Spatial object if necessary
#shp_sp <- as(shp, "Spatial")
###############################################################
# Create neighbors from polygons
nb <- poly2nb(shp, row.names = shp$region_car)
nb_inla <- nb2INLA("adjacency.graph", nb)
# Model with spatial effect (CAR) + random effect (iid)
#data$region_iid
#data$region_car <- as.character(data$region) # Si data$region va de 1 à 8
#data$region_iid <- as.character(data$region) # Si data$region va de 1 à 8
formula <- log_I_obs ~ 1 +
f(region_car, model = "besag", graph = "adjacency.graph") +
f(region_iid, model = "iid")
# Model adjustment
#names(sim_data)
#names(data)
#data$log_I_obs
head(dataF)
result <- inla(
formula,
family = "gaussian",
data =dataF,
control.predictor = list(compute = TRUE),
control.compute = list(dic = TRUE, waic = TRUE)
)
# Sommary
summary(result)
summary(result)$random$region_car
car_summary <- result$summary.random$region_car
ggplot(car_summary, aes(x = ID, y = mean)) +
geom_point(color = "blue") +
geom_errorbar(aes(ymin = `0.025quant`, ymax = `0.975quant`), width = 0.2) +
labs(title = "Effet spatial estimé par région", y = "Effet CAR", x = "Région") +
theme_minimal()
Dr. Barry Thierno Souleymane (Ph.D)
Statisticien/Démographe
Consultant
Tel: 00224622474270
Håvard Rue
Jul 24, 2025, 11:54:47 AMJul 24
to Barry Thierno Souleymane, R-inla discussion group, Finn Lindgren
need more
> You received this message because you are subscribed to the Google Groups "R-
> inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to r-inla-discussion...@googlegroups.com.
> To view this discussion, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/CAO%2Bjm%2B7bHXgjzvf9A5hSjW1SvJ%3Dj26E_n6byMhLdfFxFo8aLbg%40mail.gmail.com
> .
--
Håvard Rue
hr...@r-inla.org
Håvard Rue
Jul 24, 2025, 11:56:12 AMJul 24
to Barry Thierno Souleymane, R-inla discussion group, Finn Lindgren
I would guess the error is because variable
region_car
in the inla - data argument, is not numeric and integers in the range
1,2,3,...n, where n is the size of the graph
On Thu, 2025-07-24 at 10:35 +0000, Barry Thierno Souleymane wrote:
region_car
in the inla - data argument, is not numeric and integers in the range
1,2,3,...n, where n is the size of the graph
On Thu, 2025-07-24 at 10:35 +0000, Barry Thierno Souleymane wrote:
Barry Thierno Souleymane
Jul 24, 2025, 3:29:24 PMJul 24
to Håvard Rue, R-inla discussion group, Finn Lindgren
Okay, thank you
Dr. Barry Thierno Souleymane (Ph.D)
Statisticien/Démographe
Consultant
Tel: 00224622474270
Reply all
Reply to author
Forward
0 new messages