Bivariate Normal Distribution Family Function

Maximum likelihood estimation of the five parameters of a bivariate normal distribution.

Usage

binormal(lmean1 = "identitylink", lmean2 = "identitylink",
         lsd1   = "loglink",     lsd2   = "loglink",
         lrho   = "rhobitlink",
         imean1 = NULL,       imean2 = NULL,
         isd1   = NULL,       isd2   = NULL,
         irho   = NULL,       imethod = 1,
         eq.mean = FALSE,     eq.sd = FALSE,
         zero = c("sd", "rho"), rho.arg = NA)

Arguments

lmean1, lmean2, lsd1, lsd2, lrho

Link functions applied to the means, standard deviations and rho parameters. See Links for more choices. Being positive quantities, a log link is the default for the standard deviations.

imean1, imean2, isd1, isd2, irho, imethod, zero

See CommonVGAMffArguments for more information.

eq.mean, eq.sd

Logical or formula. Constrains the means or the standard deviations to be equal.

rho.arg

If \(\rho\) is known then this argument may be assigned the (scalar) value lying in \((-1, 1)\). The default is to estimate that parameter so that \(M=5\). If known, then other arguments such as lrho and irho are ignored, and "rho" is removed from zero.

Details

For the bivariate 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 two-column matrix. The correlation parameter is rho, which lies between \(-1\) and \(1\) (thus the rhobitlink link is a reasonable choice). The fitted means are returned as the fitted values, which is in the form of a two-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

This function may be renamed to normal2() or something like that at a later date.

Author

T. W. Yee

Note

If both equal means and equal standard deviations are desired then use something like constraints = list("(Intercept)" = matrix(c(1,1,0,0,0, 0,0,1,1,0 ,0,0,0,0,1), 5, 3)) and maybe zero = NULL etc.

Examples

set.seed(123); nn <- 1000
bdata <- data.frame(x2 = runif(nn), x3 = runif(nn))
bdata <- transform(bdata, y1 = rnorm(nn, 1 + 2 * x2),
                          y2 = rnorm(nn, 3 + 4 * x2))
fit1 <- vglm(cbind(y1, y2) ~ x2,
             binormal(eq.sd = TRUE), bdata, trace = TRUE)
#> Iteration 1: loglikelihood = -2982.0414
#> Iteration 2: loglikelihood = -2842.6316
#> Iteration 3: loglikelihood = -2824.4788
#> Iteration 4: loglikelihood = -2824.2708
#> Iteration 5: loglikelihood = -2824.2708
#> Iteration 6: loglikelihood = -2824.2708
coef(fit1, matrix = TRUE)
#>                mean1    mean2 loglink(sd1) loglink(sd2) rhobitlink(rho)
#> (Intercept) 1.021178 2.929393 -0.006632272 -0.006632272      0.05229185
#> x2          2.042807 4.101541  0.000000000  0.000000000      0.00000000
constraints(fit1)
#> $`(Intercept)`
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    0    0
#> [2,]    0    1    0    0
#> [3,]    0    0    1    0
#> [4,]    0    0    1    0
#> [5,]    0    0    0    1
#> 
#> $x2
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    0    1
#> [3,]    0    0
#> [4,]    0    0
#> [5,]    0    0
#> 
summary(fit1)
#> 
#> Call:
#> vglm(formula = cbind(y1, y2) ~ x2, family = binormal(eq.sd = TRUE), 
#>     data = bdata, trace = TRUE)
#> 
#> Coefficients: 
#>                Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1  1.021178   0.062789  16.264   <2e-16 ***
#> (Intercept):2  2.929393   0.062789  46.655   <2e-16 ***
#> (Intercept):3 -0.006632   0.015817  -0.419    0.675    
#> (Intercept):4  0.052292   0.063246   0.827    0.408    
#> x2:1           2.042807   0.109326  18.685   <2e-16 ***
#> x2:2           4.101541   0.109326  37.517   <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Number of linear predictors:  5 
#> 
#> Names of linear predictors: mean1, mean2, loglink(sd1), loglink(sd2), 
#> rhobitlink(rho)
#> 
#> Log-likelihood: -2824.271 on 4994 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 6 
#> 

# Estimated P(Y1 <= y1, Y2 <= y2) under the fitted model
var1  <- loglink(2 * predict(fit1)[, "loglink(sd1)"], inv = TRUE)
var2  <- loglink(2 * predict(fit1)[, "loglink(sd2)"], inv = TRUE)
cov12 <- rhobitlink(predict(fit1)[, "rhobitlink(rho)"], inv = TRUE)
head(with(bdata, pbinorm(y1, y2,
                         mean1 = predict(fit1)[, "mean1"],
                         mean2 = predict(fit1)[, "mean2"],
                         var1 = var1, var2 = var2, cov12 = cov12)))
#> [1] 0.049938306 0.082072022 0.148273155 0.377605648 0.002533645 0.247922216