trinormal: Trivariate Normal Distribution Family Function

Description Usage Arguments Details Value Warning Author(s) See Also Examples

View source: R/family.bivariate.R

Description

Maximum likelihood estimation of the nine parameters of a trivariate normal distribution.

Usage

1
2
3
4
5
6
7
8
9
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 Links for 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 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.

imean1, imean2, imean3, isd1, isd2, isd3

See CommonVGAMffArguments for more information.

irho12, irho23, irho13, imethod, zero

See CommonVGAMffArguments for 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.

Author(s)

T. W. Yee

See Also

uninormal, binormal, rtrinorm.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
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)
coef(fit1, matrix = TRUE)
constraints(fit1)
summary(fit1)
## Not run:  # 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)
coef(fit2, matrix = TRUE)

## End(Not run)

Example output

Loading required package: stats4
Loading required package: splines
VGLM    linear loop  1 :  loglikelihood = -4546.3465
VGLM    linear loop  2 :  loglikelihood = -4297.6827
VGLM    linear loop  3 :  loglikelihood = -4235.841
VGLM    linear loop  4 :  loglikelihood = -4234.035
VGLM    linear loop  5 :  loglikelihood = -4234.0344
VGLM    linear loop  6 :  loglikelihood = -4234.0344
               mean1    mean2    mean3 loglink(sd1) loglink(sd2) loglink(sd3)
(Intercept) 1.021178 2.929393 3.931580 -0.007345242 -0.007345242 -0.007345242
x2          2.042807 4.101541 5.119168  0.000000000  0.000000000  0.000000000
            rhobitlink(rho12) rhobitlink(rho23) rhobitlink(rho13)
(Intercept)        0.04491489        0.04491489        0.04491489
x2                 0.00000000        0.00000000        0.00000000
$`(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


Call:
vglm(formula = cbind(y1, y2, y3) ~ x2, family = trinormal(eq.sd = TRUE, 
    eq.cor = TRUE), data = tdata, trace = TRUE)

Pearson residuals:
                      Min      1Q    Median     3Q   Max
mean1             -3.0864 -0.6981 -0.003713 0.7012 3.352
mean2             -2.8441 -0.6473 -0.027948 0.6673 3.430
mean3             -3.0890 -0.6359  0.013576 0.6672 2.895
loglink(sd1)      -0.7073 -0.6264 -0.358355 0.2583 7.252
loglink(sd2)      -0.7075 -0.6349 -0.403910 0.1667 7.632
loglink(sd3)      -0.7080 -0.6353 -0.409792 0.2380 6.056
rhobitlink(rho12) -4.7085 -0.3709  0.003533 0.3462 4.814
rhobitlink(rho23) -5.1662 -0.3395  0.003960 0.3526 4.658
rhobitlink(rho13) -6.9591 -0.3840 -0.003700 0.3454 4.722

Coefficients: 
               Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  1.021178   0.062744  16.275   <2e-16 ***
(Intercept):2  2.929393   0.062744  46.688   <2e-16 ***
(Intercept):3  3.931580   0.062744  62.661   <2e-16 ***
(Intercept):4 -0.007345   0.012916  -0.569    0.570    
(Intercept):5  0.044915   0.037317   1.204    0.229    
x2:1           2.042807   0.109248  18.699   <2e-16 ***
x2:2           4.101541   0.109248  37.543   <2e-16 ***
x2:3           5.119168   0.109248  46.858   <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: -4234.034 on 8992 degrees of freedom

Number of Fisher scoring iterations: 6 

Taking a modified step.
VGLM    linear loop  2 :  loglikelihood = -4487.9134
Taking a modified step.
VGLM    linear loop  3 :  loglikelihood = -4420.4091
Taking a modified step.
VGLM    linear loop  4 :  loglikelihood = -4363.3008
Taking a modified step.
VGLM    linear loop  5 :  loglikelihood = -4317.4177
Taking a modified step.
VGLM    linear loop  6 :  loglikelihood = -4283.4266
Taking a modified step.
VGLM    linear loop  7 :  loglikelihood = -4261.5633
Taking a modified step.
VGLM    linear loop  8 :  loglikelihood = -4250.4528
Taking a modified step.
VGLM    linear loop  9 :  loglikelihood = -4246.1332
Taking a modified step.
VGLM    linear loop  10 :  loglikelihood = -4239.8236
Taking a modified step.
VGLM    linear loop  11 :  loglikelihood = -4238.1585
Taking a modified step.
VGLM    linear loop  12 :  loglikelihood = -4236.4815
Taking a modified step.....
VGLM    linear loop  13 :  loglikelihood = -4235.9344
Taking a modified step.
VGLM    linear loop  14 :  loglikelihood = -4235.3495
Taking a modified step...
VGLM    linear loop  15 :  loglikelihood = -4235.2012
Taking a modified step.
VGLM    linear loop  16 :  loglikelihood = -4235.0892
Taking a modified step.....
VGLM    linear loop  17 :  loglikelihood = -4234.6041
Taking a modified step.
VGLM    linear loop  18 :  loglikelihood = -4234.5771
Taking a modified step....
VGLM    linear loop  19 :  loglikelihood = -4234.3132
Taking a modified step.
VGLM    linear loop  20 :  loglikelihood = -4234.2062
Taking a modified step..
VGLM    linear loop  21 :  loglikelihood = -4234.1684
Taking a modified step..
VGLM    linear loop  22 :  loglikelihood = -4234.1362
Taking a modified step..
VGLM    linear loop  23 :  loglikelihood = -4234.1131
Taking a modified step..
VGLM    linear loop  24 :  loglikelihood = -4234.095
Taking a modified step..
VGLM    linear loop  25 :  loglikelihood = -4234.0811
Taking a modified step..
VGLM    linear loop  26 :  loglikelihood = -4234.0707
Taking a modified step..
VGLM    linear loop  27 :  loglikelihood = -4234.0623
Taking a modified step..
VGLM    linear loop  28 :  loglikelihood = -4234.0562
Taking a modified step..
VGLM    linear loop  29 :  loglikelihood = -4234.0512
Taking a modified step..
VGLM    linear loop  30 :  loglikelihood = -4234.0476
Warning message:
In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = Xm2,  :
  convergence not obtained in 30 IRLS iterations
               mean1    mean2    mean3 loglink(sd1) loglink(sd2) loglink(sd3)
(Intercept) 1.021178 2.929393 3.931580 -0.005420121 -0.005420121 -0.005420121
x2          2.042807 4.101541 5.119168  0.000000000  0.000000000  0.000000000
            extlogitlink(rho12, min = -0.5, max = 1)
(Intercept)                               -0.6298517
x2                                         0.0000000
            extlogitlink(rho23, min = -0.5, max = 1)
(Intercept)                               -0.6298517
x2                                         0.0000000
            extlogitlink(rho13, min = -0.5, max = 1)
(Intercept)                               -0.6298517
x2                                         0.0000000

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.