oiposbinomial | R Documentation |
Fits a one-inflated positive binomial distribution by maximum likelihood estimation.
oiposbinomial(lpstr1 = "logitlink", lprob = "logitlink",
type.fitted = c("mean", "prob", "pobs1", "pstr1", "onempstr1"),
iprob = NULL, gpstr1 = ppoints(9), gprob = ppoints(9),
multiple.responses = FALSE, zero = NULL)
lpstr1 , lprob |
Link functions for the parameter |
type.fitted |
See |
iprob , gpstr1 , gprob |
For initial values;
see |
multiple.responses |
Logical.
See |
zero |
See |
These functions are based on
P(Y=y) = \phi + (1-\phi) N \mu (1-\mu)^N / (1-(1-\mu)^N),
for y=1/N
, and
P(Y=y) = (1-\phi) {N \choose Ny} \mu^{Ny}
(1-\mu)^{N(1-y)} / (1-(1-\mu)^N).
for y=2/N,\ldots,1
. That is, the response is a sample
proportion out of N
trials, and the argument size
in
roiposbinom
is N
here.
Ideally N > 2
is needed.
The parameter \phi
is the probability of a structural one,
and it satisfies 0 < \phi < 1
(usually).
The mean of Y
is
E(Y)=\phi + (1-\phi) \mu / (1-(1-\mu)^N)
and these are returned as the default fitted values.
By default, the two linear/additive predictors
for oiposbinomial()
are (logit(\phi), logit(\mu))^T
.
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions
such as vglm
and vgam
.
The response variable should have one of the formats described by
binomialff
, e.g., a factor or two column matrix
or a vector of sample proportions with the weights
argument specifying the values of N
.
To work well, one ideally needs large values of N
and
\mu
much greater than 0, i.e., the larger N
and
\mu
are, the better. If N = 1
then the model
is unidentifiable since the number of parameters is excessive.
Estimated probabilities of a structural one and an
observed one are returned, as in zipoisson
.
The one-deflated positive binomial distribution might
be fitted by setting lpstr1 = "identitylink"
, albeit,
not entirely reliably. See zipoisson
for information that can be applied here.
T. W. Yee
roiposbinom
,
posbinomial
,
binomialff
,
rbinom
.
size <- 10 # Number of trials; N in the notation above
nn <- 200
odata <- data.frame(pstr1 = logitlink( 0, inv = TRUE), # 0.50
mubin1 = logitlink(-1, inv = TRUE), # Binomial mean
svec = rep(size, length = nn),
x2 = runif(nn))
odata <- transform(odata,
mubin2 = logitlink(-1 + x2, inv = TRUE))
odata <- transform(odata,
y1 = roiposbinom(nn, svec, pr = mubin1, pstr1 = pstr1),
y2 = roiposbinom(nn, svec, pr = mubin2, pstr1 = pstr1))
with(odata, table(y1))
fit1 <- vglm(y1 / svec ~ 1, oiposbinomial, data = odata,
weights = svec, trace = TRUE, crit = "coef")
fit2 <- vglm(y2 / svec ~ x2, oiposbinomial, data = odata,
weights = svec, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1) # Useful for intercept-only models
head(fitted(fit1, type = "pobs1")) # Estimate of P(Y = 1)
head(fitted(fit1))
with(odata, mean(y1)) # Compare this with fitted(fit1)
summary(fit1)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.