Multivariate Normal Model
mvnorm.RmdThe family mvnorm_bamlss implements the multivariate normal distribution that can be used for estimation with bamlss. To illustrate the usage we simulate data from a trivariate normal.
library("bamlss")
library("mvtnorm")
## Set the seed for reproducability.
set.seed(123)
## Simulate data.
n <- 1000
## True functions.
f0 <- function(x) 2 * sin(pi * x)
f1 <- function(x) exp(2 * x)
f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 * (1 - x)^10
f3 <- function(x) sin(x * 3) - 3
f4 <- function(x) cos(x * 6) - 2
f5 <- function(x) {
eta <- sin(scale2(x, -3, 3)) - 3
eta / sqrt(1 + eta^2)
}
f6 <- function(x) {
eta <- cos(scale2(x, -3, 3)) - 2
eta / sqrt(1 + eta^2)
}
## Covariate data.
d <- data.frame(
"x0" = runif(n),
"x1" = runif(n),
"x2" = runif(n),
"x3" = runif(n)
)
s1 <- exp(f3(d$x1))
s2 <- exp(f4(d$x2))
s3 <- exp(rep(-1, n))
rho12 <- f5(d$x3)
rho13 <- rep(0.1, n)
rho23 <- f6(d$x3)
y <- matrix(0, n, 3)
## Generate response observations.
for(i in 1:n) {
V <- diag(c(s1[i], s2[i], s3[i])^2)
V[1, 2] <- rho12[i] * s1[i] * s2[i]
V[1, 3] <- rho13[i] * s1[i] * s3[i]
V[2, 3] <- rho23[i] * s2[i] * s3[i]
V[2, 1] <- V[1, 2]
V[3, 1] <- V[1, 3]
V[3, 2] <- V[2, 3]
mu <- c(f0(d$x0[i]) + f1(d$x1[i]), f2(d$x2[i]), f3(d$x3[i]))
y[i, ] <- rmvn(1, mu, V)
}
## Specify response names and add to data frame.
colnames(y) <- paste0("y", 1:3)
d <- cbind(d, as.data.frame(y))Now we can setup the model formula with
f <- list(
y1 ~ s(x0) + s(x1),
y2 ~ s(x2),
y3 ~ s(x3),
sigma1 ~ s(x1),
sigma2 ~ s(x2),
sigma3 ~ 1,
rho12 ~ s(x3),
rho13 ~ 1,
rho23 ~ s(x3)
)Note that the numbering of the correlation parameters is somewhat untrivial, but it corresponds to the position of the correlation matrix. The model can be estimated using the usual specifications, but the dimension of the multivariate normal distribution must be specified directly in the family.
set.seed(123)
b <- bamlss(f, data = d, family = mvnorm_bamlss(k = 3))After the MCMC algorithm is completed, all summary and extraction functions can be used as usual.
summary(b)##
## Call:
## bamlss(formula = f, family = mvnorm_bamlss(k = 3), data = d)
## ---
## Family: mvnorm
## Link function: mu1 = identity, mu2 = identity, mu3 = identity, sigma1 = log, sigma2 = log, sigma3 = log, rho12 = rhogit, rho13 = rhogit, rho23 = rhogit
## *---
## Formula mu1:
## ---
## y1 ~ s(x0) + s(x1)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) 4.463 4.456 4.463 4.471 4.464
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 1 1 1 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x0).tau21 1.4957 0.4989 1.2566 3.8044 0.031
## s(x0).alpha 1.0000 1.0000 1.0000 1.0000 NA
## s(x0).edf 8.8549 8.6823 8.8675 8.9566 6.663
## s(x1).tau21 1.5608 0.5454 1.2539 4.7034 0.195
## s(x1).alpha 1.0000 1.0000 1.0000 1.0000 NA
## s(x1).edf 8.8365 8.6602 8.8473 8.9575 8.206
## ---
## Formula mu2:
## ---
## y2 ~ s(x2)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) 3.276 3.260 3.276 3.295 3.272
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 1 1 1 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x2).tau21 581.594 196.200 472.467 1570.624 138.045
## s(x2).alpha 1.000 1.000 1.000 1.000 NA
## s(x2).edf 8.999 8.999 8.999 9.000 8.998
## ---
## Formula mu3:
## ---
## y3 ~ s(x3)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) -2.325 -2.347 -2.325 -2.304 -2.32
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 1 1 1 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x3).tau21 0.7537 0.2247 0.5956 2.3505 0.104
## s(x3).alpha 1.0000 1.0000 1.0000 1.0000 NA
## s(x3).edf 7.0921 6.0282 7.0858 8.2334 5.160
## ---
## Formula sigma1:
## ---
## sigma1 ~ s(x1)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) -2.086 -2.128 -2.087 -2.044 -2.085
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 0.9757 0.8250 0.9999 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x1).tau21 1.3533 0.1459 0.8480 5.0489 2.804
## s(x1).alpha 0.8541 0.3320 0.9356 1.0000 NA
## s(x1).edf 5.6462 3.7602 5.6178 7.5910 6.999
## ---
## Formula sigma2:
## ---
## sigma2 ~ s(x2)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) -1.777 -1.815 -1.776 -1.742 -1.774
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 0.9865 0.8944 1.0000 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x2).tau21 73.8518 22.2223 55.4773 219.8731 85.182
## s(x2).alpha 0.8510 0.3495 0.9390 1.0000 NA
## s(x2).edf 8.8633 8.6955 8.8754 8.9686 8.916
## ---
## Formula sigma3:
## ---
## sigma3 ~ 1
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) -1.0301 -1.0693 -1.0298 -0.9938 -1.022
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 0.9621 0.7437 0.9991 1
## ---
## Formula rho12:
## ---
## rho12 ~ s(x3)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) -0.8469 -0.8874 -0.8501 -0.7702 -0.866
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 6.278e-01 1.070e-61 8.759e-01 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x3).tau21 1.284e+01 3.322e+00 9.790e+00 4.524e+01 24.314
## s(x3).alpha 2.554e-01 4.625e-67 7.836e-02 1.000e+00 NA
## s(x3).edf 7.995e+00 6.998e+00 8.001e+00 8.726e+00 8.569
## ---
## Formula rho13:
## ---
## rho13 ~ 1
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) 0.03327 -0.02900 0.03953 0.05855 0.037
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 6.586e-01 5.487e-60 9.194e-01 1
## ---
## Formula rho23:
## ---
## rho23 ~ s(x3)
## -
## Parametric coefficients:
## Mean 2.5% 50% 97.5% parameters
## (Intercept) -0.8673 -0.9094 -0.8686 -0.8193 -0.869
## -
## Acceptance probability:
## Mean 2.5% 50% 97.5%
## alpha 6.547e-01 2.169e-64 9.121e-01 1
## -
## Smooth terms:
## Mean 2.5% 50% 97.5% parameters
## s(x3).tau21 1.067e+00 8.240e-02 5.090e-01 7.421e+00 0.035
## s(x3).alpha 2.871e-01 4.360e-186 7.989e-02 1.000e+00 NA
## s(x3).edf 5.174e+00 3.581e+00 5.058e+00 7.802e+00 3.008
## ---
## Sampler summary:
## -
## DIC = -2986.723 logLik = 1526.518 pd = 66.3138
## runtime = 147.739
## ---
## Optimizer summary:
## -
## AICc = -2971.752 edf = 65.5196 logLik = 1556.064
## logPost = 1193.908 nobs = 1000 runtime = 110.476
plot(b, ask = FALSE)
Umlauf, Nikolaus, Nadja Klein, Achim Zeileis, and Thorsten Simon. 2021. bamlss: Bayesian Additive Models for Location Scale and Shape (and Beyond). https://CRAN.R-project.org/package=bamlss.