MB | R Documentation |
Various utilities to coerce and manipulate MB
objects
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, ...)
dep |
Primary argument to |
m |
Vector containing the relative sizes of the various marginal binomial distributions |
x |
Object of class |
... |
Further arguments to |
pnames |
In function |
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 \mjdeqn \prod_i=1^t m_i\choose x_i\, z_ip_i^x_iq_i^z_i\theta_i^x_iz_i \prod_i < j\phi_ij^x_ix_j Equation 20 of the vignette
Thus if \mjeqn\theta=\phi=1theta=phi=1 the system reduces to a product of independent binomial distributions with probability \mjseqnp_i and size \mjseqnm_i for \mjeqni=1,...,t1,...,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 \mjeqnm_1=2,m_2=3,m_3=1m1=2,m2=3,m3=1. So \mjeqnm_1=2m1=2 means that \mjseqni=1 corresponds to a binomial distribution with size 2 [that is, the observation is in the set \mjeqn{0,1,2}0,1,2]; and \mjeqnm_2=3m2=3 means that \mjseqni=2 corresponds to a binomial with size 3 [ie the set \mjeqn{0,1,2,3}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 \mjeqnx_1=1,x_2=0,x_3=0x1=1,x2=0,x3=0. Thus, because
\mjeqnx_i+z_i=m_ixi+zi=mi, we have \mjeqnz_1=1,z_2=3,z_3=1z1=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 \mjeqnm_1, m_2, m_3m1,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 \mjeqnm_1=2m1=2), the second pair total 3 (because \mjeqnm_2=3m2=3), and the third pair total 1 (because \mjeqnm_3=1m3=1). Each pair of columns has only a single degree of freedom, because \mjeqnm_imi 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.
Robin K. S. Hankin
MM
, Lindsey
, danaher
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))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.