trinormal.RdMaximum likelihood estimation of the nine parameters of a trivariate normal distribution.
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)Link functions applied to the means and standard deviations.
See Links for more choices.
Being positive quantities, a log link is the default for the
standard deviations.
Link functions applied to the correlation parameters.
See Links for more choices.
By default the correlation parameters are allowed to have
a value between -1 and 1, but that may be problematic
when eq.cor = TRUE because they should have a value
between -0.5 and 1.
See CommonVGAMffArguments for more information.
See CommonVGAMffArguments for more information.
Logical. Constrain the means or the standard deviations or correlation parameters to be equal?
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.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm,
and vgam.
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.
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{}