VGAMextra (version 0.0-6)

MVNcov: Multivariate Normal Distribution Family Function

Description

Maximum likelihood estimation of the Multivariate Normal distribution. The covariances (not correlation coefficients) are included in the parameter vector.

Usage

MVNcov(zero = c("var", "cov"),
             lmean = "identitylink",
             lvar  = "loglink",
             lcov  = "identitylink")

Value

An object of class "vglmff"

(see vglmff-class) to be used by VGLM/VGAM modelling functions, e.g., vglm or vgam.

Arguments

zero

Integer or character--string vector. Which linear predictors are intercept--only. Details at zero or CommonVGAMffArguments.

lmean, lvar, lcov

VGLM--link functions applied to the means, variances and covariances.

Author

Victor Miranda.

Details

For the \(K\)--dimensional normal distribution, this fits a linear model to the \(K\) means \(\mu_j\) \(j = 1, \ldots, K\), which are the first entries in the linear predictor. The variances \(\sigma^2_j\) \(j = 1, \ldots, K\) and then the covariances \(cov_{ij}\) \(i = 1, \ldots, K, j = i + 1, \ldots, K\), are next aligned. The fitted means are returned as the fitted values.

The log--likelihood is computed via dmultinorm, an implementation of the multivariate Normal density.

The score and expected information matrices are internally computed at each Fisher scoring step, using its vectorized form.

The response should be an \(K\)--column matrix. The covariances may be any real number so that the identitylink is a reasonable choice. For further details on VGLM/VGAM--link functions, see Links.

See Also

dmultinorm, binormal, CommonVGAMffArguments, vglm.

Examples

# K = 3.
set.seed(180227)
nn  <- 85
mvndata <- data.frame(x2 = runif(nn), x3 = runif(nn))
mvndata <- transform(mvndata, 
                     y = rbinorm(nn, mean1 = 2 - 2 * x2 + 1 * x3,
                          mean2 = 2 - 1.5 * x3,
                          var1 = exp(1.0), var2 = exp(-0.75),
                          cov12 = 0.5 * exp(1.0) * exp(-0.75)))
mvndata <- transform(mvndata, y3 = rnorm(nn, mean = 2 + x2, sd = exp(1.5)))
colnames(mvndata) <- c("x2", "x3", "y1", "y2", "y3")

mvnfit <- vglm(cbind(y1, y2, y3) ~ x2 + x3, MVNcov, data = mvndata, trace = TRUE)
(mvncoef  <- coef(mvnfit, mat = TRUE))

## Check variances and covariances: var1, var2 and var3.
exp(mvncoef[c(10, 13, 16)])  # True are var1 = exp(1.0) = 2.718, 
                             # var2 = exp(-0.75) = 0.472
                             # and var3 = exp(1.5)^2 = 20.08554
vcov(mvnfit)
constraints(mvnfit)
summary(mvnfit)