# packages
library(INLA)
#> Carico il pacchetto richiesto: Matrix
#> Carico il pacchetto richiesto: foreach
#> Carico il pacchetto richiesto: parallel
#> Carico il pacchetto richiesto: sp
#> This is INLA_21.02.23 built 2021-02-22 21:11:05 UTC.
#> - See www.r-inla.org/contact-us for how to get help.
# simulate from the following model
# yi ~ N(b_0 + b_z * z, 1), i = 1, ..., 200
# z ~ N(0, 1); b_0 = 1; b_z = -1
set.seed(08032021)
n <- 200; b_0 <- 1; b_z <- -1
z <- rnorm(n)
mu <- b_0 + b_z * z
y <- rnorm(n, mu)
summary(inla(
formula = y ~ z,
data = list(y = y, z = z)
))
#>
#> Call:
#> "inla(formula = y ~ z, data = list(y = y, z = z))"
#> Time used:
#> Pre = 1.47, Running = 0.231, Post = 0.294, Total = 2
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> (Intercept) 0.746 0.071 0.606 0.746 0.886 0.746 0
#> z -0.982 0.068 -1.116 -0.982 -0.848 -0.982 0
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant
#> Precision for the Gaussian observations 1.00 0.10 0.816 0.999
#> 0.975quant mode
#> Precision for the Gaussian observations 1.21 0.993
#>
#> Expected number of effective parameters(stdev): 2.00(0.00)
#> Number of equivalent replicates : 99.94
#>
#> Marginal log-Likelihood: -302.67
# If I add an iid Gausian random effect, then I need to set the precision of the
# gaussian family to a very large and fixed value, otherwise I think there are
# confouding effects. For example
ID <- 1:n
summary(inla(
formula = y ~ z + f(ID),
data = list(y = y, z = z, ID = ID),
control.family = list(
hyper = list(
prec = list(
initial = log(10e5),
fixed = TRUE
)
)
)
))
#>
#> Call:
#> c("inla(formula = y ~ z + f(ID), data = list(y = y, z = z, ID = ID), ",
#> " control.family = list(hyper = list(prec = list(initial = log(1e+06),
#> ", " fixed = TRUE))))")
#> Time used:
#> Pre = 1.46, Running = 0.375, Post = 0.215, Total = 2.05
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> (Intercept) 0.746 0.071 0.606 0.746 0.886 0.746 0
#> z -0.982 0.068 -1.116 -0.982 -0.848 -0.982 0
#>
#> Random effects:
#> Name Model
#> ID IID model
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant 0.975quant mode
#> Precision for ID 1.00 0.10 0.816 0.999 1.21 0.993
#>
#> Expected number of effective parameters(stdev): 200.00(0.00)
#> Number of equivalent replicates : 1.00
#>
#> Marginal log-Likelihood: -302.67
# Clear ws
rm(list = ls())
# Simulate data from the following model:
# 1) Regression model
# y_i ~ Pois(lambda_i), i = 1, ..., 500
# log(lambda_i) = b_0 + b_z * z + b_x * x
# 2) Exposure model
# x = a_0 + a_z * z + eps
# 3) Error model
# w = x + u
# where
# b_0 = 1, b_z = -1 and b_x = 1;
# a_0 = 1; a_z = -1
# z ~ N(0, 1); eps ~ N(0, 1); u ~ N(0, 1)
# Constants
set.seed(08032021)
n <- 500
b_0 <- 1; b_z <- -1; b_x <- 1
a_0 <- 1; a_z <- -1
# Exposure model
z <- rnorm(n)
eps <- rnorm(n)
x <- a_0 + a_z * z + eps
# Error model
u <- rnorm(n)
w <- x + u
# Regression model
lambda <- b_0 + b_z * z + b_x * x
y <- rpois(n, exp(lambda))
# Build the Y matrix
Y <- matrix(NA_real_, nrow = 3 * n, ncol = 3)
Y[1:n, 1] <- y
Y[n + 1:n, 2] <- 0
Y[2 * n + 1:n, 3] <- w
# Reformat the data
my_data <- list(
# Regression model
Y = Y,
b_0 = c(rep(1, n), rep(NA, 2 * n)),
b_z = c(z, rep(NA, 2 * n)),
b_x = c(1:n, rep(NA, 2 * n)),
weight_x = c(rep(1, n), rep(-1, n), rep(1, n)),
idx_x = c(rep(NA, n), 1:n, 1:n),
# Exposure model
a_0 = c(rep(NA, n), rep(1, n), rep(NA, n)),
a_z = c(rep(NA, n), z, rep(NA, n))
)
# Run the INLA model
summary(inla(
formula = Y ~ -1 +
f(b_x, copy = "idx_x", hyper = list(beta = list(param = c(0, 0.1), fixed = FALSE))) +
f(idx_x, weight_x, model = "iid", values = 1:500, hyper = list(
prec = list(initial = -15, fixed = TRUE)
)) +
b_0 + b_z + a_0 + a_z,
family = c("poisson", "gaussian", "gaussian"),
data = my_data
))
#>
#> Call:
#> c("inla(formula = Y ~ -1 + f(b_x, copy = \"idx_x\", hyper = list(beta =
#> list(param = c(0, ", " 0.1), fixed = FALSE))) + f(idx_x, weight_x,
#> model = \"iid\", ", " values = 1:500, hyper = list(prec = list(initial
#> = -15, fixed = TRUE))) + ", " b_0 + b_z + a_0 + a_z, family =
#> c(\"poisson\", \"gaussian\", ", " \"gaussian\"), data = my_data)")
#> Time used:
#> Pre = 2.1, Running = 4.44, Post = 0.491, Total = 7.03
#> Fixed effects:
#> mean sd 0.025quant 0.5quant 0.975quant mode kld
#> b_0 0.937 0.080 0.780 0.939 1.087 0.940 0
#> b_z -0.937 0.067 -1.064 -0.939 -0.802 -0.942 0
#> a_0 1.053 0.062 0.931 1.053 1.175 1.053 0
#> a_z -0.989 0.060 -1.106 -0.989 -0.871 -0.989 0
#>
#> Random effects:
#> Name Model
#> idx_x IID model
#> b_x Copy
#>
#> Model hyperparameters:
#> mean sd 0.025quant 0.5quant
#> Precision for the Gaussian observations[2] 1.01 0.095 0.836 1.00
#> Precision for the Gaussian observations[3] 1.08 0.092 0.909 1.07
#> Beta for b_x 1.02 0.044 0.934 1.02
#> 0.975quant mode
#> Precision for the Gaussian observations[2] 1.21 0.993
#> Precision for the Gaussian observations[3] 1.27 1.067
#> Beta for b_x 1.11 1.017
#>
#> Expected number of effective parameters(stdev): 504.01(0.001)
#> Number of equivalent replicates : 2.98
#>
#> Marginal log-Likelihood: -6852.87
# On the other hand
summary(inla(
formula = Y ~ -1 +
f(b_x, copy = "idx_x", hyper = list(beta = list(param = c(0, 0.1), fixed = FALSE))) +
f(idx_x, weight_x, model = "iid", values = 1:500, hyper = list(
prec = list(initial = 15, fixed = TRUE)
)) +
b_0 + b_z + a_0 + a_z,
family = c("poisson", "gaussian", "gaussian"),
data = my_data,
verbose = TRUE
))
#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.
#> If this does not help, please contact the developers at <[email protected]>.
Created
March 8, 2021 15:25
-
-
Save agila5/90b13440dff7c2cc727c97ca04476a19 to your computer and use it in GitHub Desktop.
INLA-ME
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment