Multivariate multiplicative binomial distribution

Description

Various utilities to coerce and manipulate MB objects

Usage

MB(dep, m, pnames=character(0))
## S3 method for class 'MB'
as.array(x, ...)
## S4 method for signature 'MB'
getM(x)
## S3 method for class 'gunter_MB'
print(x, ...)

Arguments

Details

Function MB() returns an object of class MB. This is essentially a matrix with one row corresponding to a single observation; repeated rows indicate identical observations as shown below. Observational data is typically in this form. The idea is that the user can coerce to a gunter_MB object, which is then analyzable by Lindsey().

The multivariate multiplicative binomial distribution is defined by

Equation 20 of the vignette

Thus if θ=φ=1 the system reduces to a product of independent binomial distributions with probability p_i and size m_i for 1,...,t.

There follows a short R transcript showing the MB class in use, with annotation.

The first step is to define an m vector:

R> m <- c(2,3,1)
 

This means that m1=2,m2=3,m3=1. So m1=2 means that i=1 corresponds to a binomial distribution with size 2 [that is, the observation is in the set {0,1,2}]; and m2=3 means that i=2 corresponds to a binomial with size 3 [ie the set {0,1,2,3}].

Now we need some observations:

R> a <- matrix(c(1,0,0, 1,0,0, 1,1,1, 2,3,1, 2,0,1),5,3,byrow=T)
R> a
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    1    0    0
[3,]    1    1    1
[4,]    2    3    1
[5,]    2    0    1 

In matrix a, the first observation, viz c(1,0,0) is interpreted as x1=1,x2=0,x3=0. Thus, because x_i+z_i=m_i, we have z1=1,z2=3,z3=1. Now we can create an object of class MB, using function MB():

R>  mx <- MB(a, m, letters[1:3])   

The third argument gives names to the observations corresponding to the columns of a. The values of m1,m2,m3 may be extracted using getM():

R> getM(mx)
a b c 
2 3 1 
R> 

The getM() function returns a named vector, with names given as the third argument to MB().

Now we illustrate the print method:

R> mx
     a na     b nb     c nc    
[1,] 1 1      0 3      0 1     
[2,] 1 1      0 3      0 1     
[3,] 1 1      1 2      1 0     
[4,] 2 0      3 0      1 0     
[5,] 2 0      0 3      1 0     
R> 

See how the columns are in pairs: the first pair total 2 (because m1=2), the second pair total 3 (because m2=3), and the third pair total 1 (because m3=1). Each pair of columns has only a single degree of freedom, because m_i is known.

Also observe how the column names are in pairs. The print method puts these in place. Take the first two columns. These are named ‘a’ and ‘na’: this is intented to mean ‘a’ and ‘not a’.

We can now coerce to a gunter_MB:

R> (gx <- gunter(mx))
$tbl
   a b c
1  0 0 0
2  1 0 0
3  2 0 0
[snip]
24 2 3 1

$d
 [1] 0 2 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1

$m
a b c 
2 3 1 

Take the second line of the element tbl of gx, as an example. This reads c(1,0,0) corresponding to the observations of a,b,c respectively, and the second line of element d [“d” for “data”], viz 2, shows that this observation occurred twice (and in fact these were the first two lines of a).

Now we can coerce object mx to an array:

R> (ax <- as.array(mx))
, , c = 0

   b
a   0 1 2 3
  0 0 0 0 0
  1 0 0 2 0
  2 0 0 0 0

, , c = 1

   b
a   0 1 2 3
  0 0 1 0 0
  1 0 0 0 0
  2 1 1 0 0
>

(actually, ax is an Oarray object). The location of an element in ax corresponds to an observation of abc, and the entry corresponds to the number of times that observation was made. For example, ax[1,2,0]=2 shows that c(1,2,0) occurred twice (the first two lines of a).

The Lindsey Poisson device is applicable: see help(danaher) for an application to the bivariate case and help(Lindsey) for an example where a table is created from scratch.

Author(s)

Robin K. S. Hankin

See Also

MM, Lindsey, danaher

Examples

a <- matrix(c(1,0,0, 1,0,0, 1,1,1, 2,3,1, 2,0,1),5,3,byrow=TRUE)
m <- c(2,3,1)
mx <- MB(a, m, letters[1:3])   # mx is of class 'MB'; column headings
                   #  mean "a" and "not a".
ax <- as.array(mx)
gx <- gunter(ax)
ax2 <- as.array(gx)

data(danaher)
summary(Lindsey_MB(danaher))

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(MM)
Loading required package: magic
Loading required package: abind
Loading required package: partitions
Loading required package: emulator
Loading required package: mvtnorm
Loading required package: Oarray
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MM/MB.Rd_%03d_medium.png", width=480, height=480)
> ### Name: MB
> ### Title: Multivariate multiplicative binomial distribution
> ### Aliases: MB MB-class as.array.MB as.array.gunter_MB print.gunter_MB
> ###   counts,MB-method counts getM,MB-method getM
> 
> ### ** Examples
> 
> 
> a <- matrix(c(1,0,0, 1,0,0, 1,1,1, 2,3,1, 2,0,1),5,3,byrow=TRUE)
> m <- c(2,3,1)
> mx <- MB(a, m, letters[1:3])   # mx is of class 'MB'; column headings
>                    #  mean "a" and "not a".
> ax <- as.array(mx)
> gx <- gunter(ax)
> ax2 <- as.array(gx)
> 
> data(danaher)
> summary(Lindsey_MB(danaher))

Call:
glm(formula = d ~ (.), family = poisson, data = x, offset = Off)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.87863  -0.57513   0.03813   0.50323   1.58913  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)     5.51168    0.06062  90.916  < 2e-16 ***
bacon          -1.65283    0.16666  -9.918  < 2e-16 ***
eggs           -1.04872    0.08102 -12.944  < 2e-16 ***
`bacon:nbacon` -0.51486    0.06266  -8.217  < 2e-16 ***
`eggs:neggs`   -0.35546    0.04058  -8.759  < 2e-16 ***
`bacon:eggs`    0.30061    0.05980   5.027 4.98e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3046.440  on 24  degrees of freedom
Residual deviance:   18.666  on 19  degrees of freedom
AIC: 108.26

Number of Fisher Scoring iterations: 5

> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>