binormal: Bivariate Normal Distribution Family Function

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

View source: R/family.bivariate.R

Description

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

Usage

1
2
3
4
5
6
7
8
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"))

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.

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.

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.

Author(s)

T. W. Yee

See Also

uninormal, trinormal, pbinorm, bistudentt.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
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), data = bdata, trace = TRUE)
coef(fit1, matrix = TRUE)
constraints(fit1)
summary(fit1)

# Estimated P(Y1 <= y1, Y2 <= y2) under the fitted model
var1  <- loglink(2 * predict(fit1)[, "loglink(sd1)"], inverse = TRUE)
var2  <- loglink(2 * predict(fit1)[, "loglink(sd2)"], inverse = TRUE)
cov12 <- rhobitlink(predict(fit1)[, "rhobitlink(rho)"], inverse = TRUE)
head(with(bdata, pbinorm(y1, y2,
                         mean1 = predict(fit1)[, "mean1"],
                         mean2 = predict(fit1)[, "mean2"],
                         var1 = var1, var2 = var2, cov12 = cov12)))

Example output

Loading required package: stats4
Loading required package: splines
VGLM    linear loop  1 :  loglikelihood = -3012.1567
VGLM    linear loop  2 :  loglikelihood = -2851.7724
VGLM    linear loop  3 :  loglikelihood = -2824.3618
VGLM    linear loop  4 :  loglikelihood = -2823.7626
VGLM    linear loop  5 :  loglikelihood = -2823.7624
VGLM    linear loop  6 :  loglikelihood = -2823.7624
               mean1    mean2 loglink(sd1) loglink(sd2) rhobitlink(rho)
(Intercept) 1.021178 2.929393  0.009053151  -0.02282574      0.05231844
x2          2.042807 4.101541  0.000000000   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    0    1

$x2
     [,1] [,2]
[1,]    1    0
[2,]    0    1
[3,]    0    0
[4,]    0    0
[5,]    0    0


Call:
vglm(formula = cbind(y1, y2) ~ x2, family = binormal(eq.sd = TRUE), 
    data = bdata, trace = TRUE)

Pearson residuals:
                    Min      1Q    Median     3Q   Max
mean1           -3.0223 -0.6832  0.004237 0.6924 3.299
mean2           -2.8869 -0.6539 -0.030958 0.6704 3.481
loglink(sd1)    -0.7071 -0.6281 -0.375149 0.2344 7.000
loglink(sd2)    -0.7072 -0.6334 -0.396139 0.1919 7.886
rhobitlink(rho) -4.7126 -0.3756  0.002839 0.3409 4.735

Coefficients: 
               Estimate Std. Error z value Pr(>|z|)    
(Intercept):1  1.021178   0.063781  16.011   <2e-16 ***
(Intercept):2  2.929393   0.061780  47.416   <2e-16 ***
(Intercept):3  0.009053   0.022361   0.405    0.686    
(Intercept):4 -0.022826   0.022361  -1.021    0.307    
(Intercept):5  0.052318   0.063246   0.827    0.408    
x2:1           2.042807   0.111054  18.395   <2e-16 ***
x2:2           4.101541   0.107570  38.129   <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: -2823.762 on 4993 degrees of freedom

Number of Fisher scoring iterations: 6 

No Hauck-Donner effect found in any of the estimates

[1] 0.050752063 0.084538836 0.146966302 0.380193747 0.002858084 0.244247322

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