Trivariate Normal Distribution Family Function
trinormal.RdMaximum likelihood estimation of the nine parameters of a trivariate normal distribution.
Usage
trinormal(zero = c("sd", "rho"), eq.mean = FALSE,
eq.sd = FALSE, eq.cor = FALSE,
lmean1 = "identitylink", lmean2 = "identitylink",
lmean3 = "identitylink",
lsd1 = "loglink", lsd2 = "loglink", lsd3 = "loglink",
lrho12 = "rhobitlink", lrho23 = "rhobitlink", lrho13 = "rhobitlink",
imean1 = NULL, imean2 = NULL, imean3 = NULL,
isd1 = NULL, isd2 = NULL, isd3 = NULL,
irho12 = NULL, irho23 = NULL, irho13 = NULL, imethod = 1)Arguments
- lmean1, lmean2, lmean3, lsd1, lsd2, lsd3
Link functions applied to the means and standard deviations. See
Linksfor more choices. Being positive quantities, a log link is the default for the standard deviations.- lrho12, lrho23, lrho13
Link functions applied to the correlation parameters. See
Linksfor more choices. By default the correlation parameters are allowed to have a value between -1 and 1, but that may be problematic wheneq.cor = TRUEbecause they should have a value between -0.5 and 1.- imean1, imean2, imean3, isd1, isd2, isd3
See
CommonVGAMffArgumentsfor more information.- irho12, irho23, irho13, imethod, zero
See
CommonVGAMffArgumentsfor more information.- eq.mean, eq.sd, eq.cor
Logical. Constrain the means or the standard deviations or correlation parameters to be equal?
Details
For the trivariate normal distribution,
this fits a linear model (LM) to the means, and
by default,
the other parameters are intercept-only.
The response should be a three-column matrix.
The three correlation parameters are prefixed by rho,
and the default gives them values between \(-1\) and \(1\)
however, this may be problematic when the correlation parameters
are constrained to be equal, etc..
The fitted means are returned as the fitted values, which is in
the form of a three-column matrix.
Fisher scoring is implemented.
Value
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm,
and vgam.
Warning
The default parameterization does not make the estimated
variance-covariance matrix positive-definite.
In order for the variance-covariance matrix to be positive-definite
the quantity
1 - rho12^2 - rho13^2 - rho23^2 + 2 * rho12 * rho13 * rho23
must be positive, and if eq.cor = TRUE then
this means that the rhos must be between -0.5 and 1.
Examples
if (FALSE) set.seed(123); nn <- 1000
tdata <- data.frame(x2 = runif(nn), x3 = runif(nn))
tdata <- transform(tdata, y1 = rnorm(nn, 1 + 2 * x2),
y2 = rnorm(nn, 3 + 4 * x2),
y3 = rnorm(nn, 4 + 5 * x2))
fit1 <- vglm(cbind(y1, y2, y3) ~ x2, data = tdata,
trinormal(eq.sd = TRUE, eq.cor = TRUE), trace = TRUE)
#> Iteration 1: loglikelihood = -4634.3328
#> Iteration 2: loglikelihood = -4373.5692
#> Iteration 3: loglikelihood = -4305.3938
#> Iteration 4: loglikelihood = -4303.299
#> Iteration 5: loglikelihood = -4303.2983
#> Iteration 6: loglikelihood = -4303.2983
coef(fit1, matrix = TRUE)
#> mean1 mean2 mean3 loglink(sd1) loglink(sd2) loglink(sd3)
#> (Intercept) 1.014695 2.970145 3.980051 0.01582878 0.01582878 0.01582878
#> x2 2.055752 4.056072 4.990266 0.00000000 0.00000000 0.00000000
#> rhobitlink(rho12) rhobitlink(rho23) rhobitlink(rho13)
#> (Intercept) -0.05128255 -0.05128255 -0.05128255
#> x2 0.00000000 0.00000000 0.00000000
constraints(fit1)
#> $`(Intercept)`
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 1 0 0 0 0
#> [2,] 0 1 0 0 0
#> [3,] 0 0 1 0 0
#> [4,] 0 0 0 1 0
#> [5,] 0 0 0 1 0
#> [6,] 0 0 0 1 0
#> [7,] 0 0 0 0 1
#> [8,] 0 0 0 0 1
#> [9,] 0 0 0 0 1
#>
#> $x2
#> [,1] [,2] [,3]
#> [1,] 1 0 0
#> [2,] 0 1 0
#> [3,] 0 0 1
#> [4,] 0 0 0
#> [5,] 0 0 0
#> [6,] 0 0 0
#> [7,] 0 0 0
#> [8,] 0 0 0
#> [9,] 0 0 0
#>
summary(fit1)
#>
#> Call:
#> vglm(formula = cbind(y1, y2, y3) ~ x2, family = trinormal(eq.sd = TRUE,
#> eq.cor = TRUE), data = tdata, trace = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 1.01470 0.06342 15.998 <2e-16 ***
#> (Intercept):2 2.97014 0.06342 46.830 <2e-16 ***
#> (Intercept):3 3.98005 0.06342 62.753 <2e-16 ***
#> (Intercept):4 0.01583 0.01292 1.225 0.220
#> (Intercept):5 -0.05128 0.03555 -1.442 0.149
#> x2:1 2.05575 0.11051 18.603 <2e-16 ***
#> x2:2 4.05607 0.11051 36.704 <2e-16 ***
#> x2:3 4.99027 0.11051 45.157 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Number of linear predictors: 9
#>
#> Names of linear predictors: mean1, mean2, mean3, loglink(sd1), loglink(sd2),
#> loglink(sd3), rhobitlink(rho12), rhobitlink(rho23), rhobitlink(rho13)
#>
#> Log-likelihood: -4303.298 on 8992 degrees of freedom
#>
#> Number of Fisher scoring iterations: 6
#>
# Try this when eq.sd = TRUE, eq.cor = TRUE:
fit2 <-
vglm(cbind(y1, y2, y3) ~ x2, data = tdata, stepsize = 0.25,
trinormal(eq.sd = TRUE, eq.cor = TRUE,
lrho12 = "extlogitlink(min = -0.5)",
lrho23 = "extlogitlink(min = -0.5)",
lrho13 = "extlogitlink(min = -0.5)"), trace = TRUE)
#> Taking a modified step.
#> Iteration 2 : loglikelihood = -4572.0611
#> Taking a modified step.
#> Iteration 3 : loglikelihood = -4500.0551
#> Taking a modified step.
#> Iteration 4 : loglikelihood = -4438.0946
#> Taking a modified step.
#> Iteration 5 : loglikelihood = -4387.3534
#> Taking a modified step.
#> Iteration 6 : loglikelihood = -4349.4727
#> Taking a modified step.
#> Iteration 7 : loglikelihood = -4329.3141
#> Taking a modified step.
#> Iteration 8 : loglikelihood = -4324.0094
#> Taking a modified step...................
coef(fit2, matrix = TRUE)
#> mean1 mean2 mean3 loglink(sd1) loglink(sd2) loglink(sd3)
#> (Intercept) 1.014695 2.970145 3.980051 0.08877413 0.08877413 0.08877413
#> x2 2.055752 4.056072 4.990266 0.00000000 0.00000000 0.00000000
#> extlogitlink(rho12, min = -0.5, max = 1)
#> (Intercept) -0.9916732
#> x2 0.0000000
#> extlogitlink(rho23, min = -0.5, max = 1)
#> (Intercept) -0.9916732
#> x2 0.0000000
#> extlogitlink(rho13, min = -0.5, max = 1)
#> (Intercept) -0.9916732
#> x2 0.0000000
# \dontrun{}