Estimate the four parameters of the (bivariate) \(N_1\)--Poisson copula mixed data type model by maximum likelihood estimation.
N1poisson(lmean = "identitylink", lsd = "loglink",
lvar = "loglink", llambda = "loglink", lapar = "rhobitlink",
zero = c(if (var.arg) "var" else "sd", "apar"),
doff = 5, nnodes = 20, copula = "gaussian",
var.arg = FALSE, imethod = 1, isd = NULL,
ilambda = NULL, iapar = NULL)
An object of class "vglmff"
(see vglmff-class
).
The object is used by modelling functions
such as vglm
and vgam
.
Details at CommonVGAMffArguments
.
See Links
for more link function choices.
The second response is primarily controlled by
the parameter \(\lambda_2\).
Initial values.
Details at CommonVGAMffArguments
.
Details at CommonVGAMffArguments
.
Numeric of unit length, the denominator offset
\(\delta>0\).
A monotonic transformation
\(\Delta^* = \lambda_2^{2/3} /
(|\delta| + \lambda_2^{2/3})\)
is taken to map the Poisson mean onto the
unit interval.
This argument is \(\delta\).
The default reflects the property that the normal
approximation to the Poisson work wells for
\(\lambda_2 \geq 10\) or thereabouts, hence
that value is mapped to the origin by
qnorm
.
That's because 10**(2/3)
is approximately 5.
It is known that the \(\lambda_2\) rate
parameter raised to
the power of \(2/3\) is a transformation
that approximates the normal density more
closely.
Alternatively,
delta
may be assigned a single
negative value. If so, then
\(\Delta^* = \log(1 + \lambda_2)
/ [|\delta| + \log(1 + \lambda_2)]\)
is used.
For this, doff = -log1p(10)
is
suggested.
Details at N1binomial
.
See uninormal
.
T. W. Yee
The bivariate response comprises
\(Y_1\) from a linear model
having parameters
mean
and sd
for
\(\mu_1\) and \(\sigma_1\),
and the Poisson count
\(Y_2\) having parameter
lambda
for its mean \(\lambda_2\).
The
joint probability density/mass function is
\(P(y_1, Y_2 = y_2) = \phi_1(y_1; \mu_1, \sigma_1)
\exp(-h^{-1}(\Delta))
[h^{-1}(\Delta)]^{y_2} / y_2!\)
where \(\Delta\) adjusts \(\lambda_2\)
according to the association parameter
\(\alpha\).
The quantity \(\Delta\) is
\(\Phi((\Phi^{-1}(h(\lambda_2)) -
\alpha Z_1) / \sqrt{1 - \alpha^2})\)
where \(h\) maps
\(\lambda_2\) onto the unit interval.
The quantity \(Z_1\) is \((Y_1-\mu_1) / \sigma_1\).
Thus there is an underlying bivariate normal
distribution, and a copula is used to bring the
two marginal distributions together.
Here,
\(-1 < \alpha < 1\), and
\(\Phi\) is the
cumulative distribution function
pnorm
of a standard univariate normal.
The first marginal
distribution is a normal distribution
for the linear model.
The second column of the response must
have nonnegative integer values.
When \(\alpha = 0\)
then \(\Delta=\Delta^*\).
Together, this family function combines
uninormal
and
poissonff
.
If the response are correlated then
a more efficient joint analysis
should result.
The second marginal distribution allows for overdispersion relative to an ordinary Poisson distribution---a property due to \(\alpha\).
This VGAM family function cannot handle multiple responses. Only a two-column matrix is allowed. The two-column fitted value matrix has columns \(\mu_1\) and \(\lambda_2\).
rN1pois
,
N1binomial
,
binormalcop
,
uninormal
,
poissonff
,
dpois
.
apar <- rhobitlink(0.3, inverse = TRUE)
nn <- 1000; mymu <- 1; sdev <- exp(1)
lambda <- loglink(1, inverse = TRUE)
mat <- rN1pois(nn, mymu, sdev, lambda, apar)
npdata <- data.frame(y1 = mat[, 1], y2 = mat[, 2])
with(npdata, var(y2) / mean(y2)) # Overdispersion
fit1 <- vglm(cbind(y1, y2) ~ 1, N1poisson,
npdata, trace = TRUE)
coef(fit1, matrix = TRUE)
Coef(fit1)
head(fitted(fit1))
summary(fit1)
confint(fit1)
Run the code above in your browser using DataLab