MVNcov {VGAMextra}

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")
 

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.

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.

Value

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

Note

Unlike other implementations, e.g., binormal from VGAM in terms of \rho and standard deviations, MVNcov estimates the variances and covariances, modeled as intercept–only. See argument zero, whose default is c("var", "cov"), to change this.

Thus far, there is no guarantee that the estimated var–cov matrix will be positive–definite. Proper procedures to validate this will be incorporated shortly, such as the @validparams slot.

Although the function has been tested on K \leq 5 data sets, it is recommended that K \leq 3, unless the data are nice and n is sufficiently large.

Author(s)

Victor Miranda.

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)