Nothing
# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.
dzanegbin <- function(x, size, # prob = NULL,
munb, # = NULL,
pobs0 = 0,
log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(pobs0), length(size))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
if (!is.Numeric(size, positive = TRUE))
stop("argument 'size' must be in (0,Inf)")
index0 <- x == 0
if (log.arg) {
ans[ index0] <- log(pobs0[index0])
if (any(!index0))
ans[!index0] <- log1p(-pobs0[!index0]) +
dgaitdnbinom(x[!index0],
size[!index0], truncate = 0,
munb.p = munb[!index0],
log = TRUE)
} else {
ans[ index0] <- pobs0[index0]
if (any(!index0))
ans[!index0] <- (1 - pobs0[!index0]) *
dgaitdnbinom(x[!index0], size[!index0],
munb.p = munb[!index0],
truncate = 0)
}
ans
} # dzanegbin
pzanegbin <- function(q, size, # prob = NULL,
munb, # = NULL,
pobs0 = 0) {
LLL <- max(length(q), length(pobs0), length(size))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
qindex <- (q > 0)
ans[ qindex] <- pobs0[qindex] + (1 - pobs0[qindex]) *
pgaitdnbinom(q[qindex], size[qindex],
munb.p = munb[qindex],
truncate = 0)
ans[q < 0] <- 0
ans[q == 0] <- pobs0[q == 0]
ans <- pmax(0, ans)
ans <- pmin(1, ans)
ans
} # pzanegbin
qzanegbin <- function(p, size, # prob = NULL,
munb, # = NULL,
pobs0 = 0) {
LLL <- max(length(p), length(pobs0), length(size))
if (length(p) != LLL) p <- rep_len(p, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be between 0 and 1 inclusive")
ans <- p
ans[p <= pobs0] <- 0
pindex <- (p > pobs0)
ans[pindex] <-
qgaitdnbinom((p[pindex] - pobs0[pindex]) / (1 -
pobs0[pindex]),
size[pindex],
munb.p = munb[pindex],
truncate = 0)
ans
} # qzanegbin
rzanegbin <- function(n, size, # prob = NULL,
munb, # = NULL,
pobs0 = 0) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
ans <- rgaitdnbinom(n = use.n, size, # prob,
munb.p = munb,
truncate = 0)
if (length(pobs0) != use.n)
pobs0 <- rep_len(pobs0, use.n)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be between 0 and 1 inclusive")
ifelse(runif(use.n) < pobs0, 0, ans)
} # rzanegbin
dzapois <- function(x, lambda, pobs0 = 0, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(lambda), length(pobs0))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
index0 <- (x == 0)
if (log.arg) {
ans[ index0] <- log(pobs0[index0])
ans[!index0] <- log1p(-pobs0[!index0]) +
dgaitdpois(x[!index0], lambda[!index0],
log = TRUE, truncate = 0)
} else {
ans[ index0] <- pobs0[index0]
ans[!index0] <- (1 - pobs0[!index0]) *
dgaitdpois(x[!index0], lambda[!index0],
truncate = 0)
}
ans
}
pzapois <- function(q, lambda, pobs0 = 0) {
LLL <- max(length(q), length(lambda), length(pobs0))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
ans[q > 0] <- pobs0[q > 0] +
(1 - pobs0[q > 0]) * pgaitdpois(q[q > 0], lambda[q > 0],
truncate = 0)
ans[q < 0] <- 0
ans[q == 0] <- pobs0[q == 0]
ans <- pmax(0, ans)
ans <- pmin(1, ans)
ans
}
qzapois <- function(p, lambda, pobs0 = 0) {
LLL <- max(length(p), length(lambda), length(pobs0))
if (length(p) != LLL) p <- rep_len(p, LLL)
if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be between 0 and 1 inclusive")
ans <- p
ind4 <- (p > pobs0)
ans[!ind4] <- 0
ans[ ind4] <- qgaitdpois((p[ind4] - pobs0[ind4]) / (
1 - pobs0[ind4]),
lambda[ind4], truncate = 0)
ans
}
rzapois <- function(n, lambda, pobs0 = 0) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
ans <- rgaitdpois(use.n, lambda, truncate = 0)
if (length(pobs0) != use.n)
pobs0 <- rep_len(pobs0, use.n)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must in [0,1]")
ifelse(runif(use.n) < pobs0, 0, ans)
}
dzipois <- function(x, lambda, pstr0 = 0, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(lambda), length(pstr0))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- x + lambda + pstr0
index0 <- (x == 0)
if (log.arg) {
ans[ index0] <- log(pstr0[ index0] + (1 - pstr0[ index0]) *
dpois(x[ index0], lambda[ index0]))
ans[!index0] <- log1p(-pstr0[!index0]) +
dpois(x[!index0], lambda[!index0],
log = TRUE)
} else {
ans[ index0] <- pstr0[ index0] + (1 - pstr0[ index0]) *
dpois(x[ index0], lambda[ index0])
ans[!index0] <- (1 - pstr0[!index0]) *
dpois(x[!index0], lambda[!index0])
}
deflat.limit <- -1 / expm1(lambda)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
pzipois <- function(q, lambda, pstr0 = 0) {
LLL <- max(length(pstr0), length(lambda), length(q))
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)
if (length(q) != LLL) q <- rep_len(q, LLL)
ans <- ppois(q, lambda)
ans <- ifelse(q < 0, 0, pstr0 + (1 - pstr0) * ans)
deflat.limit <- -1 / expm1(lambda)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
qzipois <- function(p, lambda, pstr0 = 0) {
LLL <- max(length(p), length(lambda), length(pstr0))
if (length(p) != LLL) p <- rep_len(p, LLL)
if (length(lambda) != LLL) lambda <- rep_len(lambda, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- rep_len(NA_real_, LLL)
deflat.limit <- -1 / expm1(lambda)
ans[p <= pstr0] <- 0
pindex <- (pstr0 < p) & (deflat.limit <= pstr0)
ans[pindex] <-
qpois((p[pindex] - pstr0[pindex]) / (1 - pstr0[pindex]),
lambda = lambda[pindex])
ans[pstr0 < deflat.limit] <- NaN
ans[1 < pstr0] <- NaN
ans[lambda < 0] <- NaN
ans[p < 0] <- NaN
ans[1 < p] <- NaN
ans
} # qzipois
rzipois <- function(n, lambda, pstr0 = 0) {
qzipois(runif(n), lambda, pstr0 = pstr0)
}
zapoisson <-
function(lpobs0 = "logitlink", llambda = "loglink",
type.fitted = c("mean", "lambda", "pobs0", "onempobs0"),
imethod = 1,
ipobs0 = NULL, ilambda = NULL, ishrinkage = 0.95,
probs.y = 0.35,
zero = NULL) {
lpobs.0 <- as.list(substitute(lpobs0))
epobs.0 <- link2list(lpobs.0)
lpobs.0 <- attr(epobs.0, "function.name")
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0",
"onempobs0"))[1]
new("vglmff",
blurb = c("Zero-altered Poisson ",
"(Bernoulli and positive-Poisson conditional ",
"model)\n\n",
"Links: ",
namesof("pobs0", lpobs.0, epobs.0, tag = FALSE), ", ",
namesof("lambda", llambda, elambda, tag = FALSE), "\n",
"Mean: (1 - pobs0) * lambda / (1 - exp(-lambda))"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("pobs0", "lambda"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero, .type.fitted = type.fitted ))),
rqresslot = eval(substitute(
function(mu, y, w, eta, extra = NULL) {
TFvec <- c(TRUE, FALSE)
phimat <- eta2theta(eta[, TFvec, drop = FALSE],
.lpobs.0 , earg = .epobs.0 )
lambda <- eta2theta(eta[, !TFvec, drop = FALSE],
.llambda , earg = .elambda )
scrambleseed <- runif(1) # To scramble the seed
qnorm(runif(length(y),
pzapois(y - 1, pstr0 = phimat, lambda = lambda),
pzapois(y , pstr0 = phimat, lambda = lambda)))
},
list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$skip.these <- skip.these <-
matrix(as.logical(y0), n, NOS)
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("pobs0", ncoly, skip1 = TRUE)
mynames2 <- param.names("lambda", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lpobs.0 , .epobs.0 , tag = FALSE),
namesof(mynames2, .llambda , .elambda , tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
if (!length(etastart)) {
lambda.init <- Init.mu(y = y, w = w, imethod = .imethod ,
imu = .ilambda , # x = x,
ishrinkage = .ishrinkage ,
pos.only = TRUE,
probs.y = .probs.y )
etastart <-
cbind(theta2eta(if (length( .ipobs0 )) .ipobs0 else
(0.5 + w * y0) / (1 + w),
.lpobs.0 , .epobs.0 ),
theta2eta(lambda.init, .llambda , .elambda ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
}
}), list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda,
.ipobs0 = ipobs0, .ilambda = ilambda,
.ishrinkage = ishrinkage, .probs.y = probs.y,
.imethod = imethod,
.type.fitted = type.fitted ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0", "onempobs0"))[1]
M1 <- 2
NOS <- ncol(eta) / M1
pobs.0 <- cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lpobs.0 , earg = .epobs.0 ))
lambda <- cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.llambda , earg = .elambda ))
ans <- switch(type.fitted,
"mean" = (1 - pobs.0) * lambda / (-expm1(-lambda)),
"lambda" = lambda,
"pobs0" = pobs.0, # P(Y=0)
"onempobs0" = 1 - pobs.0) # P(Y>0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
last = eval(substitute(expression({
temp.names <- c(rep_len( .lpobs.0 , NOS),
rep_len( .llambda , NOS))
temp.names <- temp.names[interleave.VGAM(M1*NOS, M1 = M1)]
misc$link <- temp.names
names(misc$link) <-
c(mynames1, mynames2)[interleave.VGAM(M1*NOS, M1 = M1)]
misc$earg <- vector("list", M1 * NOS)
names(misc$earg) <- names(misc$link)
for (ii in 1:NOS) {
misc$earg[[M1*ii-1]] <- .epobs.0
misc$earg[[M1*ii ]] <- .elambda
}
}), list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
pobs0 <- cbind(eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
.lpobs.0, earg = .epobs.0))
lambda <- cbind(eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
.llambda, earg = .elambda ))
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dzapois(x = y, pobs0 = pobs0,
lambda = lambda, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
vfamily = c("zapoisson"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
TFvec <- c(TRUE, FALSE)
phimat <- eta2theta(eta[, TFvec, drop = FALSE],
.lpobs.0 , earg = .epobs.0 )
lambda <- eta2theta(eta[, !TFvec, drop = FALSE],
.llambda , earg = .elambda )
okay1 <- all(is.finite(lambda)) && all(0 < lambda) &&
all(is.finite(phimat)) &&
all(0 < phimat & phimat < 1)
okay1
}, list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
pobs0 <- eta2theta(eta[, c(TRUE, FALSE)],
.lpobs.0 , earg = .epobs.0 )
lambda <- eta2theta(eta[, c(FALSE, TRUE)],
.llambda , earg = .elambda )
rzapois(nsim * length(lambda), lambda, pobs0 = pobs0)
}, list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
deriv = eval(substitute(expression({
M1 <- 2
NOS <- ncol(eta) / M1 # extra$NOS
y0 <- extra$y0
skip <- extra$skip.these
TFvec <- c(TRUE, FALSE)
phimat <- eta2theta(eta[, TFvec, drop = FALSE],
.lpobs.0 , earg = .epobs.0 )
lambda <- eta2theta(eta[, !TFvec, drop = FALSE],
.llambda , earg = .elambda )
dl.dlambda <- y / lambda + 1 / expm1(-lambda)
dl.dphimat <- -1 / (1 - phimat) # For y > 0 obsns
for (spp. in 1:NOS) {
dl.dphimat[skip[, spp.], spp.] <-
1 / phimat[skip[, spp.], spp.]
dl.dlambda[skip[, spp.], spp.] <- 0
}
dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda )
mu.phi0 <- phimat
temp3 <- if ( .lpobs.0 == "logitlink") {
c(w) * (y0 - mu.phi0)
} else {
c(w) * dtheta.deta(mu.phi0, .lpobs.0 , .epobs.0 ) *
dl.dphimat
}
ans <- cbind(temp3,
c(w) * dl.dlambda * dlambda.deta)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
ans
}), list( .lpobs.0 = lpobs.0, .llambda = llambda,
.epobs.0 = epobs.0, .elambda = elambda ))),
weight = eval(substitute(expression({
wz <- matrix(0.0, n, M1 * NOS)
temp5 <- expm1(lambda)
ned2l.dlambda2 <- (1 - phimat) * (temp5 + 1) *
(1 / lambda - 1 / temp5) / temp5
wz[, NOS+(1:NOS)] <- c(w) * ned2l.dlambda2 * dlambda.deta^2
tmp100 <- mu.phi0 * (1 - mu.phi0)
tmp200 <- if ( .lpobs.0 == "logitlink" &&
is.empty.list( .epobs.0 )) {
cbind(c(w) * tmp100)
} else {
cbind(c(w) * (1 / tmp100) *
dtheta.deta(mu.phi0, .lpobs.0 , .epobs.0 )^2)
}
if (FALSE)
for (ii in 1:NOS) {
index200 <- abs(tmp200[, ii]) < .Machine$double.eps
if (any(index200)) {
tmp200[index200, ii] <- 10.0 * .Machine$double.eps^(3/4)
}
}
wz[, 1:NOS] <- tmp200
wz <- wz[, interleave.VGAM(ncol(wz), M1 = M1)]
wz
}), list( .lpobs.0 = lpobs.0,
.epobs.0 = epobs.0 ))))
} # zapoisson
zapoissonff <-
function(llambda = "loglink", lonempobs0 = "logitlink",
type.fitted = c("mean", "lambda",
"pobs0", "onempobs0"),
imethod = 1,
ilambda = NULL, ionempobs0 = NULL,
ishrinkage = 0.95,
probs.y = 0.35,
zero = "onempobs0") {
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lonempobs0 <- as.list(substitute(lonempobs0))
eonempobs0 <- link2list(lonempobs0)
lonempobs0 <- attr(eonempobs0, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0", "onempobs0"))[1]
new("vglmff",
blurb = c("Zero-altered Poisson ",
"(Bernoulli and positive-Poisson ",
"conditional model)\n\n",
"Links: ",
namesof("lambda", llambda, earg = elambda,
tag = FALSE), ", ",
namesof("onempobs0", lonempobs0, earg = eonempobs0,
tag = FALSE), "\n",
"Mean: onempobs0 * lambda / (1 - exp(-lambda))"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("lambda", "onempobs0"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted
))),
rqresslot = eval(substitute(
function(mu, y, w, eta, extra = NULL) {
NOS <- extra$NOS
M1 <- 2
lambda <- cbind(eta2theta(eta[, M1*(1:NOS)-1,
drop = FALSE],
.llambda , .elambda ))
onempobs0 <- cbind(eta2theta(eta[, M1*(1:NOS)-0,
drop = FALSE],
.lonempobs0 , .eonempobs0 ))
scrambleseed <- runif(1) # To scramble the seed
qnorm(runif(length(y),
pzapois(y - 1, pobs0 = 1 - onempobs0, lambda),
pzapois(y , pobs0 = 1 - onempobs0, lambda)))
}, list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.integer.y = TRUE,
Is.nonnegative.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$skip.these <- skip.these <- matrix(as.logical(y0),
n, NOS)
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("lambda", ncoly, skip1 = TRUE)
mynames2 <- param.names("onempobs0", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .llambda, .elambda , tag = FALSE),
namesof(mynames2, .lonempobs0 , .eonempobs0 , tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
if (!length(etastart)) {
lambda.init <- Init.mu(y = y, w = w, imethod = .imethod ,
imu = .ilambda, # x = x,
ishrinkage = .ishrinkage,
pos.only = TRUE,
probs.y = .probs.y )
etastart <-
cbind(theta2eta(lambda.init, .llambda , .elambda ),
theta2eta(1 - (0.5 + w * y0) / (1 + w),
.lonempobs0 , earg = .eonempobs0 ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
}
}), list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda,
.ilambda = ilambda,
.ishrinkage = ishrinkage, .probs.y = probs.y,
.type.fitted = type.fitted,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0", "onempobs0"))[1]
M1 <- 2
NOS <- ncol(eta) / M1
lambda <-
cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.llambda , .elambda ))
onempobs0 <-
cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lonempobs0 , .eonempobs0 ))
ans <- switch(type.fitted,
"mean" = onempobs0 * lambda / (-expm1(-lambda)),
"lambda" = lambda,
"pobs0" = 1 - onempobs0, # P(Y=0)
"onempobs0" = onempobs0) # P(Y>0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
last = eval(substitute(expression({
misc$expected <- TRUE
misc$multipleResponses <- TRUE
temp.names <- c(rep_len( .llambda , NOS),
rep_len( .lonempobs0 , NOS))
temp.names <- temp.names[interleave.VGAM(M1*NOS, M1 = M1)]
misc$link <- temp.names
names(misc$link) <-
c(mynames1, mynames2)[interleave.VGAM(M1*NOS, M1 = M1)]
misc$earg <- vector("list", M1 * NOS)
names(misc$earg) <- names(misc$link)
for (ii in 1:NOS) {
misc$earg[[M1*ii-1]] <- .elambda
misc$earg[[M1*ii ]] <- .eonempobs0
}
}), list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
NOS <- extra$NOS
M1 <- 2
lambda <- cbind(eta2theta(eta[, M1*(1:NOS)-1,
drop = FALSE],
.llambda , .elambda ))
onempobs0 <- cbind(eta2theta(eta[, M1*(1:NOS)-0,
drop = FALSE],
.lonempobs0 , .eonempobs0 ))
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzapois(y, lambda, pobs0 = 1 - onempobs0,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
vfamily = c("zapoissonff"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
TFvec <- c(TRUE, FALSE)
lambda <- eta2theta(eta[, TFvec, drop=FALSE],
.llambda , e= .elambda )
onempobs0 <- eta2theta(eta[, !TFvec, drop=FALSE],
.lonempobs0 , e= .eonempobs0 )
okay1 <- all(is.finite(lambda)) && all(0 < lambda) &&
all(is.finite(onempobs0)) &&
all(0 < onempobs0 & onempobs0 < 1)
okay1
}, list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempobs0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempobs0 ,
earg = .eonempobs0 )
rzapois(nsim * length(lambda), lambda = lambda,
pobs0 = 1 - onempobs0)
}, list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
deriv = eval(substitute(expression({
M1 <- 2
NOS <- ncol(eta) / M1 # extra$NOS
y0 <- extra$y0
skip <- extra$skip.these
lambda <- cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.llambda, earg = .elambda ))
omphimat <- cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lonempobs0, earg = .eonempobs0 ))
phimat <- 1 - omphimat
dl.dlambda <- y / lambda + 1 / expm1(-lambda)
dl.dPHImat <- +1 / (omphimat) # For y > 0 obsns
for (spp. in 1:NOS) {
dl.dPHImat[skip[, spp.], spp.] <-
-1 / phimat[skip[, spp.], spp.]
dl.dlambda[skip[, spp.], spp.] <- 0
}
dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda )
mu.phi0 <- omphimat
temp3 <- if ( FALSE && .lonempobs0 == "logitlink") {
} else {
c(w) * dtheta.deta(mu.phi0, .lonempobs0 , earg = .eonempobs0 ) *
dl.dPHImat
}
ans <- cbind(c(w) * dl.dlambda * dlambda.deta,
temp3)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
ans
}), list( .lonempobs0 = lonempobs0, .llambda = llambda,
.eonempobs0 = eonempobs0, .elambda = elambda ))),
weight = eval(substitute(expression({
wz <- matrix(0.0, n, M1 * NOS)
temp5 <- expm1(lambda)
ned2l.dlambda2 <- (1 - phimat) * (temp5 + 1) *
(1 / lambda - 1 / temp5) / temp5
wz[, 0 * NOS + (1:NOS)] <- c(w) * ned2l.dlambda2 *
dlambda.deta^2
tmp100 <- mu.phi0 * (1.0 - mu.phi0)
tmp200 <- if ( .lonempobs0 == "logitlink" &&
is.empty.list( .eonempobs0 )) {
cbind(c(w) * tmp100)
} else {
cbind(c(w) * (1 / tmp100) *
dtheta.deta(mu.phi0, link = .lonempobs0,
earg = .eonempobs0)^2)
}
wz[, 1 * NOS + (1:NOS)] <- tmp200
wz <- wz[, interleave.VGAM(ncol(wz), M1 = M1)]
wz
}), list( .lonempobs0 = lonempobs0,
.eonempobs0 = eonempobs0 ))))
} # End of zapoissonff
zanegbinomial.control <-
function(save.weights = TRUE,
summary.HDEtest = FALSE, # Overwrites summary() default.
...) {
list(save.weights = save.weights,
summary.HDEtest = summary.HDEtest)
}
zanegbinomial <-
function(
zero = "size",
type.fitted = c("mean", "munb", "pobs0"),
mds.min = 1e-3,
nsimEIM = 500,
cutoff.prob = 0.999, # higher is better for large 'size'
eps.trig = 1e-7,
max.support = 4000, # 20160127; I have changed this
max.chunk.MB = 30, # max.memory = Inf is allowed
lpobs0 = "logitlink", lmunb = "loglink", lsize = "loglink",
imethod = 1,
ipobs0 = NULL,
imunb = NULL,
iprobs.y = NULL,
gprobs.y = (0:9)/10, # 20160709; grid for
isize = NULL,
gsize.mux = exp(c(-30, -20, -15, -10, -6:3))) {
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(eps.trig, length.arg = 1,
positive = TRUE) || eps.trig > 0.001)
stop("argument 'eps.trig' must be positive and",
" smaller in value")
if (!is.Numeric(nsimEIM, length.arg = 1,
positive = TRUE, integer.valued = TRUE))
stop("argument 'nsimEIM' must be a positive integer")
if (nsimEIM <= 30)
warning("argument 'nsimEIM' should be greater than 30, say")
if (length(ipobs0) &&
(!is.Numeric(ipobs0, positive = TRUE) ||
max(ipobs0) >= 1))
stop("If given, argument 'ipobs0' must contain",
"values in (0,1) only")
if (length(isize) && !is.Numeric(isize, positive = TRUE))
stop("If given, argument 'isize' must contain ",
"positive values only")
lpobs0 <- as.list(substitute(lpobs0))
epobs0 <- link2list(lpobs0)
lpobs0 <- attr(epobs0, "function.name")
lmunb <- as.list(substitute(lmunb))
emunb <- link2list(lmunb)
lmunb <- attr(emunb, "function.name")
lsize <- as.list(substitute(lsize))
esize <- link2list(lsize)
lsize <- attr(esize, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0"))[1]
ipobs0.small <- 1/64 # A number easily represented exactly
new("vglmff",
blurb = c("Zero-altered negative binomial (Bernoulli and\n",
"positive-negative binomial conditional model)\n\n",
"Links: ",
namesof("pobs0", lpobs0, epobs0, tag = FALSE), ", ",
namesof("munb", lmunb, emunb, tag = FALSE), ", ",
namesof("size", lsize, esize, tag = FALSE), "\n",
"Mean: (1 - pobs0) * munb / (1 - (size / (size + ",
"munb))^size)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 3,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
mds.min = .mds.min ,
imethod = .imethod ,
multipleResponses = TRUE,
parameters.names = c("pobs0", "munb", "size"),
nsimEIM = .nsimEIM ,
eps.trig = .eps.trig ,
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero, .imethod = imethod,
.nsimEIM = nsimEIM, .eps.trig = eps.trig,
.type.fitted = type.fitted,
.mds.min = mds.min
))),
initialize = eval(substitute(expression({
M1 <- 3
temp16 <-
w.y.check(w = w, y = y,
Is.integer.y = TRUE,
Is.nonnegative.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp16$w
y <- temp16$y
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
M <- M1 * ncoly
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("pobs0", NOS, skip1 = TRUE)
mynames2 <- param.names("munb", NOS, skip1 = TRUE)
mynames3 <- param.names("size", NOS, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lpobs0 , .epobs0 , tag = FALSE),
namesof(mynames2, .lmunb , .emunb , tag = FALSE),
namesof(mynames3, .lsize , .esize , tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$skip.these <-
skip.these <- matrix(as.logical(y0), n, NOS)
gprobs.y <- .gprobs.y
imunb <- .imunb # Default in NULL
if (length(imunb))
imunb <- matrix(imunb, n, NOS, byrow = TRUE)
if (!length(etastart)) {
munb.init <-
size.init <- matrix(NA_real_, n, NOS)
gprobs.y <- .gprobs.y
if (length( .iprobs.y ))
gprobs.y <- .iprobs.y
gsize.mux <- .gsize.mux # gsize.mux is on a relative scale
for (jay in 1:NOS) { # For each response 'y_jay'... do:
TFvec <- y[, jay] > 0 # Important to exclude the 0s
posyvec <- y[TFvec, jay]
munb.init.jay <- if ( .imethod == 1 ) {
quantile(posyvec, probs = gprobs.y) - 1/2 # + 1/16
} else {
weighted.mean(posyvec, w = w[TFvec, jay]) - 1/2
}
if (length(imunb))
munb.init.jay <- imunb[, jay]
gsize <- gsize.mux * 0.5 * (mean(munb.init.jay) +
weighted.mean(posyvec, w = w[TFvec, jay]))
if (length( .isize ))
gsize <- .isize # isize is on an absolute scale
try.this <-
grid.search2(munb.init.jay, gsize,
objfun = posNBD.Loglikfun2,
y = posyvec, # x = x[TFvec, , drop = FALSE],
w = w[TFvec, jay],
ret.objfun = TRUE) # Last value is the loglik
munb.init[, jay] <- try.this["Value1"]
size.init[, jay] <- try.this["Value2"]
} # for (jay ...)
pobs0.init <- matrix(if (length( .ipobs0 ))
.ipobs0 else -1,
nrow = n, ncol = NOS, byrow = TRUE)
for (jay in 1:NOS) {
if (any(pobs0.init[, jay] < 0)) {
index.y0 <- (y[, jay] < 0.5)
pobs0.init[, jay] <-
max(min(weighted.mean(index.y0, w[, jay]),
1 - .ipobs0.small ),
.ipobs0.small )
}
}
etastart <- cbind(theta2eta(pobs0.init, .lpobs0 , .epobs0 ),
theta2eta(munb.init, .lmunb , .emunb ),
theta2eta(size.init, .lsize , .esize ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
} # End of if (!length(etastart))
}),
list( .lpobs0 = lpobs0, .lmunb = lmunb, .lsize = lsize,
.epobs0 = epobs0, .emunb = emunb, .esize = esize,
.ipobs0 = ipobs0, .isize = isize,
.ipobs0.small = ipobs0.small,
.imunb = imunb,
.gprobs.y = gprobs.y, .gsize.mux = gsize.mux,
.imethod = imethod,
.type.fitted = type.fitted, .iprobs.y = iprobs.y ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0"))[1]
M1 <- 3
NOS <- ncol(eta) / M1
phi0 <- eta2theta(eta[, M1*(1:NOS)-2], .lpobs0 , .epobs0 )
munb <- eta2theta(eta[, M1*(1:NOS)-1], .lmunb , .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS) ], .lsize , .esize )
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
prob0 <- tempk^kmat # p(0) from negative binomial
oneminusf0 <- 1 - prob0
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
prob0[big.size] <- exp(-munb[big.size]) # The limit
oneminusf0[big.size] <- -expm1(-munb[big.size])
}
ans <- switch(type.fitted,
"mean" = (1 - phi0) * munb / oneminusf0,
"munb" = munb,
"pobs0" = phi0) # P(Y=0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lpobs0 = lpobs0, .lsize = lsize, .lmunb = lmunb,
.epobs0 = epobs0, .emunb = emunb, .esize = esize,
.mds.min = mds.min ))),
last = eval(substitute(expression({
misc$link <-
c(rep_len( .lpobs0 , NOS),
rep_len( .lmunb , NOS),
rep_len( .lsize , NOS))[interleave.VGAM(M1 * NOS,
M1 = M1)]
temp.names <- c(mynames1,
mynames2,
mynames3)[interleave.VGAM(M1*NOS, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M1*NOS)
names(misc$earg) <- temp.names
for (ii in 1:NOS) {
misc$earg[[M1*ii - 2]] <- .epobs0
misc$earg[[M1*ii - 1]] <- .emunb
misc$earg[[M1*ii ]] <- .esize
}
misc$nsimEIM <- .nsimEIM
misc$ipobs0 <- .ipobs0
misc$isize <- .isize
misc$multipleResponses <- TRUE
}), list( .lpobs0 = lpobs0, .lmunb = lmunb, .lsize = lsize,
.epobs0 = epobs0, .emunb = emunb, .esize = esize,
.ipobs0 = ipobs0, .isize = isize,
.nsimEIM = nsimEIM ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
M1 <- 3
NOS <- ncol(eta) / M1
phi0 <- eta2theta(eta[, M1*(1:NOS)-2], .lpobs0 , .epobs0 )
munb <- eta2theta(eta[, M1*(1:NOS)-1], .lmunb , .emunb )
size <- eta2theta(eta[, M1*(1:NOS) ], .lsize , .esize )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzanegbin(x = y, pobs0 = phi0,
munb = munb, size = size,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
},
list( .lpobs0 = lpobs0, .lmunb = lmunb, .lsize = lsize,
.epobs0 = epobs0, .emunb = emunb, .esize = esize ))),
vfamily = c("zanegbinomial"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
phi0 <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
.lpobs0 , earg = .epobs0 )
munb <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
.lsize , earg = .esize )
rzanegbin(nsim * length(munb),
pobs0 = phi0, munb = munb, size = kmat)
}, list( .lpobs0 = lpobs0, .lmunb = lmunb, .lsize = lsize,
.epobs0 = epobs0, .emunb = emunb, .esize = esize ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
M1 <- 3
NOS <- ncol(eta) / M1
phi0 <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lpobs0 , earg = .epobs0 )
munb <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lmunb , earg = .emunb )
size <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lsize , earg = .esize )
okay1 <- all(is.finite(munb)) && all(0 < munb) &&
all(is.finite(size)) && all(0 < size) &&
all(is.finite(phi0)) && all(0 < phi0 & phi0 < 1)
smallval <- .mds.min # .munb.div.size
overdispersion <- if (okay1)
all(smallval < munb / size) else FALSE
if (!overdispersion)
warning("parameter 'size' has very large values; ",
"try fitting a zero-altered Poisson ",
"model instead.")
okay1 && overdispersion
}, list( .lpobs0 = lpobs0, .lmunb = lmunb, .lsize = lsize,
.epobs0 = epobs0, .emunb = emunb, .esize = esize,
.mds.min = mds.min))),
deriv = eval(substitute(expression({
M1 <- 3
NOS <- ncol(eta) / M1
y0 <- extra$y0
phi0 <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lpobs0 , earg = .epobs0 )
munb <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lsize , earg = .esize )
skip <- extra$skip.these
dphi0.deta <- dtheta.deta(phi0, .lpobs0 , .epobs0 )
dmunb.deta <- dtheta.deta(munb, .lmunb , .emunb )
dsize.deta <- dtheta.deta(kmat, .lsize , .esize )
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
if (FALSE)
warning("parameter 'size' has very large values; ",
"try fitting a zero-altered Poisson ",
"model instead")
kmat[big.size] <- munb[big.size] / smallval
}
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
tempm <- munb / (kmat + munb)
prob0 <- tempk^kmat
oneminusf0 <- 1 - prob0
AA16 <- tempm + log(tempk)
df0.dmunb <- -tempk * prob0
df0.dkmat <- prob0 * AA16
df02.dmunb2 <- prob0 * tempk * (1 +
1 / kmat) / (1 + munb / kmat)
df02.dkmat2 <- prob0 * ((tempm^2) / kmat + AA16^2)
df02.dkmat.dmunb <- -prob0 * (tempm / kmat +
AA16) / (1 + munb / kmat)
if (any(big.size)) {
prob0[big.size] <- exp(-munb[big.size]) # The limit
oneminusf0[big.size] <- -expm1(-munb[big.size])
df0.dmunb[big.size] <- -tempk[big.size] *
prob0[big.size]
df0.dkmat[big.size] <- prob0[big.size] *
AA16[big.size]
df02.dmunb2[big.size] <- prob0[big.size] *
tempk[big.size] *
(1 + 1/kmat[big.size]) / (1 + smallval)
df02.dkmat2[big.size] <- prob0[big.size] *
((tempm[big.size])^2 / kmat[big.size] +
AA16[big.size]^2)
df02.dkmat.dmunb[big.size] <- -prob0[big.size] *
(tempm[big.size]/kmat[big.size] +
AA16[big.size]) / (1+smallval)
}
mymu <- munb / oneminusf0 # E(Y) of Pos-NBD
dl.dphi0 <- -1 / (1 - phi0)
dl.dmunb <- y / munb - (1 + y/kmat) / (1 + munb/kmat) +
df0.dmunb / oneminusf0
dl.dsize <- digamma(y + kmat) - digamma(kmat) -
(y - munb) / (munb + kmat) + log(tempk) +
df0.dkmat / oneminusf0
if (any(big.size)) {
}
dl.dphi0[y == 0] <- 1 / phi0[y == 0] # Do it in one line
skip <- extra$skip.these
for (spp. in 1:NOS) {
dl.dsize[skip[, spp.], spp.] <-
dl.dmunb[skip[, spp.], spp.] <- 0
}
dl.deta23 <- c(w) * cbind(dl.dmunb * dmunb.deta,
dl.dsize * dsize.deta)
dl.deta1 <- if ( .lpobs0 == "logitlink") {
c(w) * (y0 - phi0)
} else {
c(w) * dl.dphi0 * dphi0.deta
}
ans <- cbind(dl.deta1, dl.deta23)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
ans
}), list( .lpobs0 = lpobs0 , .lmunb = lmunb , .lsize = lsize ,
.epobs0 = epobs0 , .emunb = emunb , .esize = esize,
.mds.min = mds.min ))),
weight = eval(substitute(expression({
wz <- matrix(0, n, M + M-1) # tridiagonal
max.support <- .max.support
max.chunk.MB <- .max.chunk.MB
mu.phi0 <- phi0 # pobs0 # phi0
tmp100 <- mu.phi0 * (1 - mu.phi0)
wz[, (1:NOS)*M1 - 2] <-
if ( .lpobs0 == "logitlink" && is.empty.list( .epobs0 )) {
cbind(c(w) * tmp100)
} else {
cbind(c(w) * (1 / tmp100) *
dtheta.deta(mu.phi0, .lpobs0 , .epobs0 )^2)
}
ned2l.dmunb2 <- mymu / munb^2 -
((1 + mymu/kmat) / kmat) / (1 + munb/kmat)^2 -
df02.dmunb2 / oneminusf0 -
(df0.dmunb / oneminusf0)^2
wz[, M1*(1:NOS) - 1] <- c(w) * (1 - phi0) *
ned2l.dmunb2 * dmunb.deta^2
ned2l.dmunbsize <- (munb - mymu) / (munb + kmat)^2 -
df02.dkmat.dmunb / oneminusf0 -
df0.dmunb * df0.dkmat / oneminusf0^2
wz[, M + M1*(1:NOS) - 1] <- c(w) * (1 - phi0) *
ned2l.dmunbsize * dmunb.deta * dsize.deta
ind2 <- matrix(FALSE, n, NOS) # Used for SFS
for (jay in 1:NOS) {
eff.p <- sort(c( .cutoff.prob , 1 - .cutoff.prob ))
Q.mins <- 1
Q.maxs <- qgaitdnbinom(p = eff.p[2],
truncate = 0, # prob = phi0,
kmat[, jay],
munb.p = munb[, jay]) + 10
eps.trig <- .eps.trig
Q.MAXS <- pmax(10, ceiling(1 / sqrt(eps.trig)))
Q.maxs <- pmin(Q.maxs, Q.MAXS)
ind1 <- if (max.chunk.MB > 0)
(Q.maxs - Q.mins < max.support) else
FALSE
if ((NN <- sum(ind1)) > 0) {
Object.Size <- NN * 8 * max(Q.maxs - Q.mins) / (2^20)
n.chunks <- if (intercept.only) 1 else
max(1, ceiling( Object.Size / max.chunk.MB))
chunk.rows <- ceiling(NN / n.chunks)
ind2[, jay] <- ind1 # Save this
wind2 <- which(ind1)
upr.ptr <- 0
lwr.ptr <- upr.ptr + 1
while (lwr.ptr <= NN) {
upr.ptr <- min(upr.ptr + chunk.rows, NN)
sind2 <- wind2[lwr.ptr:upr.ptr]
aaa <-
wz[sind2, M1*jay] <-
EIM.posNB.specialp(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only)
if (FALSE)
wz2[sind2, M1*jay] <-
EIM.posNB.speciald(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.min = min(Q.mins2[sind2]),
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only) # *
if (any(eim.kk.TF <- wz[sind2, M1*jay] <= 0 |
is.na(wz[sind2, M1*jay]))) {
ind2[sind2[eim.kk.TF], jay] <- FALSE
}
lwr.ptr <- upr.ptr + 1
} # while
} # if
} # end of for (jay in 1:NOS)
for (jay in 1:NOS) {
run.varcov <- 0
ii.TF <- !ind2[, jay] # Not assigned above
if (any(ii.TF)) {
kkvec <- kmat[ii.TF, jay]
muvec <- munb[ii.TF, jay]
for (ii in 1:( .nsimEIM )) {
ysim <- rzanegbin(sum(ii.TF),
munb = muvec, size = kkvec,
pobs0 = phi0[ii.TF, jay])
dl.dk <- digamma(ysim + kkvec) - digamma(kkvec) -
(ysim - muvec) / (muvec + kkvec) +
log1p(-muvec / (kkvec + muvec)) +
df0.dkmat[ii.TF, jay] / oneminusf0[ii.TF, jay]
dl.dk[ysim == 0] <- 0
run.varcov <- run.varcov + dl.dk^2
} # end of for loop
run.varcov <- c(run.varcov / .nsimEIM )
ned2l.dk2 <- if (intercept.only) mean(run.varcov) else
run.varcov
wz[ii.TF, M1*jay] <- ned2l.dk2 # * (dsize.deta[ii.TF, jay])^2
}
} # jay
wz[, M1*(1:NOS) ] <- wz[, M1*(1:NOS) ] * dsize.deta^2
save.weights <- !all(ind2)
wz[, M1*(1:NOS) ] <- c(w) * (1 - phi0) *
wz[, M1*(1:NOS) ]
wz
}), list( .lpobs0 = lpobs0,
.epobs0 = epobs0,
.cutoff.prob = cutoff.prob, .eps.trig = eps.trig,
.max.support = max.support,
.max.chunk.MB = max.chunk.MB,
.nsimEIM = nsimEIM ))))
} # zanegbinomial()
zanegbinomialff.control <- function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
zanegbinomialff <-
function(
lmunb = "loglink",
lsize = "loglink", lonempobs0 = "logitlink",
type.fitted = c("mean", "munb", "pobs0", "onempobs0"),
isize = NULL, ionempobs0 = NULL,
zero = c("size", "onempobs0"),
mds.min = 1e-3,
iprobs.y = NULL, # 0.35,
gprobs.y = (0:9)/10, # 20160709; grid for finding munb.init
cutoff.prob = 0.999, # higher is better for large 'size'
eps.trig = 1e-7,
max.support = 4000, # 20160127; I have changed this
max.chunk.MB = 30, # max.memory = Inf is allowed
gsize.mux = exp(c(-30, -20, -15, -10, -6:3)),
imethod = 1,
imunb = NULL,
nsimEIM = 500) {
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
if (!is.Numeric(eps.trig, length.arg = 1,
positive = TRUE) || eps.trig > 0.001)
stop("argument 'eps.trig' must be positive and ",
"smaller in value")
if (!is.Numeric(nsimEIM, length.arg = 1,
positive = TRUE, integer.valued = TRUE))
stop("argument 'nsimEIM' must be a positive integer")
if (nsimEIM <= 30)
warning("argument 'nsimEIM' should be greater than 30, say")
if (length(ionempobs0) &&
(!is.Numeric(ionempobs0, positive = TRUE) ||
max(ionempobs0) >= 1))
stop("If given, argument 'ionempobs0' must contain ",
"values in (0,1) only")
if (length(isize) && !is.Numeric(isize, positive = TRUE))
stop("If given, argument 'isize' must contain ",
"positive values only")
lmunb <- as.list(substitute(lmunb))
emunb <- link2list(lmunb)
lmunb <- attr(emunb, "function.name")
lsize <- as.list(substitute(lsize))
esize <- link2list(lsize)
lsize <- attr(esize, "function.name")
lonempobs0 <- as.list(substitute(lonempobs0))
eonempobs0 <- link2list(lonempobs0)
lonempobs0 <- attr(eonempobs0, "function.name")
ipobs0.small <- 1/64 # A number easily represented exactly
type.fitted <- match.arg(type.fitted,
c("mean", "munb",
"pobs0", "onempobs0"))[1]
new("vglmff",
blurb = c("Zero-altered negative binomial (Bernoulli and\n",
"positive-negative binomial conditional model)\n\n",
"Links: ",
namesof("munb", lmunb, emunb, tag = FALSE), ", ",
namesof("size", lsize, esize, tag = FALSE), ", ",
namesof("onempobs0", lonempobs0, earg = eonempobs0,
tag = FALSE), "\n",
"Mean: onempobs0 * munb / (1 - (size / (size + ",
"munb))^size)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 3,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
mds.min = .mds.min ,
multipleResponses = TRUE,
nsimEIM = .nsimEIM ,
parameters.names = c("munb", "size", "onempobs0"),
eps.trig = .eps.trig ,
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.nsimEIM = nsimEIM, .eps.trig = eps.trig,
.type.fitted = type.fitted,
.mds.min = mds.min
))),
initialize = eval(substitute(expression({
M1 <- 3
temp16 <-
w.y.check(w = w, y = y,
Is.integer.y = TRUE,
Is.nonnegative.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp16$w
y <- temp16$y
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
M <- M1 * ncoly
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("munb", NOS, skip1 = TRUE)
mynames2 <- param.names("size", NOS, skip1 = TRUE)
mynames3 <- param.names("onempobs0", NOS, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lmunb , .emunb , tag = FALSE),
namesof(mynames2, .lsize , .esize , tag = FALSE),
namesof(mynames3, .lonempobs0 , .eonempobs0 ,
tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$skip.these <-
skip.these <- matrix(as.logical(y0), n, NOS)
gprobs.y <- .gprobs.y
imunb <- .imunb # Default in NULL
if (length(imunb))
imunb <- matrix(imunb, n, NOS, byrow = TRUE)
if (!length(etastart)) {
munb.init <-
size.init <- matrix(NA_real_, n, NOS)
gprobs.y <- .gprobs.y
if (length( .iprobs.y ))
gprobs.y <- .iprobs.y
gsize.mux <- .gsize.mux # gsize.mux is on a relative scale
for (jay in 1:NOS) { # For each response 'y_jay'... do:
TFvec <- y[, jay] > 0 # Important to exclude the 0s
posyvec <- y[TFvec, jay]
munb.init.jay <- if ( .imethod == 1 ) {
quantile(posyvec, probs = gprobs.y) - 1/2 # + 1/16
} else {
weighted.mean(posyvec, w = w[TFvec, jay]) - 1/2
}
if (length(imunb))
munb.init.jay <- imunb[, jay]
gsize <- gsize.mux * 0.5 * (mean(munb.init.jay) +
weighted.mean(posyvec, w = w[TFvec, jay]))
if (length( .isize ))
gsize <- .isize # isize is on an absolute scale
try.this <-
grid.search2(munb.init.jay, gsize,
objfun = posNBD.Loglikfun2,
y = posyvec, # x = x[TFvec, , drop = FALSE],
w = w[TFvec, jay],
ret.objfun = TRUE) # Last value is the loglik
munb.init[, jay] <- try.this["Value1"]
size.init[, jay] <- try.this["Value2"]
} # for (jay ...)
pobs0.init <- matrix(if (length( .ionempobs0 ))
1 - .ionempobs0 else -1,
n, NOS, byrow = TRUE)
for (jay in 1:NOS) {
if (any(pobs0.init[, jay] < 0)) {
index.y0 <- y[, jay] < 0.5
pobs0.init[, jay] <- max(min(mean(index.y0),
1 - .ipobs0.small ),
.ipobs0.small )
}
}
etastart <-
cbind(theta2eta(munb.init , .lmunb , .emunb ),
theta2eta(size.init , .lsize , .esize ),
theta2eta(1 - pobs0.init, .lonempobs0 , .eonempobs0 ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
} # End of if (!length(etastart))
}),
list( .lonempobs0 = lonempobs0, .lmunb = lmunb, .lsize = lsize,
.eonempobs0 = eonempobs0, .emunb = emunb, .esize = esize,
.ionempobs0 = ionempobs0, .imunb = imunb, .isize = isize,
.gprobs.y = gprobs.y, .gsize.mux = gsize.mux,
.ipobs0.small = ipobs0.small,
.imethod = imethod,
.iprobs.y = iprobs.y, .type.fitted = type.fitted ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0",
"onempobs0"))[1]
M1 <- 3
NOS <- ncol(eta) / M1
munb <- eta2theta(eta[, M1*(1:NOS)-2], .lmunb , .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS)-1], .lsize , .esize )
onempobs0 <- eta2theta(eta[, M1*(1:NOS) ], .lonempobs0 ,
earg = .eonempobs0 )
tempk <- 1 / (1 + munb / kmat) # kmat/(kmat+munb); NBD p(0)
prob0 <- tempk^kmat # p(0) from negative binomial
oneminusf0 <- 1 - prob0
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
prob0[big.size] <- exp(-munb[big.size]) # Lim as
oneminusf0[big.size] <- -expm1(-munb[big.size])
}
ans <- switch(type.fitted,
"mean" = onempobs0 * munb / oneminusf0,
"munb" = munb,
"pobs0" = 1 - onempobs0, # P(Y=0)
"onempobs0" = onempobs0) # P(Y>0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
},
list( .lonempobs0 = lonempobs0, .lsize = lsize, .lmunb = lmunb,
.eonempobs0 = eonempobs0, .emunb = emunb, .esize = esize,
.mds.min = mds.min ))),
last = eval(substitute(expression({
misc$link <-
c(rep_len( .lmunb , NOS),
rep_len( .lsize , NOS),
rep_len( .lonempobs0 , NOS))[
interleave.VGAM(M1*NOS, M1 = M1)]
temp.names <- c(mynames1,
mynames2,
mynames3)[interleave.VGAM(M1*NOS, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M1*NOS)
names(misc$earg) <- temp.names
for (ii in 1:NOS) {
misc$earg[[M1*ii-2]] <- .emunb
misc$earg[[M1*ii-1]] <- .esize
misc$earg[[M1*ii ]] <- .eonempobs0
}
misc$nsimEIM <- .nsimEIM
misc$imethod <- .imethod
misc$ionempobs0 <- .ionempobs0
misc$isize <- .isize
misc$multipleResponses <- TRUE
}),
list( .lonempobs0 = lonempobs0, .lmunb = lmunb, .lsize = lsize,
.eonempobs0 = eonempobs0, .emunb = emunb, .esize = esize,
.ionempobs0 = ionempobs0, .isize = isize,
.nsimEIM = nsimEIM,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
M1 <- 3
NOS <- ncol(eta) / M1
munb <- eta2theta(eta[, M1*(1:NOS)-2], .lmunb , .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS)-1], .lsize , .esize )
onempobs0 <- eta2theta(eta[, M1*(1:NOS) ], .lonempobs0 ,
earg = .eonempobs0 )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzanegbin(x = y, pobs0 = 1 - onempobs0,
munb = munb, size = kmat, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
},
list( .lonempobs0 = lonempobs0,
.lmunb = lmunb, .lsize = lsize,
.eonempobs0 = eonempobs0,
.emunb = emunb, .esize = esize ))),
vfamily = c("zanegbinomialff"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
munb <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
.lsize , earg = .esize )
onempobs0 <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
.lonempobs0 ,
earg = .eonempobs0 )
rzanegbin(nsim * length(munb),
pobs0 = 1 - onempobs0, munb = munb, size = kmat)
},
list( .lonempobs0 = lonempobs0,
.lmunb = lmunb, .lsize = lsize,
.eonempobs0 = eonempobs0,
.emunb = emunb, .esize = esize ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
M1 <- 3
NOS <- ncol(eta) / M1
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
size <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
onempobs0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempobs0 , earg = .eonempobs0 )
okay1 <- all(is.finite(munb)) && all(0 < munb) &&
all(is.finite(size)) && all(0 < size) &&
all(is.finite(onempobs0)) &&
all(0 < onempobs0 & onempobs0 < 1)
smallval <- .mds.min # .munb.div.size
overdispersion <- if (okay1)
all(munb / size > smallval) else FALSE
if (!overdispersion)
warning("parameter 'size' has very large values; ",
"try fitting a zero-altered Poisson ",
"model instead.")
okay1 && overdispersion
},
list( .lonempobs0 = lonempobs0, .lmunb = lmunb, .lsize = lsize,
.eonempobs0 = eonempobs0, .emunb = emunb, .esize = esize,
.mds.min = mds.min))),
deriv = eval(substitute(expression({
M1 <- 3
NOS <- ncol(eta) / M1
y0 <- extra$y0
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
onempobs0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempobs0 , earg = .eonempobs0 )
skip <- extra$skip.these
phi0 <- 1 - onempobs0
dmunb.deta <- dtheta.deta(munb, .lmunb , earg = .emunb )
dsize.deta <- dtheta.deta(kmat, .lsize , earg = .esize )
donempobs0.deta <- dtheta.deta(onempobs0, .lonempobs0 ,
earg = .eonempobs0 )
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
if (FALSE)
warning("parameter 'size' has very large values; ",
"try fitting a zero-altered Poisson ",
"model instead")
kmat[big.size] <- munb[big.size] / smallval
}
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
tempm <- munb / (kmat + munb)
prob0 <- tempk^kmat
oneminusf0 <- 1 - prob0
AA16 <- tempm + log(tempk)
df0.dmunb <- -tempk * prob0
df0.dkmat <- prob0 * AA16
df02.dmunb2 <- prob0 * tempk * (1 +
1 / kmat) / (1 + munb / kmat)
df02.dkmat2 <- prob0 * ((tempm^2) / kmat + AA16^2)
df02.dkmat.dmunb <- -prob0 * (tempm / kmat +
AA16) / (1 + munb / kmat)
mymu <- munb / oneminusf0 # E(Y) of Pos-NBD
dl.dmunb <- y / munb - (1 + y/kmat) / (1 + munb/kmat) +
df0.dmunb / oneminusf0
dl.dsize <- digamma(y + kmat) - digamma(kmat) -
(y - munb) / (munb + kmat) + log(tempk) +
df0.dkmat / oneminusf0
dl.donempobs0 <- +1 / (onempobs0)
if (any(big.size)) {
}
dl.donempobs0[y == 0] <-
-1 / (1 - onempobs0[y == 0]) # Do it in 1 line
skip <- extra$skip.these
for (spp. in 1:NOS) {
dl.dsize[skip[, spp.], spp.] <-
dl.dmunb[skip[, spp.], spp.] <- 0
}
dl.deta12 <- c(w) * cbind(dl.dmunb * dmunb.deta,
dl.dsize * dsize.deta)
dl.deta3 <- if ( .lonempobs0 == "logitlink") {
-c(w) * (y0 - phi0)
} else {
-c(w) * dl.donempobs0 * donempobs0.deta
}
ans <- cbind(dl.deta12, dl.deta3)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
ans
}),
list( .lonempobs0 = lonempobs0 ,
.lmunb = lmunb , .lsize = lsize,
.eonempobs0 = eonempobs0 ,
.emunb = emunb , .esize = esize,
.mds.min = mds.min ))),
weight = eval(substitute(expression({
wz <- matrix(0, n, M + M-1) # tridiagonal
max.support <- .max.support
max.chunk.MB <- .max.chunk.MB
tmp100 <- onempobs0 * (1 - onempobs0)
wz[, (1:NOS)*M1 ] <-
if ( .lonempobs0 == "logitlink" &&
is.empty.list( .eonempobs0 )) {
cbind(c(w) * tmp100)
} else {
cbind(c(w) * (1 / tmp100) *
dtheta.deta(onempobs0, link = .lonempobs0 ,
earg = .eonempobs0 )^2)
}
ned2l.dmunb2 <- mymu / munb^2 -
((1 + mymu/kmat) / kmat) / (1 + munb/kmat)^2 -
df02.dmunb2 / oneminusf0 -
(df0.dmunb / oneminusf0)^2
wz[, M1*(1:NOS) - 2] <- c(w) * (1 - phi0) *
ned2l.dmunb2 * dmunb.deta^2
ned2l.dmunbsize <- (munb - mymu) / (munb + kmat)^2 -
df02.dkmat.dmunb / oneminusf0 -
df0.dmunb * df0.dkmat / oneminusf0^2
wz[, M + M1*(1:NOS) - 2] <- c(w) * (1 - phi0) *
ned2l.dmunbsize * dmunb.deta * dsize.deta
ind2 <- matrix(FALSE, n, NOS) # Used for SFS
for (jay in 1:NOS) {
eff.p <- sort(c( .cutoff.prob , 1 - .cutoff.prob ))
Q.mins <- 1
Q.maxs <- qgaitdnbinom(p = eff.p[2], truncate = 0, # prob = phi0,
kmat[, jay], munb.p = munb[, jay]) + 10
eps.trig <- .eps.trig
Q.MAXS <- pmax(10, ceiling(1 / sqrt(eps.trig)))
Q.maxs <- pmin(Q.maxs, Q.MAXS)
ind1 <- if (max.chunk.MB > 0)
(Q.maxs - Q.mins < max.support) else FALSE
if ((NN <- sum(ind1)) > 0) {
Object.Size <- NN * 8 * max(Q.maxs - Q.mins) / (2^20)
n.chunks <- if (intercept.only) 1 else
max(1, ceiling( Object.Size / max.chunk.MB))
chunk.rows <- ceiling(NN / n.chunks)
ind2[, jay] <- ind1 # Save this
wind2 <- which(ind1)
upr.ptr <- 0
lwr.ptr <- upr.ptr + 1
while (lwr.ptr <= NN) {
upr.ptr <- min(upr.ptr + chunk.rows, NN)
sind2 <- wind2[lwr.ptr:upr.ptr]
wz[sind2, M1*jay - 1] <-
EIM.posNB.specialp(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only)
if (FALSE)
wz2[sind2, M1*jay - 1] <-
EIM.posNB.speciald(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.min = min(Q.mins2[sind2]),
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only) # *
if (any(eim.kk.TF <- wz[sind2, M1*jay - 1] <= 0 |
is.na(wz[sind2, M1*jay - 1]))) {
ind2[sind2[eim.kk.TF], jay] <- FALSE
}
lwr.ptr <- upr.ptr + 1
} # while
} # if
} # end of for (jay in 1:NOS)
for (jay in 1:NOS) {
run.varcov <- 0
ii.TF <- !ind2[, jay] # Not assigned above
if (any(ii.TF)) {
kkvec <- kmat[ii.TF, jay]
muvec <- munb[ii.TF, jay]
for (ii in 1:( .nsimEIM )) {
ysim <- rzanegbin(sum(ii.TF), munb = muvec,
size = kkvec,
pobs0 = phi0[ii.TF, jay])
dl.dk <- digamma(ysim + kkvec) - digamma(kkvec) -
(ysim - muvec) / (muvec + kkvec) +
log1p(-muvec / (kkvec + muvec)) +
df0.dkmat[ii.TF, jay] / oneminusf0[ii.TF, jay]
dl.dk[ysim == 0] <- 0
run.varcov <- run.varcov + dl.dk^2
} # end of for loop
run.varcov <- c(run.varcov / .nsimEIM )
ned2l.dk2 <- if (intercept.only)
mean(run.varcov) else run.varcov
wz[ii.TF, M1*jay - 1] <- ned2l.dk2 # *
}
} # jay
wz[, M1*(1:NOS) - 1] <- wz[, M1*(1:NOS) - 1] * dsize.deta^2
save.weights <- !all(ind2)
wz[, M1*(1:NOS) - 1] <- c(w) * (1 - phi0) *
wz[, M1*(1:NOS) - 1]
wz
}),
list( .lonempobs0 = lonempobs0,
.eonempobs0 = eonempobs0,
.cutoff.prob = cutoff.prob, .eps.trig = eps.trig,
.max.support = max.support,
.max.chunk.MB = max.chunk.MB,
.nsimEIM = nsimEIM ))))
} # End of zanegbinomialff()
zipoisson <-
function(lpstr0 = "logitlink", llambda = "loglink",
type.fitted = c("mean", "lambda", "pobs0",
"pstr0", "onempstr0"),
ipstr0 = NULL, ilambda = NULL,
gpstr0 = NULL, # (1:9) / 10,
imethod = 1,
ishrinkage = 0.95, probs.y = 0.35,
parallel = FALSE, # Added 20171223
zero = NULL) {
ipstr00 <- ipstr0
gpstr00 <- gpstr0
ipstr0.small <- 1/64 # A number easily represented exactly
lpstr0 <- as.list(substitute(lpstr0))
epstr00 <- link2list(lpstr0)
lpstr00 <- attr(epstr00, "function.name")
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0",
"pstr0", "onempstr0"))[1]
if (length(ipstr00))
if (!is.Numeric(ipstr00, positive = TRUE) ||
any(ipstr00 >= 1))
stop("argument 'ipstr0' values must be inside the ",
"interval (0, 1)")
if (length(ilambda))
if (!is.Numeric(ilambda, positive = TRUE))
stop("argument 'ilambda' values must be positive")
new("vglmff",
blurb = c("Zero-inflated Poisson\n\n",
"Links: ",
namesof("pstr0", lpstr00, earg = epstr00 ), ", ",
namesof("lambda", llambda, earg = elambda ), "\n",
"Mean: (1 - pstr0) * lambda"),
constraints = eval(substitute(expression({
constraints <- cm.VGAM(matrix(1, M, 1), x = x,
bool = .parallel ,
constraints = constraints)
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .parallel = parallel, .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
hadof = TRUE,
multipleResponses = TRUE,
parallel = .parallel ,
parameters.names = c("pstr0", "lambda"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted,
.parallel = parallel
))),
rqresslot = eval(substitute(
function(mu, y, w, eta, extra = NULL) {
TFvec <- c(TRUE, FALSE)
phimat <- eta2theta(eta[, TFvec], .lpstr00 , .epstr00 )
lambda <- eta2theta(eta[, !TFvec], .llambda , .elambda )
scrambleseed <- runif(1) # To scramble the seed
qnorm(runif(length(y),
pzipois(y - 1, pstr0 = phimat, lambda),
pzipois(y , pstr0 = phimat, lambda)))
}, list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda ))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("pstr0", ncoly, skip1 = TRUE)
mynames2 <- param.names("lambda", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lpstr00 , .epstr00 , tag = FALSE),
namesof(mynames2, .llambda , .elambda , tag = FALSE))[
interleave.VGAM(M, M1 = M1)]
if (!length(etastart)) {
matL <- Init.mu(y = y, w = w, imethod = .imethod ,
imu = .ilambda ,
ishrinkage = .ishrinkage ,
pos.only = TRUE, # x = x,
probs.y = .probs.y )
matP <- matrix(if (length( .ipstr00 )) .ipstr00 else 0,
n, ncoly, byrow = TRUE)
phi.grid <- .gpstr00 # seq(0.02, 0.98, len = 21)
ipstr0.small <- .ipstr0.small # Easily represented exactly
if (!length( .ipstr00 ))
for (jay in 1:ncoly) {
zipois.Loglikfun <-
function(phival, y, x, w, extraargs) {
sum(c(w) * dzipois(y, pstr0 = phival,
lambda = extraargs$lambda,
log = TRUE))
}
Phi.init <- if (length(phi.grid)) {
grid.search(phi.grid, objfun = zipois.Loglikfun,
y = y[, jay], w = w[, jay], # x = x,
extraargs = list(lambda = matL[, jay]))
} else {
pmax(ipstr0.small,
weighted.mean(y[, jay] == 0, w[, jay]) -
dpois(0, matL[, jay]))
}
if (mean(Phi.init == ipstr0.small) > 0.95 &&
.lpstr00 != "identitylink")
warning("from the initial values only, the data ",
"appears to have little or no 0-inflation,",
" and possibly 0-deflation.")
matP[, jay] <- Phi.init
} # for (jay)
etastart <- cbind(theta2eta(matP, .lpstr00 , .epstr00 ),
theta2eta(matL, .llambda , .elambda ))[,
interleave.VGAM(M, M1 = M1)]
mustart <- NULL # Since etastart has been computed.
} # End of !length(etastart)
}), list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda,
.ipstr00 = ipstr00, .ilambda = ilambda,
.gpstr00 = gpstr00,
.imethod = imethod, .probs.y = probs.y,
.ipstr0.small = ipstr0.small,
.type.fitted = type.fitted,
.ishrinkage = ishrinkage ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0",
"pstr0", "onempstr0"))[1]
phimat <- eta2theta(eta[, c(TRUE, FALSE)],
.lpstr00 , .epstr00 )
lambda <- eta2theta(eta[, c(FALSE, TRUE)],
.llambda , .elambda )
ans <- switch(type.fitted,
"mean" = (1 - phimat) * lambda,
"lambda" = lambda,
"pobs0" = phimat + (1-phimat)*exp(-lambda), # P(Y=0)
"pstr0" = phimat,
"onempstr0" = 1 - phimat)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda
))),
last = eval(substitute(expression({
M1 <- extra$M1
misc$link <-
c(rep_len( .lpstr00 , ncoly),
rep_len( .llambda , ncoly))[interleave.VGAM(M, M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M,
M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- .epstr00
misc$earg[[M1*ii ]] <- .elambda
}
misc$M1 <- M1
misc$imethod <- .imethod
misc$expected <- TRUE
misc$multipleResponses <- TRUE
if (FALSE) {
misc$pobs0 <- phimat + (1 - phimat) * exp(-lambda) # P(Y=0)
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pobs0) <- dimnames(y)
misc$pstr0 <- phimat
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pstr0) <- dimnames(y)
}
}), list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
phimat <- eta2theta(eta[, c(TRUE, FALSE)],
.lpstr00 , .epstr00 )
lambda <- eta2theta(eta[, c(FALSE, TRUE)],
.llambda , .elambda )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) * dzipois(y, pstr0 = phimat,
lambda = lambda, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda ))),
vfamily = c("zipoisson"),
hadof = eval(substitute(
function(eta, extra = list(),
linpred.index = 1, w = 1,
dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2),
deriv = 1, ...) {
M1 <- 2
n <- NROW(eta)
M <- NCOL(eta)
phimat <- eta2theta(eta[, c(TRUE, FALSE)],
.lpstr00 , .epstr00 )
lambda <- eta2theta(eta[, c(FALSE, TRUE)],
.llambda , .elambda )
which.param <- ifelse(linpred.index %% M1 == 1,
"phi", "lambda")
which.y <- ceiling(linpred.index / M1)
prob0 <- exp(-lambda)
pobs0 <- phimat + (1 - phimat) * prob0
if (deriv == 0) {
ned2l.dphimat2 <- -expm1(-lambda) / ((1 - phimat) * pobs0)
ned2l.dphimatlambda <- -exp(-lambda) / pobs0
ned2l.dlambda2 <- (1 - phimat) / lambda -
phimat * (1 - phimat) * exp(-lambda) / pobs0
wz <- array(c(c(w) * ned2l.dphimat2,
c(w) * ned2l.dlambda2,
c(w) * ned2l.dphimatlambda),
dim = c(n, M / M1, 3))
return(arwz2wz(wz, M = M, M1 = M1, full.arg = TRUE))
}
if (which.param == "phi") {
NED2l.dphimat2 <-
+expm1(-lambda) * (1 - 2 * pobs0) / ((1 - phimat) * pobs0)^2
NED2l.dphimatlambda <- -exp(-lambda) *
expm1(-lambda) / pobs0^2
NED2l.dlambda2 <- -1 / lambda -
exp(-lambda) * ((1 - phimat)^2 * exp(-lambda) -
phimat^2) / pobs0^2
} else {
NED2l.dphimat2 <- exp(-lambda) / ((1 - phimat) * pobs0^2)
NED2l.dphimatlambda <- phimat * exp(-lambda) / pobs0^2
NED2l.dlambda2 <- -(1 - phimat) / lambda^2 +
phimat^2 * (1 - phimat) * exp(-lambda) / pobs0^2
}
if (deriv == 2)
NED2l.dphimat2 <-
NED2l.dphimatlambda <-
NED2l.dlambda2 <- matrix(NA_real_, n, M)
WZ <- switch(as.character(deriv),
"1" =
array(c(c(w) * retain.col(NED2l.dphimat2, which.y),
c(w) * retain.col(NED2l.dlambda2, which.y),
c(w) * retain.col(NED2l.dphimatlambda, which.y)),
dim = c(n, M / M1, 3)),
"2" =
array(c(c(w) * retain.col(NED2l.dphimat2, which.y),
c(w) * retain.col(NED2l.dlambda2, which.y),
c(w) * retain.col(NED2l.dphimatlambda, which.y)),
dim = c(n, M / M1, 3)),
stop("argument 'deriv' must be 0 or 1 or 2"))
return(arwz2wz(WZ, M = M, M1 = M1, full.arg = TRUE))
}, list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
phimat <- eta2theta(eta[, c(TRUE, FALSE)],
.lpstr00 , .epstr00 )
lambda <- eta2theta(eta[, c(FALSE, TRUE)],
.llambda , .elambda )
okay1 <- all(is.finite(lambda)) && all(0 < lambda) &&
all(is.finite(phimat)) && all(phimat < 1)
deflat.limit <- -1 / expm1(lambda)
okay2.deflat <- TRUE
if (okay1 && !(okay2.deflat <- all(deflat.limit < phimat)))
warning("parameter 'pstr0' is too negative even ",
"allowing for 0-deflation.")
okay1 && okay2.deflat
}, list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
phimat <- eta2theta(eta[, c(TRUE, FALSE)],
.lpstr00 , .epstr00 )
lambda <- eta2theta(eta[, c(FALSE, TRUE)],
.llambda , .elambda )
rzipois(nsim * length(lambda), lambda = lambda,
pstr0 = phimat)
}, list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda ))),
deriv = eval(substitute(expression({
M1 <- 2
phimat <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
.lpstr00 , earg = .epstr00 )
lambda <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
.llambda , earg = .elambda )
prob0 <- exp(-lambda)
pobs0 <- phimat + (1 - phimat) * prob0
index0 <- as.matrix(y == 0)
dl.dphimat <- -expm1(-lambda) / pobs0
dl.dphimat[!index0] <- -1 / (1 - phimat[!index0])
dl.dlambda <- -(1 - phimat) * exp(-lambda) / pobs0
dl.dlambda[!index0] <- (y[!index0] -
lambda[!index0]) / lambda[!index0]
dphimat.deta <- dtheta.deta(phimat, .lpstr00 , .epstr00 )
dlambda.deta <- dtheta.deta(lambda, .llambda , .elambda )
ans <- c(w) * cbind(dl.dphimat * dphimat.deta,
dl.dlambda * dlambda.deta)
ans <- ans[, interleave.VGAM(M, M1 = M1)]
if ( .llambda == "loglink" && is.empty.list( .elambda ) &&
any(lambda[!index0] < .Machine$double.eps)) {
for (spp. in 1:(M / M1)) {
ans[!index0[, spp.], M1 * spp.] <-
w[!index0[, spp.]] *
(y[!index0[, spp.], spp.] -
lambda[!index0[, spp.], spp.])
}
}
ans
}), list( .lpstr00 = lpstr00, .llambda = llambda,
.epstr00 = epstr00, .elambda = elambda ))),
weight = eval(substitute(expression({
ned2l.dphimat2 <- -expm1(-lambda) / ((1 - phimat) * pobs0)
ned2l.dphimatlambda <- -exp(-lambda) / pobs0
ned2l.dlambda2 <- (1 - phimat) / lambda -
phimat * (1 - phimat) * exp(-lambda) / pobs0
wz <- array(c(c(w) * ned2l.dphimat2 * dphimat.deta^2,
c(w) * ned2l.dlambda2 * dlambda.deta^2,
c(w) * ned2l.dphimatlambda *
dphimat.deta * dlambda.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .llambda = llambda, .elambda = elambda ))))
} # zipoisson
zipoissonff <-
function(llambda = "loglink",
lonempstr0 = "logitlink",
type.fitted = c("mean", "lambda", "pobs0",
"pstr0", "onempstr0"),
ilambda = NULL, ionempstr0 = NULL,
gonempstr0 = NULL, # (1:9) / 10,
imethod = 1,
ishrinkage = 0.95, probs.y = 0.35,
zero = "onempstr0") {
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0",
"pstr0", "onempstr0"))[1]
llambda <- as.list(substitute(llambda))
elambda <- link2list(llambda)
llambda <- attr(elambda, "function.name")
lonempstr0 <- as.list(substitute(lonempstr0))
eonempstr0 <- link2list(lonempstr0)
lonempstr0 <- attr(eonempstr0, "function.name")
ipstr0.small <- 1/64 # A number easily represented exactly
if (length(ilambda))
if (!is.Numeric(ilambda, positive = TRUE))
stop("'ilambda' values must be positive")
if (length(ionempstr0))
if (!is.Numeric(ionempstr0, positive = TRUE) ||
any(ionempstr0 >= 1))
stop("'ionempstr0' values must be inside ",
"the interval (0, 1)")
new("vglmff",
blurb = c("Zero-inflated Poisson\n\n",
"Links: ",
namesof("lambda", llambda, elambda), ", ",
namesof("onempstr0", lonempstr0, eonempstr0), "\n",
"Mean: onempstr0 * lambda"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
hadof = TRUE,
multipleResponses = TRUE,
parameters.names = c("lambda", "onempstr0"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted
))),
rqresslot = eval(substitute(
function(mu, y, w, eta, extra = NULL) {
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
scrambleseed <- runif(1) # To scramble the seed
qnorm(runif(length(y),
pzipois(y - 1, pstr0 = 1 - onempstr0, lambda),
pzipois(y , pstr0 = 1 - onempstr0, lambda)))
}, list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
ncoly <- ncol(y)
extra$ncoly <- ncoly
extra$M1 <- M1
M <- M1 * ncoly
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("lambda", ncoly, skip1 = TRUE)
mynames2 <- param.names("onempstr0", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .llambda , .elambda , tag = FALSE),
namesof(mynames2, .lonempstr0 , .eonempstr0 , tag = FALSE))[
interleave.VGAM(M, M1 = M1)]
if (!length(etastart)) {
matL <- Init.mu(y = y, w = w, imethod = .imethod ,
imu = .ilambda , ishrinkage = .ishrinkage ,
pos.only = TRUE, # x = x,
probs.y = .probs.y )
matP <- matrix(if (length( .ionempstr0 ))
.ionempstr0 else 0,
n, ncoly, byrow = TRUE)
phi0.grid <- .gonempstr0
ipstr0.small <- .ipstr0.small # Easily
if (!length( .ionempstr0 ))
for (jay in 1:ncoly) {
zipois.Loglikfun <-
function(phival, y, x, w, extraargs) {
sum(c(w) * dzipois(y, pstr0 = phival,
lambda = extraargs$lambda,
log = TRUE))
}
Phi0.init <- if (length(phi0.grid)) {
grid.search(phi0.grid,
objfun = zipois.Loglikfun,
y = y[, jay], x = x, w = w[, jay],
extraargs = list(lambda = matL[, jay]))
} else {
pmax(ipstr0.small,
weighted.mean(y[, jay] == 0, w[, jay]) -
dpois(0, matL[, jay]))
}
if (mean(Phi0.init == ipstr0.small) > 0.95 &&
.lonempstr0 != "identitylink")
warning("from the initial values only, the data ",
"appears to have little or no 0-inflation,",
" and possibly 0-deflation.")
matP[, jay] <- Phi0.init
} # for (jay)
etastart <-
cbind(theta2eta( matL, .llambda , .elambda ),
theta2eta(1 - matP, .lonempstr0 , .eonempstr0 ))[,
interleave.VGAM(M, M1 = M1)]
mustart <- NULL # Since etastart has been computed.
}
}), list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda,
.ionempstr0 = ionempstr0, .ilambda = ilambda,
.gonempstr0 = gonempstr0,
.type.fitted = type.fitted, .probs.y = probs.y,
.ipstr0.small = ipstr0.small,
.imethod = imethod, .ishrinkage = ishrinkage ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <-
if (length(extra$type.fitted)) extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "lambda", "pobs0",
"pstr0", "onempstr0"))[1]
M1 <- 2
NOS <- ncoly <- ncol(eta) / M1
lambda <- eta2theta(eta[, M1*(1:ncoly) - 1], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, M1*(1:ncoly) ], .lonempstr0 ,
earg = .eonempstr0 )
ans <- switch(type.fitted,
"mean" = onempstr0 * lambda,
"lambda" = lambda,
"pobs0" = 1 + onempstr0 * expm1(-lambda), # P(Y=0)
"pstr0" = 1 - onempstr0,
"onempstr0" = onempstr0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
last = eval(substitute(expression({
M1 <- extra$M1
misc$link <-
c(rep_len( .llambda , ncoly),
rep_len( .lonempstr0 , ncoly))[interleave.VGAM(M,
M1 = M1)]
temp.names <- c(mynames1, mynames2)[interleave.VGAM(M,
M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M1 * ncoly)
names(misc$earg) <- temp.names
for (ii in 1:ncoly) {
misc$earg[[M1*ii-1]] <- .elambda
misc$earg[[M1*ii ]] <- .eonempstr0
}
misc$M1 <- M1
misc$imethod <- .imethod
if (FALSE) {
misc$pobs0 <- (1 - onempstr0) +
onempstr0 * exp(-lambda) # P(Y=0)
misc$pobs0 <- as.matrix(misc$pobs0)
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pobs0) <- dimnames(y)
misc$pstr0 <- (1 - onempstr0)
misc$pstr0 <- as.matrix(misc$pstr0)
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pstr0) <- dimnames(y)
}
}), list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <- c(w) *
dzipois(y, pstr0 = 1 - onempstr0, lambda, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
vfamily = c("zipoissonff"),
hadof = eval(substitute(
function(eta, extra = list(),
linpred.index = 1, w = 1,
dim.wz = c(NROW(eta), NCOL(eta) * (NCOL(eta)+1)/2),
deriv = 1, ...) {
M1 <- 2
n <- NROW(eta)
M <- NCOL(eta)
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
namevec <- c("lambda", "onempstr0")
whichj <- 2 - (linpred.index %% M1) # \in 1:M1
which.param <- namevec[whichj]
which.y <- ceiling(linpred.index / M1)
if (deriv == 0) {
denom <- 1 + onempstr0 * expm1(-lambda)
ned2l.dlambda2 <- ( onempstr0) / lambda -
onempstr0 * (1 - onempstr0) * exp(-lambda) / denom
ned2l.donempstr0.2 <- -expm1(-lambda) / ((onempstr0) *
denom)
ned2l.dphilambda <- +exp(-lambda) / denom
wz <- array(c(c(w) * ned2l.dlambda2,
c(w) * ned2l.donempstr0.2,
c(w) * ned2l.dphilambda),
dim = c(n, M / M1, 3))
return(arwz2wz(wz, M = M, M1 = M1, full.arg = TRUE))
}
d3.11 <- eval(
deriv3( ~ onempstr0 / lambda -
onempstr0 * (1 - onempstr0) * exp(-lambda) /
(1 + onempstr0 * expm1(-lambda)),
which.param, hessian = deriv == 2))
d3.22 <- eval(
deriv3( ~ -expm1(-lambda) / (onempstr0 *
(1 + onempstr0 * expm1(-lambda))),
which.param, hessian = deriv == 2))
d3.12 <- eval(
deriv3( ~ exp(-lambda) / (1 + onempstr0 * expm1(-lambda)),
which.param, hessian = deriv == 2))
Dl.dlambda <- matrix(attr(d3.11, "gradient"), n, M/M1)
Dl.donempstr0 <- matrix(attr(d3.22, "gradient"), n, M/M1)
Dl.dphilambda <- matrix(attr(d3.12, "gradient"), n, M/M1)
if (deriv == 2) {
NED2l.dlambda2 <-
matrix(attr(d3.11, "hessian"), n, M/M1)
NED2l.donempstr0.2 <-
matrix(attr(d3.22, "hessian"), n, M/M1)
NED2l.dphilambda <-
matrix(attr(d3.12, "hessian"), n, M/M1)
}
WZ <- switch(as.character(deriv),
"1" =
array(c(c(w) * retain.col(Dl.dlambda, which.y),
c(w) * retain.col(Dl.donempstr0, which.y),
c(w) * retain.col(Dl.dphilambda, which.y)),
dim = c(n, M / M1, 3)),
"2" =
array(c(c(w) * retain.col(NED2l.dlambda2, which.y),
c(w) * retain.col(NED2l.donempstr0.2, which.y),
c(w) * retain.col(NED2l.dphilambda, which.y)),
dim = c(n, M / M1, 3)),
stop("argument 'deriv' must be 0 or 1 or 2"))
return(arwz2wz(WZ, M = M, M1 = M1, full.arg = TRUE))
}, list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
rzipois(nsim * length(lambda), lambda, pstr0 = 1 - onempstr0)
},
list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
okay1 <- all(is.finite(lambda )) && all(0 < lambda ) &&
all(is.finite(onempstr0)) && all(0 < onempstr0)
deflat.limit <- -1 / expm1(lambda)
okay2.deflat <- TRUE
if (okay1 &&
!(okay2.deflat <- all(onempstr0 < 1 - deflat.limit)))
warning("parameter 'onempstr0' is too positive even ",
"allowing for 0-deflation.")
okay1 && okay2.deflat
}, list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
deriv = eval(substitute(expression({
M1 <- 2
ncoly <- ncol(eta) / M1 # extra$ncoly
lambda <- eta2theta(eta[, c(TRUE, FALSE)], .llambda ,
earg = .elambda )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
dlambda.deta <- dtheta.deta(lambda , .llambda ,
earg = .elambda )
donempstr0.deta <- dtheta.deta(onempstr0, .lonempstr0 ,
earg = .eonempstr0 )
denom <- 1 + onempstr0 * expm1(-lambda)
ind0 <- (y == 0)
dl.dlambda <- -onempstr0 * exp(-lambda) / denom
dl.dlambda[!ind0] <- (y[!ind0] -
lambda[!ind0]) / lambda[!ind0]
dl.donempstr0 <- expm1(-lambda) / denom
dl.donempstr0[!ind0] <- 1 / onempstr0[!ind0]
ans <- c(w) * cbind(dl.dlambda * dlambda.deta,
dl.donempstr0 * donempstr0.deta)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
if ( .llambda == "loglink" && is.empty.list( .elambda ) &&
any(lambda[!ind0] < .Machine$double.eps)) {
for (spp. in 1:ncoly) {
ans[!ind0[, spp.], M1 * spp.] <-
w[!ind0[, spp.]] *
(y[!ind0[, spp.], spp.] - lambda[!ind0[, spp.], spp.])
}
}
ans
}), list( .lonempstr0 = lonempstr0, .llambda = llambda,
.eonempstr0 = eonempstr0, .elambda = elambda ))),
weight = eval(substitute(expression({
ned2l.dlambda2 <- ( onempstr0) / lambda -
onempstr0 * (1 - onempstr0) * exp(-lambda) / denom
ned2l.donempstr0.2 <- -expm1(-lambda) / ((onempstr0) * denom)
ned2l.dphilambda <- +exp(-lambda) / denom
wz <- array(c(c(w) * ned2l.dlambda2 * dlambda.deta^2,
c(w) * ned2l.donempstr0.2 * donempstr0.deta^2,
c(w) * ned2l.dphilambda * donempstr0.deta *
dlambda.deta),
dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .llambda = llambda ))))
} # zipoissonff
zibinomial <-
function(lpstr0 = "logitlink", lprob = "logitlink",
type.fitted = c("mean", "prob", "pobs0",
"pstr0", "onempstr0"),
ipstr0 = NULL,
zero = NULL, # 20130917; was originally zero = 1,
multiple.responses = FALSE, imethod = 1) {
if (as.logical(multiple.responses))
stop("argument 'multiple.responses' must be FALSE")
lpstr0 <- as.list(substitute(lpstr0))
epstr0 <- link2list(lpstr0)
lpstr0 <- attr(epstr0, "function.name")
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0",
"pstr0", "onempstr0"))[1]
if (is.Numeric(ipstr0))
if (!is.Numeric(ipstr0, positive = TRUE) || any(ipstr0 >= 1))
stop("'ipstr0' values must be inside the interval (0,1)")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
new("vglmff",
blurb = c("Zero-inflated binomial\n\n",
"Links: ",
namesof("pstr0", lpstr0, earg = epstr0), ", ",
namesof("prob" , lprob , earg = eprob ), "\n",
"Mean: (1 - pstr0) * prob"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = NA,
type.fitted = .type.fitted ,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("pstr0", "prob"),
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted
))),
initialize = eval(substitute(expression({
if (!all(w == 1))
extra$orig.w <- w
if (NCOL(y) == 1) {
if (is.factor(y))
y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + w * y) / (1.0 + w)
no.successes <- y
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(no.successes - round(no.successes)) > 1.0e-8))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(y - round(y)) > 1.0e-8))
stop("Count data must be integer-valued")
y <- round(y)
nvec <- y[, 1] + y[, 2]
y <- ifelse(nvec > 0, y[, 1] / nvec, 0)
w <- w * nvec
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + nvec * y) / (1 + nvec)
} else {
stop("for the binomialff family, response 'y' must be a ",
"vector of 0 and 1's\n",
"or a factor ",
"(first level = fail, other levels = success),\n",
"or a 2-column matrix where col 1 is the no. of ",
"successes and col 2 is the no. of failures")
}
if ( .imethod == 1)
mustart <- (mustart + y) / 2
if ( .imethod == 2)
mustart <- mean(mustart) + 0 * y
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
predictors.names <-
c(namesof("pstr0", .lpstr0 , earg = .epstr0 , tag = FALSE),
namesof("prob" , .lprob , earg = .eprob , tag = FALSE))
extra$w <- w # Needed for @linkinv
phi.init <- if (length( .ipstr0 )) .ipstr0 else {
prob0.est <- sum(w[y == 0]) / sum(w)
if ( .imethod == 1) {
(prob0.est - (1 - mustart)^w) / (1 - (1 - mustart)^w)
} else {
prob0.est
}
}
phi.init[phi.init <= 0.05] <- 0.05 # Last resort
phi.init[phi.init >= 0.95] <- 0.95 # Last resort
if ( length(mustart) && !length(etastart))
mustart <- cbind(rep_len(phi.init, n),
mustart) # 1st coln not a real mu
}), list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob,
.ipstr0 = ipstr0,
.type.fitted = type.fitted,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
pstr0 <- eta2theta(eta[, 1], .lpstr0 , earg = .epstr0 )
mubin <- eta2theta(eta[, 2], .lprob , earg = .eprob )
orig.w <- if (length(tmp3 <- extra$orig.w)) tmp3 else
rep_len(1, nrow(eta))
priorw <- extra$w
nvec <- priorw / orig.w
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "pstr0",
"onempstr0"))[1]
ans <- switch(type.fitted,
"mean" = (1 - pstr0) * mubin,
"prob" = mubin,
"pobs0" = pstr0 + (1-pstr0)*(1-mubin)^nvec, # P(Y=0)
"pstr0" = pstr0,
"onempstr0" = 1 - pstr0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob ))),
last = eval(substitute(expression({
misc$link <- c("pstr0" = .lpstr0 , "prob" = .lprob )
misc$earg <- list("pstr0" = .epstr0 , "prob" = .eprob )
misc$imethod <- .imethod
}), list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob,
.imethod = imethod ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
cbind(theta2eta(mu[, 1], .lpstr0 , earg = .epstr0 ),
theta2eta(mu[, 2], .lprob , earg = .eprob ))
}, list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
pstr0 <- eta2theta(eta[, 1], .lpstr0 , earg = .epstr0 )
mubin <- eta2theta(eta[, 2], .lprob , earg = .eprob )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
dzibinom(x = round(w * y), size = w, prob = mubin,
log = TRUE, pstr0 = pstr0)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob ))),
vfamily = c("zibinomial"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
pstr0 <- eta2theta(eta[, 1], .lpstr0 , earg = .epstr0 )
probb <- eta2theta(eta[, 2], .lprob , earg = .eprob )
size <- extra$w
okay1 <- all(is.finite(probb)) && all(0 < probb) &&
all(is.finite(pstr0)) && all(pstr0 < 1)
prob0 <- (1 - probb)^size
Prob0.check <- dbinom(0, size = size, prob = probb)
deflat.limit <- -prob0 / (1 - prob0)
okay2.deflat <- TRUE
if (okay1 && !(okay2.deflat <- all(deflat.limit < pstr0)))
warning("parameter 'pstr0' is too negative even ",
"allowing for 0-deflation.")
okay1 && okay2.deflat
}, list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob ))),
deriv = eval(substitute(expression({
phi <- eta2theta(eta[, 1], .lpstr0 , earg = .epstr0 )
mubin <- eta2theta(eta[, 2], .lprob , earg = .eprob )
prob0 <- (1 - mubin)^w # Actually q^w
pobs0 <- phi + (1 - phi) * prob0
index <- (y == 0)
dl.dphi <- (1 - prob0) / pobs0
dl.dphi[!index] <- -1 / (1 - phi[!index])
dl.dmubin <- -w * (1 - phi) * (1 - mubin)^(w - 1) / pobs0
dl.dmubin[!index] <- w[!index] *
( y[!index] / mubin[!index] -
(1 - y[!index]) / (1 - mubin[!index]))
dphi.deta <- dtheta.deta(phi, .lpstr0 , earg = .epstr0 )
dmubin.deta <- dtheta.deta(mubin, .lprob , earg = .eprob )
ans <- cbind(dl.dphi * dphi.deta,
dl.dmubin * dmubin.deta)
if ( .lprob == "logitlink") {
ans[!index, 2] <- w[!index] * (y[!index] - mubin[!index])
}
ans
}), list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob ))),
weight = eval(substitute(expression({
wz <- matrix(NA_real_, nrow = n, ncol = dimm(M))
ned2l.dphi2 <- (1 - prob0) / ((1 - phi) * pobs0)
ned2l.dphimubin <- -w * ((1 - mubin)^(w - 1)) / pobs0
ned2l.dmubin2 <- (w * (1 - phi) / (mubin * (1 - mubin)^2)) *
(1 - mubin - w * mubin *
(1 - mubin)^w * phi / pobs0)
wz[,iam(1, 1, M)] <- ned2l.dphi2 * dphi.deta^2
wz[,iam(2, 2, M)] <- ned2l.dmubin2 * dmubin.deta^2
wz[,iam(1, 2, M)] <- ned2l.dphimubin * dphi.deta *
dmubin.deta
if (TRUE) {
ind6 <- (wz[, iam(2, 2, M)] < .Machine$double.eps)
if (any(ind6))
wz[ind6, iam(2, 2, M)] <- .Machine$double.eps
}
wz
}), list( .lpstr0 = lpstr0, .lprob = lprob,
.epstr0 = epstr0, .eprob = eprob ))))
} # zibinomial
zibinomialff <-
function(lprob = "logitlink", lonempstr0 = "logitlink",
type.fitted = c("mean", "prob", "pobs0",
"pstr0", "onempstr0"),
ionempstr0 = NULL,
zero = "onempstr0",
multiple.responses = FALSE, imethod = 1) {
if (as.logical(multiple.responses))
stop("argument 'multiple.responses' must be FALSE")
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
lonempstr0 <- as.list(substitute(lonempstr0))
eonempstr0 <- link2list(lonempstr0)
lonempstr0 <- attr(eonempstr0, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0",
"pstr0", "onempstr0"))[1]
if (is.Numeric(ionempstr0))
if (!is.Numeric(ionempstr0, positive = TRUE) ||
any(ionempstr0 >= 1))
stop("'ionempstr0' values must be inside ",
"the interval (0,1)")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
new("vglmff",
blurb = c("Zero-inflated binomial\n\n",
"Links: ",
namesof("prob" , lprob , eprob ), ", ",
namesof("onempstr0", lonempstr0, eonempstr0), "\n",
"Mean: onempstr0 * prob"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = NA,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("prob", "onempstr0"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted
))),
initialize = eval(substitute(expression({
if (!all(w == 1))
extra$orig.w <- w
if (NCOL(y) == 1) {
if (is.factor(y))
y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + w * y) / (1.0 + w)
no.successes <- y
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(no.successes - round(no.successes)) > 1.0e-8))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(y - round(y)) > 1.0e-8))
stop("Count data must be integer-valued")
y <- round(y)
nvec <- y[, 1] + y[, 2]
y <- ifelse(nvec > 0, y[, 1] / nvec, 0)
w <- w * nvec
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + nvec * y) / (1 + nvec)
} else {
stop("for the binomialff family, response 'y' ",
"must be a vector of 0 and 1's\n",
"or a factor ",
"(first level = fail, other levels = success),\n",
"or a 2-column matrix where col 1 is the no. of ",
"successes and col 2 is the no. of failures")
}
if ( .imethod == 1)
mustart <- (mustart + y) / 2
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
predictors.names <-
c(namesof("prob" ,
.lprob , .eprob , tag = FALSE),
namesof("onempstr0",
.lonempstr0 , .eonempstr0 , tag = FALSE))
extra$w <- w # Needed for @linkinv
onemphi.init <- if (length( .ionempstr0 ))
.ionempstr0 else {
prob0.est <- sum(w[y == 0]) / sum(w)
if ( .imethod == 1) {
1 - (prob0.est - (1 - mustart)^w) / (1 -
(1 - mustart)^w)
} else {
1 - prob0.est
}
}
onemphi.init[onemphi.init <= -0.10] <- 0.10 # Lots of
onemphi.init[onemphi.init <= 0.05] <- 0.15 # Last resort
onemphi.init[onemphi.init >= 0.80] <- 0.80 # Last resort
if ( length(mustart) && !length(etastart))
mustart <- cbind(mustart,
rep_len(onemphi.init, n)) # 1st coln not a real mu
}), list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob,
.ionempstr0 = ionempstr0,
.type.fitted = type.fitted,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
mubin <- eta2theta(eta[, 1], .lprob , .eprob )
onempstr0 <- eta2theta(eta[, 2], .lonempstr0 , .eonempstr0 )
orig.w <- if (length(tmp3 <- extra$orig.w)) tmp3 else
rep_len(1, nrow(eta))
priorw <- extra$w
nvec <- priorw / orig.w
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0",
"pstr0", "onempstr0"))[1]
ans <- switch(type.fitted,
"mean" = (onempstr0) * mubin,
"prob" = mubin,
"pobs0" = 1 - onempstr0 +
(onempstr0)*(1 - mubin)^nvec, # P(Y=0)
"pstr0" = 1 - onempstr0,
"onempstr0" = onempstr0)
label.cols.y(ans, colnames.y = extra$colnames.y,
NOS = NOS)
}, list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob ))),
last = eval(substitute(expression({
misc$link <-
c("prob" = .lprob , "onempstr0" = .lonempstr0 )
misc$earg <-
list("prob" = .eprob , "onempstr0" = .eonempstr0 )
misc$imethod <- .imethod
misc$pobs0 <- phi + (1 - phi) * (1 - mubin)^w # [1] # P(Y=0)
misc$pstr0 <- phi
}), list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob,
.imethod = imethod ))),
linkfun = eval(substitute(function(mu, extra = NULL) {
cbind(theta2eta(mu[, 1], .lprob , earg = .eprob ),
theta2eta(mu[, 2], .lonempstr0 , earg = .eonempstr0 ))
}, list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
mubin <- eta2theta(eta[, 1], .lprob , .eprob )
onempstr0 <- eta2theta(eta[, 2], .lonempstr0 , .eonempstr0 )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
dzibinom(x = round(w * y), size = w, prob = mubin,
log = TRUE, pstr0 = 1 - onempstr0)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob ))),
vfamily = c("zibinomialff"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
probb <- eta2theta(eta[, 1], .lprob , .eprob )
onempstr0 <- eta2theta(eta[, 2], .lonempstr0 , .eonempstr0 )
size <- extra$w
okay1 <- all(is.finite(probb)) && all(0 < probb) &&
all(is.finite(onempstr0)) && all(0 < onempstr0)
prob0 <- (1 - probb)^size
Prob0.check <- dbinom(0, size = size, prob = probb)
deflat.limit <- -prob0 / (1 - prob0)
okay2.deflat <- TRUE
if (okay1 &&
!(okay2.deflat <- all(onempstr0 < 1 - deflat.limit)))
warning("parameter 'onempstr0' is too positive even ",
"allowing for 0-deflation.")
okay1 && okay2.deflat
}, list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob ))),
deriv = eval(substitute(expression({
mubin <- eta2theta(eta[, 1], .lprob , .eprob )
onempstr0 <- eta2theta(eta[, 2], .lonempstr0 , .eonempstr0 )
omphi <- onempstr0
phi <- 1 - onempstr0
prob0 <- (1 - mubin)^w # Actually q^w
pobs0 <- phi + (omphi) * prob0
index <- (y == 0)
dl.domphi <- -(1 - prob0) / pobs0 # Note "-"
dl.domphi[!index] <- +1 / (omphi[!index]) # Note "+"
dl.dmubin <- -w * (omphi) * (1 - mubin)^(w - 1) / pobs0
dl.dmubin[!index] <- w[!index] *
( y[!index] / mubin[!index] -
(1 - y[!index]) / (1 - mubin[!index]))
dmubin.deta <- dtheta.deta(mubin, .lprob , .eprob )
domphi.deta <- dtheta.deta(omphi, .lonempstr0 , .eonempstr0 )
ans <- cbind(dl.dmubin * dmubin.deta,
dl.domphi * domphi.deta)
if ( .lprob == "logitlink") {
ans[!index, 1] <- w[!index] * (y[!index] - mubin[!index])
}
ans
}), list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob ))),
weight = eval(substitute(expression({
wz <- matrix(NA_real_, nrow = n, ncol = dimm(M))
ned2l.domphi2 <- (1 - prob0) / ((omphi) * pobs0)
ned2l.domphimubin <- +w * ((1 - mubin)^(w - 1)) / pobs0
ned2l.dmubin2 <- (w * (omphi) / (mubin * (1 - mubin)^2)) *
(1 - mubin - w * mubin *
(1 - mubin)^w * phi / pobs0)
wz[,iam(1, 1, M)] <- ned2l.dmubin2 * dmubin.deta^2
wz[,iam(2, 2, M)] <- ned2l.domphi2 * domphi.deta^2
wz[,iam(1, 2, M)] <- ned2l.domphimubin * domphi.deta *
dmubin.deta
if (TRUE) {
ind6 <- (wz[, iam(1, 1, M)] < .Machine$double.eps)
if (any(ind6))
wz[ind6, iam(1, 1, M)] <- .Machine$double.eps
}
wz
}), list( .lonempstr0 = lonempstr0, .lprob = lprob,
.eonempstr0 = eonempstr0, .eprob = eprob ))))
} # zibinomialff
dzibinom <- function(x, size, prob, pstr0 = 0, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(size),
length(prob), length(pstr0))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- dbinom(x = x, size = size, prob = prob, log = TRUE)
ans <- if (log.arg) {
ifelse(x == 0, log(pstr0 + (1-pstr0) * exp(ans)),
log1p(-pstr0) + ans)
} else {
ifelse(x == 0, pstr0 + (1-pstr0) * exp(ans) ,
(1-pstr0) * exp(ans))
}
prob0 <- (1 - prob)^size
deflat.limit <- -prob0 / (1 - prob0)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
} # dzibinom
pzibinom <- function(q, size, prob, pstr0 = 0
) {
LLL <- max(length(pstr0), length(size),
length(prob), length(q))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- pbinom(q, size, prob) # lower.tail = lower.tail,
ans <- ifelse(q < 0, 0, pstr0 + (1 - pstr0) * ans)
prob0 <- (1 - prob)^size
deflat.limit <- -prob0 / (1 - prob0)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
} # pzibinom
qzibinom <- function(p, size, prob, pstr0 = 0
) {
LLL <- max(length(p), length(size),
length(prob), length(pstr0))
p <- rep_len(p, LLL)
size <- rep_len(size, LLL)
prob <- rep_len(prob, LLL)
pstr0 <- rep_len(pstr0, LLL)
ans <- p
ans[p <= pstr0] <- 0
ans[p > pstr0] <-
qbinom((p[p > pstr0] - pstr0[p > pstr0]) / (
1 - pstr0[p > pstr0]),
size[p > pstr0],
prob[p > pstr0])
prob0 <- (1 - prob)^size
deflat.limit <- -prob0 / (1 - prob0)
ind0 <- (deflat.limit <= pstr0) & (pstr0 < 0)
if (any(ind0)) {
pobs0 <- pstr0[ind0] + (1 - pstr0[ind0]) * prob0[ind0]
ans[p[ind0] <= pobs0] <- 0
pindex <- (1:LLL)[ind0 & (p > pobs0)]
Pobs0 <- pstr0[pindex] + (1 - pstr0[pindex]) *
prob0[pindex]
ans[pindex] <- qgaitdbinom((p[pindex] - Pobs0) / (1 - Pobs0),
size[pindex], prob[pindex], truncate = 0)
}
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
} # qzibinom
rzibinom <- function(n, size, prob, pstr0 = 0) {
qzibinom(runif(n), size, prob, pstr0 = pstr0)
}
dzinegbin <-
function(x, size, prob = NULL, munb = NULL, pstr0 = 0,
log = FALSE) {
if (length(munb)) {
if (length(prob))
stop("arguments 'prob' and 'munb' both specified")
prob <- size / (size + munb)
}
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(pstr0), length(size),
length(prob), length(x))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- dnbinom(x = x, size = size, prob = prob, log = log.arg)
ans <- if (log.arg)
ifelse(x == 0, log(pstr0+(1-pstr0) * exp(ans)),
log1p(-pstr0) + ans) else
ifelse(x == 0, pstr0 + (1 - pstr0) * ans,
(1 - pstr0) * ans)
prob0 <- prob^size
deflat.limit <- -prob0 / (1 - prob0)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
pzinegbin <-
function(q, size, prob = NULL, munb = NULL, pstr0 = 0) {
if (length(munb)) {
if (length(prob))
stop("arguments 'prob' and 'munb' both specified")
prob <- size / (size + munb)
}
LLL <- max(length(pstr0), length(size),
length(prob), length(q))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- pnbinom(q = q, size = size, prob = prob)
ans <- ifelse(q < 0, 0, pstr0 + (1 - pstr0) * ans)
prob0 <- prob^size
deflat.limit <- -prob0 / (1 - prob0)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
qzinegbin <-
function(p, size, prob = NULL, munb = NULL, pstr0 = 0) {
if (length(munb)) {
if (length(prob))
stop("arguments 'prob' and 'munb' both specified")
prob <- size / (size + munb)
}
LLL <- max(length(p), length(prob), length(pstr0),
length(size))
if (length(p) != LLL) p <- rep_len(p, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
ans <- rep_len(NA_real_, LLL)
prob0 <- prob^size
deflat.limit <- -prob0 / (1 - prob0)
ans[p <= pstr0] <- 0
ind4 <- (pstr0 < p) & (deflat.limit <= pstr0)
ans[ ind4] <- qnbinom(p = (p[ind4] - pstr0[ind4]) / (
1 - pstr0[ind4]),
size = size[ind4], prob = prob[ind4])
ans[pstr0 < deflat.limit] <- NaN
ans[1 < pstr0] <- NaN
ans[p < 0] <- NaN
ans[1 < p] <- NaN
ans
}
rzinegbin <-
function(n, size, prob = NULL, munb = NULL, pstr0 = 0) {
qzinegbin(runif(n), size = size, prob = prob,
munb = munb, pstr0 = pstr0)
}
zinegbinomial.control <- function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
zinegbinomial <-
function(
zero = "size",
type.fitted = c("mean", "munb", "pobs0",
"pstr0", "onempstr0"),
mds.min = 1e-3,
nsimEIM = 500,
cutoff.prob = 0.999, # higher is better for large 'size'
eps.trig = 1e-7,
max.support = 4000, # 20160127; I have changed this
max.chunk.MB = 30, # max.memory = Inf is allowed
lpstr0 = "logitlink", lmunb = "loglink", lsize = "loglink",
imethod = 1,
ipstr0 = NULL,
imunb = NULL,
iprobs.y = NULL,
isize = NULL,
gprobs.y = (0:9)/10, # 20160709; grid for finding munb.init
gsize.mux = exp(c(-30, -20, -15, -10, -6:3))) {
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
lpstr0 <- as.list(substitute(lpstr0))
epstr0 <- link2list(lpstr0)
lpstr0 <- attr(epstr0, "function.name")
lmunb <- as.list(substitute(lmunb))
emunb <- link2list(lmunb)
lmunb <- attr(emunb, "function.name")
lsize <- as.list(substitute(lsize))
esize <- link2list(lsize)
lsize <- attr(esize, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0",
"pstr0", "onempstr0"))[1]
if (!is.Numeric(eps.trig, length.arg = 1,
positive = TRUE) || eps.trig > 0.001)
stop("argument 'eps.trig' must be positive and ",
"smaller in value")
ipstr0.small <- 1/64 # A number easily represented exactly
if (length(ipstr0) &&
(!is.Numeric(ipstr0, positive = TRUE) ||
any(ipstr0 >= 1)))
stop("argument 'ipstr0' must contain values in (0,1)")
if (length(isize) && !is.Numeric(isize, positive = TRUE))
stop("argument 'isize' must contain positive values only")
if (!is.Numeric(nsimEIM, length.arg = 1,
integer.valued = TRUE))
stop("argument 'nsimEIM' must be a positive integer")
if (nsimEIM <= 50)
warning("argument 'nsimEIM' should be greater than 50, say")
new("vglmff",
blurb = c("Zero-inflated negative binomial\n\n",
"Links: ",
namesof("pstr0", lpstr0, epstr0, tag = FALSE), ", ",
namesof("munb", lmunb, emunb, tag = FALSE), ", ",
namesof("size", lsize, esize, tag = FALSE), "\n",
"Mean: (1 - pstr0) * munb"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 3,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
mds.min = .mds.min ,
multipleResponses = FALSE,
parameters.names = c("pstr0", "munb", "size"),
eps.trig = .eps.trig ,
type.fitted = .type.fitted ,
nsimEIM = .nsimEIM ,
zero = .zero )
}, list( .zero = zero,
.nsimEIM = nsimEIM, .eps.trig = eps.trig,
.type.fitted = type.fitted,
.mds.min = mds.min))),
rqresslot = eval(substitute(
function(mu, y, w, eta, extra = NULL) {
pstr0 <- eta2theta(eta[, c(TRUE, FALSE, FALSE)], .lpstr0 ,
earg = .epstr0 )
munb <- eta2theta(eta[, c(FALSE, TRUE, FALSE)], .lmunb ,
earg = .emunb )
kmat <- eta2theta(eta[, c(FALSE, FALSE, TRUE)], .lsize ,
earg = .esize )
scrambleseed <- runif(1) # To scramble the seed
qnorm(runif(length(y),
pzinegbin(y - 1, size = kmat, munb = munb,
pstr0 = pstr0),
pzinegbin(y , size = kmat, munb = munb,
pstr0 = pstr0)))
},
list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize ))),
initialize = eval(substitute(expression({
M1 <- 3
temp16 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp16$w
y <- temp16$y
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("pstr0", NOS, skip1 = TRUE)
mynames2 <- param.names("munb", NOS, skip1 = TRUE)
mynames3 <- param.names("size", NOS, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lpstr0 , earg = .epstr0 , tag = FALSE),
namesof(mynames2, .lmunb , earg = .emunb , tag = FALSE),
namesof(mynames3, .lsize , earg = .esize , tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
gprobs.y <- .gprobs.y
imunb <- .imunb # Default in NULL
if (length(imunb))
imunb <- matrix(imunb, n, NOS, byrow = TRUE)
if (!length(etastart)) {
munb.init <-
size.init <- matrix(NA_real_, n, NOS)
gprobs.y <- .gprobs.y
if (length( .iprobs.y ))
gprobs.y <- .iprobs.y
gsize.mux <- .gsize.mux # gsize.mux is on a relative scale
for (jay in 1:NOS) { # For each response 'y_jay'... do:
TFvec <- y[, jay] > 0 # Important to exclude the 0s
posyvec <- y[TFvec, jay]
munb.init.jay <- if ( .imethod == 1 ) {
quantile(posyvec, probs = gprobs.y) - 1/2 # + 1/16
} else {
weighted.mean(posyvec, w = w[TFvec, jay]) - 1/2
}
if (length(imunb))
munb.init.jay <- imunb[, jay]
gsize <- gsize.mux * 0.5 * (mean(munb.init.jay) +
weighted.mean(posyvec, w = w[TFvec, jay]))
if (length( .isize ))
gsize <- .isize # isize is on an absolute scale
try.this <-
grid.search2(munb.init.jay, gsize,
objfun = posNBD.Loglikfun2,
y = posyvec, # x = x[TFvec, , drop = FALSE],
w = w[TFvec, jay],
ret.objfun = TRUE) # Last value is the loglik
munb.init[, jay] <- try.this["Value1"]
size.init[, jay] <- try.this["Value2"]
} # for (jay ...)
if (length( .ipstr0 )) {
pstr0.init <- matrix( .ipstr0 , n, ncoly, byrow = TRUE)
} else {
pstr0.init <- matrix(0, n, ncoly)
ipstr0.small <- .ipstr0.small # Easily represented
for (jay in 1:NOS) {
Phi.init <-
pmax(ipstr0.small,
weighted.mean(y[, jay] == 0, w[, jay]) -
dnbinom(0, mu = munb.init[, jay],
size = size.init[, jay]))
if (mean(Phi.init == ipstr0.small) > 0.95 &&
.lpstr0 != "identitylink")
warning("from the initial values only, the data appears ",
"to have little or no 0-inflation, and possibly ",
"0-deflation.")
pstr0.init[, jay] <- Phi.init
} # for (jay)
}
etastart <-
cbind(theta2eta(pstr0.init, .lpstr0 , earg = .epstr0 ),
theta2eta(munb.init, .lmunb , earg = .emunb ),
theta2eta(size.init, .lsize , earg = .esize ))
etastart <-
etastart[, interleave.VGAM(ncol(etastart), M1 = M1)]
}
}), list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize,
.ipstr0 = ipstr0, .imunb = imunb, .isize = isize,
.gprobs.y = gprobs.y, .gsize.mux = gsize.mux,
.type.fitted = type.fitted,
.iprobs.y = iprobs.y,
.ipstr0.small = ipstr0.small,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 3)
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0",
"pstr0", "onempstr0"))[1]
pstr0 <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
.lpstr0 , earg = .epstr0 )
if (type.fitted %in% c("mean", "munb", "pobs0"))
munb <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
.lmunb , earg = .emunb )
if (type.fitted %in% c("pobs0")) {
kmat <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
.lsize , earg = .esize )
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
prob0 <- tempk^kmat # p(0) from negative binomial
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
prob0[big.size] <- exp(-munb[big.size])
}
}
ans <- switch(type.fitted,
"mean" = (1 - pstr0) * munb,
"munb" = munb,
"pobs0" = pstr0 + (1 - pstr0) * prob0, # P(Y=0)
"pstr0" = pstr0,
"onempstr0" = 1 - pstr0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lpstr0 = lpstr0, .lsize = lsize, .lmunb = lmunb,
.epstr0 = epstr0, .esize = esize, .emunb = emunb,
.mds.min = mds.min ))),
last = eval(substitute(expression({
misc$link <-
c(rep_len( .lpstr0 , NOS),
rep_len( .lmunb , NOS),
rep_len( .lsize , NOS))[interleave.VGAM(M1*NOS,
M1 = M1)]
temp.names <-
c(mynames1,
mynames2,
mynames3)[interleave.VGAM(M1*NOS, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M1*NOS)
names(misc$earg) <- temp.names
for (ii in 1:NOS) {
misc$earg[[M1*ii-2]] <- .epstr0
misc$earg[[M1*ii-1]] <- .emunb
misc$earg[[M1*ii ]] <- .esize
}
misc$ipstr0 <- .ipstr0
misc$isize <- .isize
misc$max.chunk.MB <- .max.chunk.MB
misc$cutoff.prob <- .cutoff.prob
misc$imethod <- .imethod
misc$nsimEIM <- .nsimEIM
misc$expected <- TRUE
misc$multipleResponses <- TRUE
}), list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize,
.ipstr0 = ipstr0, .isize = isize,
.nsimEIM = nsimEIM, .imethod = imethod,
.cutoff.prob = cutoff.prob,
.max.chunk.MB = max.chunk.MB
))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
pstr0 <- eta2theta(eta[, c(TRUE, FALSE, FALSE)], .lpstr0 ,
earg = .epstr0 )
munb <- eta2theta(eta[, c(FALSE, TRUE, FALSE)], .lmunb ,
earg = .emunb )
kmat <- eta2theta(eta[, c(FALSE, FALSE, TRUE)], .lsize ,
earg = .esize )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzinegbin(x = y, size = kmat, munb = munb,
pstr0 = pstr0, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize ))),
vfamily = c("zinegbinomial"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
pstr0 <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
.lpstr0 , earg = .epstr0 )
munb <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
.lsize , earg = .esize )
rzinegbin(nsim * length(munb),
size = kmat, munb = munb, pstr0 = pstr0)
}, list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
M1 <- 3
NOS <- ncol(eta) / M1
pstr0 <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lpstr0 , earg = .epstr0 )
munb <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lmunb , earg = .emunb )
size <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lsize , earg = .esize )
okay1 <- all(is.finite(munb)) && all(0 < munb) &&
all(is.finite(size)) && all(0 < size) &&
all(is.finite(pstr0)) && all(pstr0 < 1)
prob <- size / (size + munb)
prob0 <- prob^size
Prob0.check <- dnbinom(0, size = size, prob = prob)
deflat.limit <- -prob0 / (1 - prob0)
okay2.deflat <- TRUE
if (okay1 && !(okay2.deflat <- all(deflat.limit < pstr0)))
warning("parameter 'pstr0' is too negative even ",
"allowing for 0-deflation.")
smallval <- .mds.min # .munb.div.size
overdispersion <- if (okay1 && okay2.deflat)
all(smallval < munb / size) else FALSE
if (!overdispersion)
warning("parameter 'size' has very large values; ",
"try fitting a zero-inflated Poisson ",
"model instead.")
okay1 && okay2.deflat && overdispersion
}, list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize,
.mds.min = mds.min))),
deriv = eval(substitute(expression({
M1 <- 3
NOS <- ncol(eta) / M1
pstr0 <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lpstr0 , earg = .epstr0 )
munb <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lsize , earg = .esize )
dpstr0.deta <- dtheta.deta(pstr0, .lpstr0 , earg = .epstr0 )
dmunb.deta <- dtheta.deta(munb , .lmunb , earg = .emunb )
dsize.deta <- dtheta.deta(kmat , .lsize , earg = .esize )
dthetas.detas <-
(cbind(dpstr0.deta,
dmunb.deta,
dsize.deta))[, interleave.VGAM(M1*NOS, M1 = M1)]
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
if (FALSE)
warning("parameter 'size' has very large values; ",
"try fitting a zero-inflated Poisson ",
"model instead")
kmat[big.size] <- munb[big.size] / smallval
}
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
tempm <- munb / (kmat + munb)
prob0 <- tempk^kmat
oneminusf0 <- 1 - prob0
AA16 <- tempm + log(tempk)
df0.dmunb <- -tempk * prob0
df0.dkmat <- prob0 * AA16
df02.dmunb2 <- prob0 * tempk * (1 + 1/kmat) / (1 + munb/kmat)
df02.dkmat2 <- prob0 * ((tempm^2) / kmat + AA16^2)
df02.dkmat.dmunb <- -prob0 *
(tempm / kmat + AA16) / (1 + munb / kmat)
AA <- pobs0 <- cbind(pstr0 + (1 - pstr0) * prob0)
dl.dpstr0 <- -1 / (1 - pstr0)
dl.dmunb <- y / munb - (1 + y/kmat) / (1 + munb/kmat)
dl.dsize <- digamma(y + kmat) - digamma(kmat) -
(y - munb) / (munb + kmat) + log(tempk)
if (any(big.size)) {
dl.dsize[big.size] <- 1e-7 # A small number
}
for (spp. in 1:NOS) {
index0 <- (y[, spp.] == 0)
if (all(index0) || all(!index0))
stop("must have some 0s AND some positive ",
"counts in the data")
pstr0. <- pstr0[index0, spp.]
tempk. <- tempk[index0, spp.] # kmat. / (kmat. + munb.)
tempm. <- tempm[index0, spp.] # munb. / (kmat. + munb.)
prob0. <- prob0[index0, spp.] # tempk.^kmat.
df0.dmunb. <- df0.dmunb[index0, spp.] # -tempk.* prob0.
df0.dkmat. <- df0.dkmat[index0, spp.] # prob0. * above
denom. <- AA[index0, spp.] # pstr0. + (1-pstr0.)*prob0.
dl.dpstr0[index0, spp.] <- (1 - prob0.) / denom.
dl.dmunb[index0, spp.] <-
(1 - pstr0.) * df0.dmunb. / denom.
dl.dsize[index0, spp.] <-
(1 - pstr0.) * df0.dkmat. / denom.
} # of spp.
dl.dthetas <-
cbind(dl.dpstr0,
dl.dmunb,
dl.dsize)[, interleave.VGAM(M1*NOS, M1 = M1)]
ans <- c(w) * dl.dthetas * dthetas.detas
ans
}), list( .lpstr0 = lpstr0, .lmunb = lmunb, .lsize = lsize,
.epstr0 = epstr0, .emunb = emunb, .esize = esize,
.mds.min = mds.min ))),
weight = eval(substitute(expression({
wz <- matrix(0, n, M + M-1 + M-2)
mymu <- munb / oneminusf0 # Is the same as 'mu', == E(Y)
max.support <- .max.support
max.chunk.MB <- .max.chunk.MB
ind2 <- matrix(FALSE, n, NOS) # Used for SFS
for (jay in 1:NOS) {
eff.p <- sort(c( .cutoff.prob , 1 - .cutoff.prob ))
Q.mins <- 1
Q.maxs <- qgaitdnbinom(p = eff.p[2], truncate = 0,
kmat[, jay],
munb.p = munb[, jay]) + 10
eps.trig <- .eps.trig
Q.MAXS <- pmax(10, ceiling(1 / sqrt(eps.trig)))
Q.maxs <- pmin(Q.maxs, Q.MAXS)
ind1 <- if (max.chunk.MB > 0)
(Q.maxs - Q.mins < max.support) else FALSE
if ((NN <- sum(ind1)) > 0) {
Object.Size <- NN * 8 * max(Q.maxs - Q.mins) / (2^20)
n.chunks <- if (intercept.only) 1 else
max(1, ceiling( Object.Size / max.chunk.MB))
chunk.rows <- ceiling(NN / n.chunks)
ind2[, jay] <- ind1 # Save this
wind2 <- which(ind1)
upr.ptr <- 0
lwr.ptr <- upr.ptr + 1
while (lwr.ptr <= NN) {
upr.ptr <- min(upr.ptr + chunk.rows, NN)
sind2 <- wind2[lwr.ptr:upr.ptr]
wz[sind2, M1*jay] <-
EIM.posNB.specialp(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only,
second.deriv = FALSE)
if (FALSE)
wz2[sind2, M1*jay] <-
EIM.posNB.speciald(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.min = min(Q.mins2[sind2]),
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only,
second.deriv = FALSE)
wz[sind2, M1*jay] <-
wz[sind2, M1*jay] * (1 - AA[sind2, jay]) -
(1-pstr0[sind2, jay]) * (df02.dkmat2[sind2, jay] -
(1-pstr0[sind2, jay]) *
(df0.dkmat[sind2, jay]^2) / AA[sind2, jay])
if (any(eim.kk.TF <- wz[sind2, M1*jay] <= 0 |
is.na(wz[sind2, M1*jay]))) {
ind2[sind2[eim.kk.TF], jay] <- FALSE
}
lwr.ptr <- upr.ptr + 1
} # while
}
} # end of for (jay in 1:NOS)
for (jay in 1:NOS) {
run.varcov <- 0
ii.TF <- !ind2[, jay] # Not assigned above
if (any(ii.TF)) {
kkvec <- kmat[ii.TF, jay]
muvec <- munb[ii.TF, jay]
PSTR0 <- pstr0[ii.TF, jay]
for (ii in 1:( .nsimEIM )) {
ysim <- rzinegbin(sum(ii.TF), pstr0 = PSTR0,
mu = muvec, size = kkvec)
index0 <- (ysim == 0)
dl.dk <- digamma(ysim + kkvec) - digamma(kkvec) -
(ysim - muvec) / (muvec + kkvec) +
log1p(-muvec / (kkvec + muvec)) # +
ans0 <- (1 - PSTR0) *
df0.dkmat[ii.TF , jay] / AA[ii.TF , jay]
dl.dk[index0] <- ans0[index0]
run.varcov <- run.varcov + dl.dk^2
} # end of for loop
run.varcov <- c(run.varcov / .nsimEIM )
ned2l.dk2 <- if (intercept.only)
mean(run.varcov) else run.varcov
wz[ii.TF, M1*jay] <- ned2l.dk2 # * (dsize.deta[ii.TF, jay])^2
}
}
wz[, M1*(1:NOS) ] <- wz[, M1*(1:NOS) ] * dsize.deta^2
save.weights <- !all(ind2)
ned2l.dpstr02 <- oneminusf0 / (AA * (1 - pstr0))
wz[, M1*(1:NOS) - 2] <- ned2l.dpstr02 * dpstr0.deta^2
ned2l.dpstr0.dmunb <- df0.dmunb / AA
wz[, M + M1*(1:NOS) - 2] <- ned2l.dpstr0.dmunb *
dpstr0.deta * dmunb.deta
ned2l.dpstr0.dsize <- df0.dkmat / AA
wz[, M + M-1 + M1*(1:NOS) - 2] <- ned2l.dpstr0.dsize *
dpstr0.deta * dsize.deta
ned2l.dmunb2 <-
(1 - AA) * (mymu / munb^2 -
((1 + mymu/kmat) / kmat) / (1 + munb/kmat)^2) -
(1-pstr0) * (df02.dmunb2 -
(1 - pstr0) * (df0.dmunb^2) / AA)
wz[, M1*(1:NOS) - 1] <- ned2l.dmunb2 * dmunb.deta^2
dAA.dmunb <- (1 - pstr0) * df0.dmunb
ned2l.dmunbsize <-
(1 - AA) * (munb - mymu) / (munb + kmat)^2 -
(1-pstr0) * (df02.dkmat.dmunb -
df0.dkmat * dAA.dmunb / AA)
wz[, M + M1*(1:NOS) - 1] <- ned2l.dmunbsize * dmunb.deta *
dsize.deta
w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = NOS)
}), list( .lpstr0 = lpstr0,
.epstr0 = epstr0, .nsimEIM = nsimEIM,
.cutoff.prob = cutoff.prob, .eps.trig = eps.trig,
.max.support = max.support,
.max.chunk.MB = max.chunk.MB ))))
} # End of zinegbinomial
zinegbinomialff.control <-
function(save.weights = TRUE, ...) {
list(save.weights = save.weights)
}
zinegbinomialff <-
function(lmunb = "loglink", lsize = "loglink",
lonempstr0 = "logitlink",
type.fitted = c("mean", "munb", "pobs0",
"pstr0", "onempstr0"),
imunb = NULL, isize = NULL, ionempstr0 = NULL,
zero = c("size", "onempstr0"),
imethod = 1,
iprobs.y = NULL, # 0.35,
cutoff.prob = 0.999, # higher is better for large 'size'
eps.trig = 1e-7,
max.support = 4000, # 20160127; I have changed this
max.chunk.MB = 30, # max.memory = Inf is allowed
gprobs.y = (0:9)/10, # 20160709; grid for finding munb.init
gsize.mux = exp((-12:6)/2),
mds.min = 1e-3,
nsimEIM = 500) {
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 2)
stop("argument 'imethod' must be 1 or 2")
lmunb <- as.list(substitute(lmunb))
emunb <- link2list(lmunb)
lmunb <- attr(emunb, "function.name")
lsize <- as.list(substitute(lsize))
esize <- link2list(lsize)
lsize <- attr(esize, "function.name")
lonempstr0 <- as.list(substitute(lonempstr0))
eonempstr0 <- link2list(lonempstr0)
lonempstr0 <- attr(eonempstr0, "function.name")
ipstr0.small <- 1/64 # A number easily represented exactly
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0",
"pstr0", "onempstr0"))[1]
if (!is.Numeric(eps.trig, length.arg = 1,
positive = TRUE) || eps.trig > 0.001)
stop("argument 'eps.trig' must be positive and ",
"smaller in value")
if (length(ionempstr0) &&
(!is.Numeric(ionempstr0, positive = TRUE) ||
any(ionempstr0 >= 1)))
stop("argument 'ionempstr0' must contain values in (0,1)")
if (length(isize) && !is.Numeric(isize, positive = TRUE))
stop("argument 'isize' must contain positive values only")
if (!is.Numeric(nsimEIM, length.arg = 1,
integer.valued = TRUE))
stop("argument 'nsimEIM' must be a positive integer")
if (nsimEIM <= 50)
warning("argument 'nsimEIM' should be greater ",
"than 50, say")
new("vglmff",
blurb = c("Zero-inflated negative binomial\n\n",
"Links: ",
namesof("munb", lmunb, emunb, tag = FALSE), ", ",
namesof("size", lsize, esize, tag = FALSE), ", ",
namesof("onempstr0", lonempstr0, eonempstr0,
tag = FALSE), "\n",
"Mean: (1 - pstr0) * munb"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 3,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 3,
Q1 = 1,
expected = TRUE,
mds.min = .mds.min ,
multipleResponses = TRUE,
parameters.names = c("munb", "size", "onempstr0"),
eps.trig = .eps.trig ,
nsimEIM = .nsimEIM ,
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.nsimEIM = nsimEIM, .eps.trig = eps.trig,
.type.fitted = type.fitted,
.mds.min = mds.min
))),
rqresslot = eval(substitute(
function(mu, y, w, eta, extra = NULL) {
M1 <- 3
NOS <- extra$NOS
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
onempstr0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempstr0 , earg = .eonempstr0 )
scrambleseed <- runif(1) # To scramble the seed
qnorm(runif(length(y),
pzinegbin(y - 1, size = kmat, munb = munb,
pstr0 = 1 - onempstr0),
pzinegbin(y , size = kmat, munb = munb,
pstr0 = 1 - onempstr0)))
},
list( .lonempstr0 = lonempstr0, .lsize = lsize, .lmunb = lmunb,
.eonempstr0 = eonempstr0, .esize = esize, .emunb = emunb
))),
initialize = eval(substitute(expression({
M1 <- 3
temp16 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp16$w
y <- temp16$y
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("munb", NOS, skip1 = TRUE)
mynames2 <- param.names("size", NOS, skip1 = TRUE)
mynames3 <- param.names("onempstr0", NOS, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lmunb , earg = .emunb , tag = FALSE),
namesof(mynames2, .lsize , earg = .esize , tag = FALSE),
namesof(mynames3, .lonempstr0 , earg = .eonempstr0 ,
tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
gprobs.y <- .gprobs.y
imunb <- .imunb # Default in NULL
if (length(imunb))
imunb <- matrix(imunb, n, NOS, byrow = TRUE)
if (!length(etastart)) {
munb.init <-
size.init <- matrix(NA_real_, n, NOS)
gprobs.y <- .gprobs.y
if (length( .iprobs.y ))
gprobs.y <- .iprobs.y
gsize.mux <- .gsize.mux # gsize.mux is on a relative scale
for (jay in 1:NOS) { # For each response 'y_jay'... do:
TFvec <- y[, jay] > 0 # Important to exclude the 0s
posyvec <- y[TFvec, jay]
munb.init.jay <- if ( .imethod == 1 ) {
quantile(posyvec, probs = gprobs.y) - 1/2 # + 1/16
} else {
weighted.mean(posyvec, w = w[TFvec, jay]) - 1/2
}
if (length(imunb))
munb.init.jay <- imunb[, jay]
gsize <- gsize.mux * 0.5 * (mean(munb.init.jay) +
weighted.mean(posyvec, w = w[TFvec, jay]))
if (length( .isize ))
gsize <- .isize # isize is on an absolute scale
try.this <-
grid.search2(munb.init.jay, gsize,
objfun = posNBD.Loglikfun2,
y = posyvec, # x = x[TFvec, , drop = FALSE],
w = w[TFvec, jay],
ret.objfun = TRUE) # Last value is the loglik
munb.init[, jay] <- try.this["Value1"]
size.init[, jay] <- try.this["Value2"]
} # for (jay ...)
if (length( .ionempstr0 )) {
onempstr0.init <- matrix( .ionempstr0 , n,
ncoly, byrow = TRUE)
} else {
onempstr0.init <- matrix(0, n, ncoly)
ipstr0.small <- .ipstr0.small # Easily represented exactly
for (jay in 1:NOS) {
Phi.init <- pmax(ipstr0.small,
weighted.mean(y[, jay] == 0, w[, jay]) -
dnbinom(0, mu = munb.init[, jay],
size = size.init[, jay]))
if (mean(Phi.init == ipstr0.small) > 0.95)
warning("from the initial values only, ",
"the data appears ",
"to have little or no 0-inflation")
onempstr0.init[, jay] <- 1 - Phi.init
} # for (jay)
}
etastart <-
cbind(theta2eta(munb.init, .lmunb , .emunb ),
theta2eta(size.init, .lsize , .esize ),
theta2eta(onempstr0.init, .lonempstr0 ,
earg = .eonempstr0 ))
etastart <-
etastart[, interleave.VGAM(ncol(etastart), M1 = M1)]
}
}),
list( .lonempstr0 = lonempstr0, .lmunb = lmunb, .lsize = lsize,
.eonempstr0 = eonempstr0, .emunb = emunb, .esize = esize,
.ionempstr0 = ionempstr0, .imunb = imunb, .isize = isize,
.gprobs.y = gprobs.y, .gsize.mux = gsize.mux,
.type.fitted = type.fitted,
.ipstr0.small = ipstr0.small,
.iprobs.y = iprobs.y,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <-
if (length(extra$type.fitted)) extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "munb", "pobs0",
"pstr0", "onempstr0"))[1]
M1 <- 3
NOS <- ncol(eta) / M1
if (type.fitted %in% c("mean", "munb", "pobs0"))
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
if (type.fitted %in% c("pobs0")) {
kmat <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
prob0 <- tempk^kmat # p(0) from negative binomial
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
prob0[big.size] <- exp(-munb[big.size])
}
}
onempstr0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempstr0 , earg = .eonempstr0 )
ans <- switch(type.fitted,
"mean" = onempstr0 * munb,
"munb" = munb,
"pobs0" = 1 - onempstr0 + onempstr0 * prob0, # P(Y=0)
"pstr0" = 1 - onempstr0,
"onempstr0" = onempstr0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
},
list( .lonempstr0 = lonempstr0, .lsize = lsize, .lmunb = lmunb,
.eonempstr0 = eonempstr0, .esize = esize, .emunb = emunb,
.type.fitted = type.fitted, .mds.min = mds.min ))),
last = eval(substitute(expression({
misc$link <-
c(rep_len( .lmunb , NOS),
rep_len( .lsize , NOS),
rep_len( .lonempstr0 , NOS))[interleave.VGAM(M1*NOS,
M1 = M1)]
temp.names <-
c(mynames1,
mynames2,
mynames3)[interleave.VGAM(M1*NOS, M1 = M1)]
names(misc$link) <- temp.names
misc$earg <- vector("list", M1*NOS)
names(misc$earg) <- temp.names
for (ii in 1:NOS) {
misc$earg[[M1*ii-2]] <- .emunb
misc$earg[[M1*ii-1]] <- .esize
misc$earg[[M1*ii ]] <- .eonempstr0
}
misc$imethod <- .imethod
misc$nsimEIM <- .nsimEIM
misc$expected <- TRUE
misc$M1 <- M1
misc$ionempstr0 <- .ionempstr0
misc$isize <- .isize
misc$multipleResponses <- TRUE
}),
list( .lonempstr0 = lonempstr0, .lmunb = lmunb, .lsize = lsize,
.eonempstr0 = eonempstr0, .emunb = emunb, .esize = esize,
.ionempstr0 = ionempstr0, .isize = isize,
.nsimEIM = nsimEIM, .imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
M1 <- 3
NOS <- extra$NOS
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
onempstr0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempstr0 , earg = .eonempstr0 )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzinegbin(x = y, size = kmat, munb = munb,
pstr0 = 1 - onempstr0, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
},
list( .lonempstr0 = lonempstr0, .lmunb = lmunb, .lsize = lsize,
.eonempstr0 = eonempstr0, .emunb = emunb, .esize = esize
))),
vfamily = c("zinegbinomialff"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
munb <- eta2theta(eta[, c(TRUE, FALSE, FALSE)],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, c(FALSE, TRUE, FALSE)],
.lsize , earg = .esize )
onempstr0 <- eta2theta(eta[, c(FALSE, FALSE, TRUE)],
.lpstr0 , earg = .epstr0 )
rzinegbin(nsim * length(munb),
size = kmat, munb = munb, pstr0 = 1 - onempstr0)
},
list( .lonempstr0 = lonempstr0, .lmunb = lmunb, .lsize = lsize,
.eonempstr0 = eonempstr0, .emunb = emunb, .esize = esize
))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
M1 <- 3
NOS <- ncol(eta) / M1
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
size <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
onempstr0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempstr0 , earg = .eonempstr0 )
okay1 <- all(is.finite(munb)) && all(0 < munb) &&
all(is.finite(size)) && all(0 < size) &&
all(is.finite(onempstr0)) && all(0 < onempstr0)
prob <- size / (size + munb)
prob0 <- prob^size
Prob0.check <- dnbinom(0, size = size, prob = prob)
deflat.limit <- -prob0 / (1 - prob0)
okay2.deflat <- TRUE
if (okay1 &&
!(okay2.deflat <- all(onempstr0 < 1 - deflat.limit)))
warning("parameter 'pstr0' is too positive even ",
"allowing for 0-deflation.")
smallval <- .mds.min # .munb.div.size
overdispersion <- if (okay1 && okay2.deflat)
all(smallval < munb / size) else FALSE
if (!overdispersion)
warning("parameter 'size' has very large values; ",
"try fitting a zero-inflated Poisson ",
"model instead.")
okay1 && okay2.deflat && overdispersion
},
list( .lonempstr0 = lonempstr0, .lmunb = lmunb, .lsize = lsize,
.eonempstr0 = eonempstr0, .emunb = emunb, .esize = esize,
.mds.min = mds.min))),
deriv = eval(substitute(expression({
M1 <- 3
NOS <- ncol(eta) / M1
munb <- eta2theta(eta[, M1*(1:NOS)-2, drop = FALSE],
.lmunb , earg = .emunb )
kmat <- eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lsize , earg = .esize )
onempstr0 <- eta2theta(eta[, M1*(1:NOS) , drop = FALSE],
.lonempstr0 , earg = .eonempstr0 )
donempstr0.deta <- dtheta.deta(onempstr0, .lonempstr0 ,
earg = .eonempstr0 )
dmunb.deta <- dtheta.deta(munb , .lmunb , earg = .emunb )
dsize.deta <- dtheta.deta(kmat , .lsize , earg = .esize )
dthetas.detas <-
(cbind(dmunb.deta,
dsize.deta,
donempstr0.deta))[, interleave.VGAM(M1*NOS,
M1 = M1)]
smallval <- .mds.min # Something like this is needed
if (any(big.size <- munb / kmat < smallval)) {
if (FALSE)
warning("parameter 'size' has very large values; ",
"try fitting a zero-inflated Poisson ",
"model instead")
kmat[big.size] <- munb[big.size] / smallval
}
tempk <- 1 / (1 + munb / kmat) # kmat / (kmat + munb)
tempm <- munb / (kmat + munb)
prob0 <- tempk^kmat
oneminusf0 <- 1 - prob0
AA16 <- tempm + log(tempk)
df0.dmunb <- -tempk * prob0
df0.dkmat <- cbind(prob0 * AA16)
df02.dmunb2 <- prob0 * tempk * (1 +
1/kmat) / (1 + munb/kmat)
df02.dkmat2 <- prob0 * ((tempm^2) / kmat + AA16^2)
df02.dkmat.dmunb <- -prob0 * (tempm / kmat +
AA16) / (1 + munb/kmat)
pstr0 <- 1 - onempstr0
AA <- pobs0 <- cbind(pstr0 + (onempstr0) * prob0)
dl.dmunb <- y / munb - (1 + y/kmat) / (1 + munb/kmat)
dl.dsize <- digamma(y + kmat) - digamma(kmat) -
(y - munb) / (munb + kmat) + log(tempk)
dl.donempstr0 <- +1 / (onempstr0)
for (spp. in 1:NOS) {
index0 <- (y[, spp.] == 0)
if (all(index0) || all(!index0))
stop("must have some 0s AND some positive ",
"counts in the data")
kmat. <- kmat[index0, spp.]
munb. <- munb[index0, spp.]
onempstr0. <- onempstr0[index0, spp.]
tempk. <- kmat. / (kmat. + munb.)
tempm. <- munb. / (kmat. + munb.)
prob0. <- tempk.^kmat.
df0.dmunb. <- -tempk.* prob0.
df0.dkmat. <- prob0. * (tempm. + log(tempk.))
denom. <- 1 - onempstr0. + (onempstr0.) * prob0.
dl.donempstr0[index0, spp.] <-
-(1 - prob0.) / denom. # note "-"
dl.dmunb[index0, spp.] <- (onempstr0.) * df0.dmunb. / denom.
dl.dsize[index0, spp.] <- (onempstr0.) * df0.dkmat. / denom.
} # of spp.
dl.dthetas <-
cbind(dl.dmunb,
dl.dsize,
dl.donempstr0)[, interleave.VGAM(M1*NOS, M1 = M1)]
c(w) * dl.dthetas * dthetas.detas
}),
list( .lonempstr0 = lonempstr0, .lmunb = lmunb, .lsize = lsize,
.eonempstr0 = eonempstr0, .emunb = emunb, .esize = esize,
.mds.min = mds.min ))),
weight = eval(substitute(expression({
wz <- matrix(0, n, M + M-1 + M-2)
mymu <- munb / oneminusf0 # Is the same as 'mu', == E(Y)
max.support <- .max.support
max.chunk.MB <- .max.chunk.MB
ind2 <- matrix(FALSE, n, NOS) # Used for SFS
for (jay in 1:NOS) {
eff.p <- sort(c( .cutoff.prob , 1 - .cutoff.prob ))
Q.mins <- 1
Q.maxs <- qgaitdnbinom(p = eff.p[2], truncate = 0,
kmat[, jay],
munb.p = munb[, jay]) + 10
eps.trig <- .eps.trig
Q.MAXS <- pmax(10, ceiling(1 / sqrt(eps.trig)))
Q.maxs <- pmin(Q.maxs, Q.MAXS)
ind1 <- if (max.chunk.MB > 0)
(Q.maxs - Q.mins < max.support) else FALSE
if ((NN <- sum(ind1)) > 0) {
Object.Size <- NN * 8 * max(Q.maxs - Q.mins) / (2^20)
n.chunks <- if (intercept.only) 1 else
max(1, ceiling( Object.Size / max.chunk.MB))
chunk.rows <- ceiling(NN / n.chunks)
ind2[, jay] <- ind1 # Save this
wind2 <- which(ind1)
upr.ptr <- 0
lwr.ptr <- upr.ptr + 1
while (lwr.ptr <= NN) {
upr.ptr <- min(upr.ptr + chunk.rows, NN)
sind2 <- wind2[lwr.ptr:upr.ptr]
wz[sind2, M1*jay - 1] <-
EIM.posNB.specialp(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only,
second.deriv = FALSE)
if (FALSE)
wz2[sind2, M1*jay - 1] <-
EIM.posNB.speciald(munb = munb[sind2, jay],
size = kmat[sind2, jay],
y.min = min(Q.mins2[sind2]),
y.max = max(Q.maxs[sind2]),
cutoff.prob = .cutoff.prob ,
prob0 = prob0[sind2, jay],
df0.dkmat = df0.dkmat[sind2, jay],
df02.dkmat2 = df02.dkmat2[sind2, jay],
intercept.only = intercept.only,
second.deriv = FALSE)
wz[sind2, M1*jay - 1] <-
wz[sind2, M1*jay - 1] * (1 - AA[sind2, jay]) -
(1-pstr0[sind2, jay]) * (df02.dkmat2[sind2, jay] -
(1-pstr0[sind2, jay]) *
(df0.dkmat[sind2, jay]^2) / AA[sind2, jay])
if (any(eim.kk.TF <- wz[sind2, M1*jay - 1] <= 0 |
is.na(wz[sind2, M1*jay - 1]))) {
ind2[sind2[eim.kk.TF], jay] <- FALSE
}
lwr.ptr <- upr.ptr + 1
} # while
}
} # end of for (jay in 1:NOS)
for (jay in 1:NOS) {
run.varcov <- 0
ii.TF <- !ind2[, jay] # Not assigned above
if (any(ii.TF)) {
kkvec <- kmat[ii.TF, jay]
muvec <- munb[ii.TF, jay]
PSTR0 <- pstr0[ii.TF, jay]
for (ii in 1:( .nsimEIM )) {
ysim <- rzinegbin(sum(ii.TF), pstr0 = PSTR0,
mu = muvec, size = kkvec)
index0 <- (ysim == 0)
dl.dk <- digamma(ysim + kkvec) - digamma(kkvec) -
(ysim - muvec) / (muvec + kkvec) +
log1p(-muvec / (kkvec + muvec)) # +
ans0 <- (1 - PSTR0) *
df0.dkmat[ii.TF , jay] / AA[ii.TF , jay]
dl.dk[index0] <- ans0[index0]
run.varcov <- run.varcov + dl.dk^2
} # end of for loop
run.varcov <- c(run.varcov / .nsimEIM )
ned2l.dk2 <- if (intercept.only)
mean(run.varcov) else run.varcov
wz[ii.TF, M1*jay - 1] <- ned2l.dk2 # *
}
}
wz[, M1*(1:NOS) - 1] <- wz[, M1*(1:NOS) - 1] * dsize.deta^2
save.weights <- !all(ind2)
ned2l.donempstr02 <- oneminusf0 / (AA * (onempstr0))
wz[, M1*(1:NOS) ] <- ned2l.donempstr02 *
donempstr0.deta^2
ned2l.donempstr0.dmunb <- -df0.dmunb / AA # Negated (1/2)
wz[, M + M-1 + M1*(1:NOS) - 2] <- ned2l.donempstr0.dmunb *
donempstr0.deta *
dmunb.deta
ned2l.donempstr0.dsize <- -df0.dkmat / AA # Negated (2/2)
wz[, M + M1*(1:NOS) - 1] <- ned2l.donempstr0.dsize *
donempstr0.deta *
dsize.deta
ned2l.dmunb2 <-
(1 - AA) * (mymu / munb^2 -
((1 + mymu/kmat) / kmat) / (1 + munb/kmat)^2) -
(1-pstr0) * (df02.dmunb2 -
(1 - pstr0) * (df0.dmunb^2) / AA)
wz[, M1*(1:NOS) - 2] <- ned2l.dmunb2 * dmunb.deta^2
dAA.dmunb <- (onempstr0) * df0.dmunb
ned2l.dmunbsize <-
(1 - AA) * (munb - mymu) / (munb + kmat)^2 -
(onempstr0) * (df02.dkmat.dmunb -
df0.dkmat * dAA.dmunb / AA)
wz[, M + M1*(1:NOS) - 2] <- ned2l.dmunbsize *
dmunb.deta *
dsize.deta
w.wz.merge(w = w, wz = wz, n = n, M = M, ndepy = NOS)
}), list( .lonempstr0 = lonempstr0,
.eonempstr0 = eonempstr0, .nsimEIM = nsimEIM,
.cutoff.prob = cutoff.prob, .eps.trig = eps.trig,
.max.support = max.support,
.max.chunk.MB = max.chunk.MB ))))
} # End of zinegbinomialff
dzigeom <- function(x, prob, pstr0 = 0, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(prob), length(pstr0))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- dgeom(x = x, prob = prob, log = TRUE)
ans <- if (log.arg) {
ifelse(x == 0, log(pstr0 + (1 - pstr0) * exp(ans)),
log1p(-pstr0) + ans)
} else {
ifelse(x == 0, pstr0 + (1 - pstr0) * exp(ans) ,
(1 - pstr0) * exp(ans))
}
prob0 <- prob
deflat.limit <- -prob0 / (1 - prob0)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
pzigeom <- function(q, prob, pstr0 = 0) {
LLL <- max(length(q), length(prob), length(pstr0))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pstr0) != LLL) pstr0 <- rep_len(pstr0, LLL)
ans <- pgeom(q, prob)
ans <- ifelse(q < 0, 0, pstr0 + (1-pstr0) * ans)
prob0 <- prob
deflat.limit <- -prob0 / (1 - prob0)
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
qzigeom <- function(p, prob, pstr0 = 0) {
LLL <- max(length(p), length(prob), length(pstr0))
ans <- p <- rep_len(p, LLL)
prob <- rep_len(prob, LLL)
pstr0 <- rep_len(pstr0, LLL)
ans[p <= pstr0] <- 0
ind1 <- (p > pstr0)
ans[ind1] <-
qgeom((p[ind1] - pstr0[ind1]) / (1 - pstr0[ind1]),
prob = prob[ind1])
prob0 <- prob
deflat.limit <- -prob0 / (1 - prob0)
ind0 <- (deflat.limit <= pstr0) & (pstr0 < 0)
if (any(ind0)) {
pobs0 <- pstr0[ind0] + (1 - pstr0[ind0]) * prob0[ind0]
ans[p[ind0] <= pobs0] <- 0
pindex <- (1:LLL)[ind0 & (p > pobs0)]
Pobs0 <- pstr0[pindex] + (1 - pstr0[pindex]) * prob0[pindex]
ans[pindex] <- 1 + qgeom((p[pindex] - Pobs0) / (1 - Pobs0),
prob = prob[pindex])
}
ans[pstr0 < deflat.limit] <- NaN
ans[pstr0 > 1] <- NaN
ans
}
rzigeom <- function(n, prob, pstr0 = 0) {
qzigeom(runif(n), prob, pstr0 = pstr0)
}
zigeometric <-
function(
lpstr0 = "logitlink",
lprob = "logitlink",
type.fitted = c("mean", "prob", "pobs0",
"pstr0", "onempstr0"),
ipstr0 = NULL, iprob = NULL,
imethod = 1,
bias.red = 0.5,
zero = NULL) {
expected <- TRUE
lpstr0 <- as.list(substitute(lpstr0))
epstr0 <- link2list(lpstr0)
lpstr0 <- attr(epstr0, "function.name")
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "pstr0", "onempstr0"))[1]
if (length(ipstr0))
if (!is.Numeric(ipstr0, positive = TRUE) ||
ipstr0 >= 1)
stop("argument 'ipstr0' is out of range")
if (length(iprob))
if (!is.Numeric(iprob, positive = TRUE) ||
iprob >= 1)
stop("argument 'iprob' is out of range")
if (!is.Numeric(bias.red, length.arg = 1, positive = TRUE) ||
bias.red > 1)
stop("argument 'bias.red' must be between 0 and 1")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
new("vglmff",
blurb = c("Zero-inflated geometric distribution,\n",
"P[Y = 0] = pstr0 + (1 - pstr0) * prob,\n",
"P[Y = y] = (1 - pstr0) * prob * (1 - prob)^y, ",
"y = 1, 2, ...\n\n",
"Link: ",
namesof("pstr0", lpstr0, earg = epstr0), ", ",
namesof("prob", lprob, earg = eprob ), "\n",
"Mean: (1 - pstr0) * (1 - prob) / prob"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("pstr0", "prob"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted ))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("pstr0", ncoly, skip1 = TRUE)
mynames2 <- param.names("prob", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lpstr0, .epstr0, tag = FALSE),
namesof(mynames2, .lprob, .eprob, tag = FALSE))[
interleave.VGAM(M1 * NOS, M1 = M1)]
if (!length(etastart)) {
prob.init <- if ( .imethod == 3)
.bias.red / (1 + y + 1/8) else
if ( .imethod == 2)
.bias.red / (1 +
matrix(colMeans(y) + 1/8,
n, ncoly, byrow = TRUE)) else
.bias.red / (1 +
matrix(colSums(y * w) / colSums(w) + 1/8,
n, ncoly, byrow = TRUE))
prob.init <- if (length( .iprob )) {
matrix( .iprob , n, ncoly, byrow = TRUE)
} else {
prob.init # Already a matrix
}
prob0.est <- psze.init <- matrix(0, n, NOS)
for (jlocal in 1:NOS) {
prob0.est[, jlocal] <-
sum(w[y[, jlocal] == 0, jlocal]) / sum(w[, jlocal])
psze.init[, jlocal] <- if ( .imethod == 3)
prob0.est[, jlocal] / 2 else
if ( .imethod == 1)
pmax(0.05, (prob0.est[, jlocal] -
median(prob.init[, jlocal]))) else
prob0.est[, jlocal] / 5
}
psze.init <- if (length( .ipstr0 )) {
matrix( .ipstr0 , n, ncoly, byrow = TRUE)
} else {
psze.init # Already a matrix
}
etastart <-
cbind(theta2eta(psze.init, .lpstr0, earg = .epstr0),
theta2eta(prob.init, .lprob , earg = .eprob ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
}
}), list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0,
.iprob = iprob, .ipstr0 = ipstr0,
.type.fitted = type.fitted,
.bias.red = bias.red,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
pstr0 <- eta2theta(eta[, c(TRUE, FALSE)], .lpstr0 , .epstr0 )
prob <- eta2theta(eta[, c(FALSE, TRUE)], .lprob , .eprob )
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0",
"pstr0", "onempstr0"))[1]
ans <- switch(type.fitted,
"mean" = (1 - pstr0) * (1 - prob) / prob,
"prob" = prob,
"pobs0" = pstr0 + (1 - pstr0) * prob, # P(Y=0)
"pstr0" = pstr0,
"onempstr0" = 1 - pstr0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0 ))),
last = eval(substitute(expression({
temp.names <- c(rep_len( .lpstr0 , NOS),
rep_len( .lprob , NOS))
temp.names <- temp.names[interleave.VGAM(M1*NOS, M1 = M1)]
misc$link <- temp.names
misc$earg <- vector("list", M1 * NOS)
names(misc$link) <-
names(misc$earg) <-
c(mynames1, mynames2)[interleave.VGAM(M1*NOS, M1 = M1)]
for (ii in 1:NOS) {
misc$earg[[M1*ii-1]] <- .epstr0
misc$earg[[M1*ii ]] <- .eprob
}
misc$imethod <- .imethod
misc$zero <- .zero
misc$bias.red <- .bias.red
misc$expected <- .expected
misc$ipstr0 <- .ipstr0
misc$pobs0 <- pobs0
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pobs0) <- dimnames(y)
misc$pstr0 <- pstr0
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pstr0) <- dimnames(y)
}), list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0,
.ipstr0 = ipstr0,
.zero = zero,
.expected = expected,
.bias.red = bias.red,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
pstr0 <- eta2theta(eta[, c(TRUE, FALSE)], .lpstr0 , .epstr0 )
prob <- eta2theta(eta[, c(FALSE, TRUE)], .lprob , .eprob )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzigeom(y, prob = prob, pstr0 = pstr0, log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0 ))),
vfamily = c("zigeometric"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
pstr0 <- eta2theta(eta[, c(TRUE, FALSE)], .lpstr0 , .epstr0 )
prob <- eta2theta(eta[, c(FALSE, TRUE)], .lprob , .eprob )
rzigeom(nsim * length(pstr0), prob = prob, pstr0 = pstr0)
}, list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0 ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
pstr0 <- eta2theta(eta[, c(TRUE, FALSE)], .lpstr0 , .epstr0 )
prob <- eta2theta(eta[, c(FALSE, TRUE)], .lprob , .eprob )
okay1 <- all(is.finite(prob )) && all(0 < prob & prob < 1) &&
all(is.finite(pstr0)) && all(pstr0 < 1)
prob0 <- prob
deflat.limit <- -prob0 / (1 - prob0)
okay2.deflat <- TRUE
if (okay1 && !(okay2.deflat <- all(deflat.limit < pstr0)))
warning("parameter 'pstr0' is too negative even ",
"allowing for 0-deflation.")
okay1 && okay2.deflat
}, list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0 ))),
deriv = eval(substitute(expression({
M1 <- 2
pstr0 <- eta2theta(eta[, c(TRUE, FALSE)], .lpstr0 , .epstr0 )
prob <- eta2theta(eta[, c(FALSE, TRUE)], .lprob , .eprob )
prob0 <- prob # P(Y == 0) from parent distribution, aka f(0)
pobs0 <- pstr0 + (1 - pstr0) * prob0 # P(Y == 0)
index0 <- (y == 0)
dl.dpstr0 <- (1 - prob0) / pobs0
dl.dpstr0[!index0] <- -1 / (1 - pstr0[!index0])
dl.dprob <- (1 - pstr0) / pobs0
dl.dprob[!index0] <- 1 / prob[!index0] -
y[!index0] / (1 - prob[!index0])
dpstr0.deta <- dtheta.deta(pstr0 , .lpstr0 , earg = .epstr0 )
dprob.deta <- dtheta.deta(prob, .lprob , earg = .eprob )
dl.deta12 <- c(w) * cbind(dl.dpstr0 * dpstr0.deta,
dl.dprob * dprob.deta)
dl.deta12 <- dl.deta12[, interleave.VGAM(ncol(dl.deta12),
M1 = M1)]
dl.deta12
}), list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0 ))),
weight = eval(substitute(expression({
if ( .expected ) {
ned2l.dprob2 <- (1 - pstr0)^2 / pobs0 +
(1 - pstr0) * ((1 - prob) / prob) *
(1 / prob + 1 / (1 - prob)^2)
ned2l.dpstr0.prob <- 1 / pobs0
ned2l.dpstr02 <- (1 - prob0) / ((1 - pstr0) * pobs0)
} else {
od2l.dprob2 <- ((1 - pstr0) / pobs0)^2
od2l.dprob2[!index0] <- 1 / (prob[!index0])^2 +
y[!index0] / (1 - prob[!index0])^2
od2l.dpstr0.prob <- (pobs0 + (1 - prob0) * (1 -
pstr0)) / pobs0^2
od2l.dpstr0.prob[!index0] <- 0
od2l.dpstr02 <- ((1 - prob0) / pobs0)^2
od2l.dpstr02[!index0] <- 1 / (1 - pstr0[!index0])^2
}
allvals <- if ( .expected )
c(c(w) * ned2l.dpstr02 * dpstr0.deta^2,
c(w) * ned2l.dprob2 * dprob.deta^2,
c(w) * ned2l.dpstr0.prob * dprob.deta *
dpstr0.deta) else
c(c(w) * od2l.dpstr02 * dpstr0.deta^2,
c(w) * od2l.dprob2 * dprob.deta^2,
c(w) * od2l.dpstr0.prob * dprob.deta *
dpstr0.deta)
wz <- array(allvals, dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}), list( .lprob = lprob, .lpstr0 = lpstr0,
.eprob = eprob, .epstr0 = epstr0,
.expected = expected ))))
}
zigeometricff <-
function(lprob = "logitlink",
lonempstr0 = "logitlink",
type.fitted = c("mean", "prob", "pobs0",
"pstr0", "onempstr0"),
iprob = NULL, ionempstr0 = NULL,
imethod = 1,
bias.red = 0.5,
zero = "onempstr0") {
expected <- TRUE
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
lonempstr0 <- as.list(substitute(lonempstr0))
eonempstr0 <- link2list(lonempstr0)
lonempstr0 <- attr(eonempstr0, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "pstr0", "onempstr0"))[1]
if (length(iprob))
if (!is.Numeric(iprob, positive = TRUE) ||
iprob >= 1)
stop("argument 'iprob' is out of range")
if (length(ionempstr0))
if (!is.Numeric(ionempstr0, positive = TRUE) ||
ionempstr0 >= 1)
stop("argument 'ionempstr0' is out of range")
if (!is.Numeric(bias.red, length.arg = 1, positive = TRUE) ||
bias.red > 1)
stop("argument 'bias.red' must be between 0 and 1")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
new("vglmff",
blurb = c("Zero-inflated geometric distribution,\n",
"P[Y = 0] = 1 - onempstr0 + onempstr0 * prob,\n",
"P[Y = y] = onempstr0 * prob * (1 - prob)^y, ",
"y = 1, 2, ...\n\n",
"Link: ",
namesof("prob", lprob, eprob ), ", ",
namesof("onempstr0", lonempstr0, eonempstr0), "\n",
"Mean: onempstr0 * (1 - prob) / prob"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("prob", "onempstr0"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted ))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("prob", ncoly, skip1 = TRUE)
mynames2 <- param.names("onempstr0", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lprob , .eprob , tag = FALSE),
namesof(mynames2, .lonempstr0 , .eonempstr0 , tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
if (!length(etastart)) {
prob.init <- if ( .imethod == 3)
.bias.red / (1 + y + 1/8) else
if ( .imethod == 2)
.bias.red / (1 +
matrix(colMeans(y) + 1/8,
n, ncoly, byrow = TRUE)) else
.bias.red / (1 +
matrix(colSums(y * w) / colSums(w) + 1/8,
n, ncoly, byrow = TRUE))
prob.init <- if (length( .iprob )) {
matrix( .iprob , n, ncoly, byrow = TRUE)
} else {
prob.init # Already a matrix
}
prob0.est <- psze.init <- matrix(0, n, NOS)
for (jlocal in 1:NOS) {
prob0.est[, jlocal] <-
sum(w[y[, jlocal] == 0, jlocal]) / sum(w[, jlocal])
psze.init[, jlocal] <- if ( .imethod == 3)
prob0.est[, jlocal] / 2 else
if ( .imethod == 1)
pmax(0.05, (prob0.est[, jlocal] -
median(prob.init[, jlocal]))) else
prob0.est[, jlocal] / 5
}
psze.init <- if (length( .ionempstr0 )) {
matrix( 1 - .ionempstr0 , n, ncoly, byrow = TRUE)
} else {
psze.init # Already a matrix
}
etastart <-
cbind(theta2eta( prob.init, .lprob , .eprob ),
theta2eta(1 - psze.init, .lonempstr0 , .eonempstr0 ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
}
}), list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0,
.iprob = iprob, .ionempstr0 = ionempstr0,
.type.fitted = type.fitted,
.bias.red = bias.red,
.imethod = imethod ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
prob <- eta2theta(eta[, c(TRUE, FALSE)], .lprob ,
earg = .eprob )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0",
"pstr0", "onempstr0"))[1]
ans <- switch(type.fitted,
"mean" = onempstr0 * (1 - prob) / prob,
"prob" = prob,
"pobs0" = 1 - onempstr0 + onempstr0 * prob, # P(Y=0)
"pstr0" = 1 - onempstr0,
"onempstr0" = onempstr0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0 ))),
last = eval(substitute(expression({
temp.names <- c(rep_len( .lprob , NOS),
rep_len( .lonempstr0 , NOS))
temp.names <- temp.names[interleave.VGAM(M1*NOS, M1 = M1)]
misc$link <- temp.names
misc$earg <- vector("list", M1 * NOS)
names(misc$link) <-
names(misc$earg) <-
c(mynames1, mynames2)[interleave.VGAM(M1*NOS, M1 = M1)]
for (ii in 1:NOS) {
misc$earg[[M1*ii-1]] <- .eprob
misc$earg[[M1*ii ]] <- .eonempstr0
}
misc$imethod <- .imethod
misc$zero <- .zero
misc$bias.red <- .bias.red
misc$expected <- .expected
misc$ionempstr0 <- .ionempstr0
misc$pobs0 <- pobs0
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$pobs0) <- dimnames(y)
misc$onempstr0 <- onempstr0
if (length(dimnames(y)[[2]]) > 0)
dimnames(misc$onempstr0) <- dimnames(y)
}), list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0,
.ionempstr0 = ionempstr0,
.zero = zero,
.expected = expected,
.bias.red = bias.red,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
prob <- eta2theta(eta[, c(TRUE, FALSE)], .lprob ,
earg = .eprob )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzigeom(x = y, prob = prob, pstr0 = 1 - onempstr0,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0 ))),
vfamily = c("zigeometricff"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
prob <- eta2theta(eta[, c(TRUE, FALSE)], .lprob ,
earg = .eprob )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
rzigeom(nsim * length(onempstr0),
prob = prob, pstr0 = 1 - onempstr0)
}, list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0 ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
prob <- eta2theta(eta[, c(TRUE, FALSE)], .lprob ,
earg = .eprob )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
okay1 <- all(is.finite(onempobs0)) && all(0 < onempobs0) &&
all(is.finite(prob )) && all(0 < prob & prob < 1)
prob0 <- prob
deflat.limit <- -prob0 / (1 - prob0)
okay2.deflat <- TRUE
if (okay1 &&
!(okay2.deflat <- all(onempstr0 < 1 - deflat.limit)))
warning("parameter 'onempstr0' is too positive even ",
"allowing for 0-deflation.")
okay1 && okay2.deflat
}, list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0 ))),
deriv = eval(substitute(expression({
M1 <- 2
prob <- eta2theta(eta[, c(TRUE, FALSE)], .lprob ,
earg = .eprob )
onempstr0 <- eta2theta(eta[, c(FALSE, TRUE)], .lonempstr0 ,
earg = .eonempstr0 )
prob0 <- prob # P(Y == 0) from the parent distribution
pobs0 <- 1 - onempstr0 + (onempstr0) * prob0 # P(Y == 0)
index0 <- (y == 0)
dl.donempstr0 <- -(1 - prob0) / pobs0 # zz
dl.donempstr0[!index0] <- 1 / (onempstr0[!index0]) # zz
dl.dprob <- (onempstr0) / pobs0
dl.dprob[!index0] <- 1 / prob[!index0] -
y[!index0] / (1 - prob[!index0])
dprob.deta <- dtheta.deta(prob , .lprob ,
earg = .eprob )
donempstr0.deta <- dtheta.deta(onempstr0 , .lonempstr0 ,
earg = .eonempstr0 )
dl.deta12 <- c(w) * cbind(dl.dprob * dprob.deta,
dl.donempstr0 * donempstr0.deta)
dl.deta12 <- dl.deta12[, interleave.VGAM(ncol(dl.deta12),
M1 = M1)]
dl.deta12
}), list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0 ))),
weight = eval(substitute(expression({
if ( .expected ) {
ned2l.dprob2 <- (onempstr0)^2 / pobs0 +
(onempstr0) * ((1 - prob) / prob) *
(1 / prob + 1 / (1 - prob)^2)
ned2l.donempstr0.prob <- -1 / pobs0
ned2l.donempstr02 <- (1 - prob0) / ((onempstr0) * pobs0)
} else {
od2l.dprob2 <- (( onempstr0) / pobs0)^2
od2l.dprob2[!index0] <- 1 / (prob[!index0])^2 +
y[!index0] / (1 - prob[!index0])^2
od2l.donempstr0.prob <- -(pobs0 + (1 - prob0) *
(onempstr0)) / pobs0^2
od2l.donempstr0.prob[!index0] <- 0
od2l.donempstr02 <- ((1 - prob0) / pobs0)^2
od2l.donempstr02[!index0] <- 1 / ( onempstr0[!index0])^2
}
allvals <- if ( .expected )
c(c(w) * ned2l.dprob2 * dprob.deta^2,
c(w) * ned2l.donempstr02 * donempstr0.deta^2,
c(w) * ned2l.donempstr0.prob *
dprob.deta *
donempstr0.deta) else
c(c(w) * od2l.dprob2 * dprob.deta^2,
c(w) * od2l.donempstr02 * donempstr0.deta^2,
c(w) * od2l.donempstr0.prob *
dprob.deta *
donempstr0.deta)
wz <- array(allvals, dim = c(n, M / M1, 3))
wz <- arwz2wz(wz, M = M, M1 = M1)
wz
}),
list( .lprob = lprob, .lonempstr0 = lonempstr0,
.eprob = eprob, .eonempstr0 = eonempstr0,
.expected = expected ))))
}
dzageom <- function(x, prob, pobs0 = 0, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(prob), length(pobs0))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
index0 <- (x == 0)
if (log.arg) {
ans[ index0] <- log(pobs0[index0])
ans[!index0] <- log1p(-pobs0[!index0]) +
dposgeom(x[!index0],
prob = prob[!index0], log = TRUE)
} else {
ans[ index0] <- pobs0[index0]
ans[!index0] <- (1-pobs0[!index0]) *
dposgeom(x[!index0],
prob = prob[!index0])
}
ans
}
pzageom <- function(q, prob, pobs0 = 0) {
LLL <- max(length(q), length(prob), length(pobs0))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
ans[q > 0] <- pobs0[q > 0] +
(1 - pobs0[q > 0]) *
pposgeom(q[q > 0], prob = prob[q > 0])
ans[q < 0] <- 0
ans[q == 0] <- pobs0[q == 0]
ans <- pmax(0, ans)
ans <- pmin(1, ans)
ans
}
qzageom <- function(p, prob, pobs0 = 0) {
LLL <- max(length(p), length(prob), length(pobs0))
if (length(p) != LLL) p <- rep_len(p, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
ans <- p
ind4 <- (p > pobs0)
ans[!ind4] <- 0.0
ans[ ind4] <- qposgeom((p[ind4] - pobs0[ind4]) / (1 -
pobs0[ind4]),
prob = prob[ind4])
ans
}
rzageom <- function(n, prob, pobs0 = 0) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
ans <- rposgeom(use.n, prob)
if (length(pobs0) != use.n)
pobs0 <- rep_len(pobs0, use.n)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be between 0 and 1 inclusive")
ifelse(runif(use.n) < pobs0, 0, ans)
}
dzabinom <- function(x, size, prob, pobs0 = 0, log = FALSE) {
if (!is.logical(log.arg <- log) || length(log) != 1)
stop("bad input for argument 'log'")
rm(log)
LLL <- max(length(x), length(size), length(prob),
length(pobs0))
if (length(x) != LLL) x <- rep_len(x, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
index0 <- (x == 0)
if (log.arg) {
ans[ index0] <- log(pobs0[index0])
ans[!index0] <- log1p(-pobs0[!index0]) +
dgaitdbinom(x[!index0], size[!index0],
prob[!index0],
truncate = 0, log = TRUE)
} else {
ans[ index0] <- pobs0[index0]
ans[!index0] <- (1-pobs0[!index0]) *
dgaitdbinom(x[!index0], size[!index0],
prob[!index0],
truncate = 0)
}
ans
} # dzabinom
pzabinom <- function(q, size, prob, pobs0 = 0) {
LLL <- max(length(q), length(size), length(prob),
length(pobs0))
if (length(q) != LLL) q <- rep_len(q, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
ans <- rep_len(0.0, LLL)
if (!is.Numeric(pobs0) ||
any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
ans[q > 0] <- pobs0[q > 0] +
(1 - pobs0[q > 0]) *
pgaitdbinom(q[q > 0], size[q > 0], prob[q > 0],
truncate = 0)
ans[q < 0] <- 0
ans[q == 0] <- pobs0[q == 0]
ans <- pmax(0, ans)
ans <- pmin(1, ans)
ans
} # pzabinom
qzabinom <- function(p, size, prob, pobs0 = 0) {
LLL <- max(length(p), length(size), length(prob),
length(pobs0))
if (length(p) != LLL) p <- rep_len(p, LLL)
if (length(size) != LLL) size <- rep_len(size, LLL)
if (length(prob) != LLL) prob <- rep_len(prob, LLL)
if (length(pobs0) != LLL) pobs0 <- rep_len(pobs0, LLL)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be in [0,1]")
ans <- p
ind4 <- (p > pobs0)
ans[!ind4] <- 0.0
ans[ ind4] <- qgaitdbinom((p[ind4] - pobs0[ind4]) / (
1 - pobs0[ind4]),
size[ind4], prob[ind4], truncate = 0)
ans
}
rzabinom <- function(n, size, prob, pobs0 = 0) {
use.n <- if ((length.n <- length(n)) > 1) length.n else
if (!is.Numeric(n, integer.valued = TRUE,
length.arg = 1, positive = TRUE))
stop("bad input for argument 'n'") else n
ans <- rgaitdbinom(use.n, size, prob, truncate = 0)
if (length(pobs0) != use.n)
pobs0 <- rep_len(pobs0, use.n)
if (!is.Numeric(pobs0) || any(pobs0 < 0) || any(pobs0 > 1))
stop("argument 'pobs0' must be between 0 and 1 inclusive")
ifelse(runif(use.n) < pobs0, 0, ans)
}
zabinomial <-
function(lpobs0 = "logitlink",
lprob = "logitlink",
type.fitted = c("mean", "prob", "pobs0"),
ipobs0 = NULL, iprob = NULL,
imethod = 1,
zero = NULL # Was zero = 2 prior to 20130917
) {
lpobs0 <- as.list(substitute(lpobs0))
epobs0 <- link2list(lpobs0)
lpobs0 <- attr(epobs0, "function.name")
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0"))[1]
if (length(ipobs0))
if (!is.Numeric(ipobs0, positive = TRUE) ||
ipobs0 >= 1)
stop("argument 'ipobs0' is out of range")
if (length(iprob))
if (!is.Numeric(iprob, positive = TRUE) ||
iprob >= 1)
stop("argument 'iprob' is out of range")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
new("vglmff",
blurb = c("Zero-altered binomial distribution ",
"(Bernoulli and positive-binomial ",
"conditional model)\n\n",
"P[Y = 0] = pobs0,\n",
"P[Y = y] = (1 - pobs0) * dposbinom(x = y, size, prob), ",
"y = 1, 2, ..., size,\n\n",
"Link: ",
namesof("pobs0", lpobs0, earg = epobs0), ", ",
namesof("prob" , lprob, earg = eprob), "\n",
"Mean: (1 - pobs0) * prob / (1 - (1 - prob)^size)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = NA,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("pobs0", "prob"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted ))),
initialize = eval(substitute(expression({
if (!all(w == 1))
extra$orig.w <- w
if (NCOL(y) == 1) {
if (is.factor(y))
y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + w * y) / (1.0 + w)
no.successes <- y
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(no.successes - round(no.successes)) > 1.0e-8))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(y - round(y)) > 1.0e-8))
stop("Count data must be integer-valued")
y <- round(y)
nvec <- y[, 1] + y[, 2]
y <- ifelse(nvec > 0, y[, 1] / nvec, 0)
w <- w * nvec
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + nvec * y) / (1 + nvec)
} else {
stop("for the binomialff family, response 'y' must be a ",
"vector of 0 and 1's\n",
"or a factor ",
"(first level = fail, other levels = success),\n",
"or a 2-column matrix where col 1 is the no. of ",
"successes and col 2 is the no. of failures")
}
if (!all(w == 1))
extra$new.w <- w
y <- as.matrix(y)
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$skip.these <- skip.these <- matrix(as.logical(y0), n, NOS)
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
predictors.names <-
c(namesof("pobs0", .lpobs0 , .epobs0 , tag = FALSE),
namesof("prob" , .lprob , .eprob , tag = FALSE))
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
phi.init <- if (length( .ipobs0 )) .ipobs0 else {
prob0.est <- sum(Size[y == 0]) / sum(Size)
if ( .imethod == 1) {
(prob0.est - (1 - mustart)^Size) / (1 -
(1 - mustart)^Size)
} else
if ( .imethod == 2) {
prob0.est
} else {
prob0.est * 0.5
}
}
phi.init[phi.init <= -0.10] <- 0.50 # Lots of
phi.init[phi.init <= 0.01] <- 0.05 # Last resort
phi.init[phi.init >= 0.99] <- 0.95 # Last resort
if (!length(etastart)) {
etastart <-
cbind(theta2eta(phi.init, .lpobs0, earg = .epobs0 ),
theta2eta( mustart, .lprob, earg = .eprob ))
mustart <- NULL
}
}), list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0,
.iprob = iprob, .ipobs0 = ipobs0,
.imethod = imethod,
.type.fitted = type.fitted ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0"))[1]
phi0 <- eta2theta(eta[, 1], .lpobs0 , earg = .epobs0 )
prob <- eta2theta(eta[, 2], .lprob , earg = .eprob )
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
ans <- switch(type.fitted,
"mean" = (1 - phi0) * prob / (1 - (1 - prob)^Size),
"prob" = prob,
"pobs0" = phi0) # P(Y=0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0 ))),
last = eval(substitute(expression({
misc$link <- c(prob = .lprob, pobs0 = .lpobs0 )
misc$earg <- list(prob = .eprob, pobs0 = .epobs0 )
misc$imethod <- .imethod
misc$zero <- .zero
misc$expected <- TRUE
}), list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0,
.zero = zero,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
pobs0 <- eta2theta(eta[, 1], .lpobs0 , earg = .epobs0 )
prob <- eta2theta(eta[, 2], .lprob , earg = .eprob )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(orig.w) * dzabinom(x = round(y * Size),
size = Size,
prob = prob, pobs0 = pobs0,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0 ))),
vfamily = c("zabinomial"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
phi0 <- eta2theta(eta[, 1], .lpobs0 , earg = .epobs0 )
prob <- eta2theta(eta[, 2], .lprob , earg = .eprob )
okay1 <- all(is.finite(phi0)) && all(0 < phi0 & phi0 < 1) &&
all(is.finite(prob)) && all(0 < prob & prob < 1)
okay1
}, list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0 ))),
deriv = eval(substitute(expression({
M1 <- 2
NOS <- if (length(extra$NOS)) extra$NOS else 1
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
phi0 <- eta2theta(eta[, 1], .lpobs0 , earg = .epobs0 )
prob <- eta2theta(eta[, 2], .lprob , earg = .eprob )
dphi0.deta <- dtheta.deta(phi0, .lpobs0, earg = .epobs0 )
dprob.deta <- dtheta.deta(prob, .lprob , earg = .eprob )
df0.dprob <- -Size * (1 - prob)^(Size - 1)
df02.dprob2 <- Size * (Size - 1) * (1 - prob)^(Size - 2)
prob0 <- (1 - prob)^(Size)
oneminusf0 <- 1 - prob0
dl.dphi0 <- -1 / (1 - phi0)
dl.dprob <- c(w) * (y / prob - (1 - y) / (1 - prob)) +
c(orig.w) * df0.dprob / oneminusf0
dl.dphi0[y == 0] <- 1 / phi0[y == 0] # Do it in one line
skip <- extra$skip.these
for (spp. in 1:NOS) {
dl.dprob[skip[, spp.], spp.] <- 0
}
ans <- cbind(c(orig.w) * dl.dphi0 * dphi0.deta,
dl.dprob * dprob.deta)
ans
}), list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0 ))),
weight = eval(substitute(expression({
wz <- matrix(0.0, n, M1)
usualmeanY <- prob
meanY <- (1 - phi0) * usualmeanY / oneminusf0
term1 <- c(Size) * (meanY / prob^2 -
meanY / (1 - prob)^2) +
c(Size) * (1 - phi0) / (1 - prob)^2
term2 <- -(1 - phi0) * df02.dprob2 / oneminusf0
term3 <- -(1 - phi0) * (df0.dprob / oneminusf0)^2
ned2l.dprob2 <- term1 + term2 + term3
wz[, iam(2, 2, M)] <- ned2l.dprob2 * dprob.deta^2
mu.phi0 <- phi0
tmp100 <- mu.phi0 * (1.0 - mu.phi0)
tmp200 <- if ( .lpobs0 == "logitlink" &&
is.empty.list( .epobs0 )) {
tmp100
} else {
(dphi0.deta^2) / tmp100
}
wz[, iam(1, 1, M)] <- tmp200
c(orig.w) * wz
}), list( .lprob = lprob, .lpobs0 = lpobs0,
.eprob = eprob, .epobs0 = epobs0 ))))
} # zabinomial
zabinomialff <-
function(lprob = "logitlink",
lonempobs0 = "logitlink",
type.fitted = c("mean", "prob", "pobs0", "onempobs0"),
iprob = NULL, ionempobs0 = NULL,
imethod = 1,
zero = "onempobs0") {
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
lonempobs0 <- as.list(substitute(lonempobs0))
eonempobs0 <- link2list(lonempobs0)
lonempobs0 <- attr(eonempobs0, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "onempobs0"))[1]
if (length(iprob))
if (!is.Numeric(iprob, positive = TRUE) ||
iprob >= 1)
stop("argument 'iprob' is out of range")
if (length(ionempobs0))
if (!is.Numeric(ionempobs0, positive = TRUE) ||
ionempobs0 >= 1)
stop("argument 'ionempobs0' is out of range")
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
new("vglmff",
blurb = c("Zero-altered binomial distribution ",
"(Bernoulli and positive-binomial ",
"conditional model)\n\n",
"P[Y = 0] = 1 - onempobs0,\n",
"P[Y = y] = onempobs0 * dposbinom(x = y, size, prob), ",
"y = 1, 2, ..., size,\n\n",
"Link: ",
namesof("prob" , lprob , eprob ), ", ",
namesof("onempobs0", lonempobs0, eonempobs0), "\n",
"Mean: onempobs0 * prob / (1 - (1 - prob)^size)"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = NA,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("prob", "onempobs0"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted ))),
initialize = eval(substitute(expression({
if (!all(w == 1))
extra$orig.w <- w
if (NCOL(y) == 1) {
if (is.factor(y))
y <- y != levels(y)[1]
nn <- rep_len(1, n)
if (!all(y >= 0 & y <= 1))
stop("response values must be in [0, 1]")
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + w * y) / (1.0 + w)
no.successes <- y
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(no.successes - round(no.successes)) > 1.0e-8))
stop("Number of successes must be integer-valued")
} else if (NCOL(y) == 2) {
if (min(y) < 0)
stop("Negative data not allowed!")
if (any(abs(y - round(y)) > 1.0e-8))
stop("Count data must be integer-valued")
y <- round(y)
nvec <- y[, 1] + y[, 2]
y <- ifelse(nvec > 0, y[, 1] / nvec, 0)
w <- w * nvec
if (!length(mustart) && !length(etastart))
mustart <- (0.5 + nvec * y) / (1 + nvec)
} else {
stop("for the binomialff family, response 'y' must be a ",
"vector of 0 and 1's\n",
"or a factor ",
"(first level = fail, other levels = success),\n",
"or a 2-column matrix where col 1 is the no. of ",
"successes and col 2 is the no. of failures")
}
if (!all(w == 1))
extra$new.w <- w
y <- as.matrix(y)
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$skip.these <-
skip.these <- matrix(as.logical(y0), n, NOS)
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
predictors.names <-
c(namesof("prob" , .lprob , .eprob , tag = FALSE),
namesof("onempobs0", .lonempobs0 , .eonempobs0 , tag = FALSE))
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
phi.init <- if (length( .ionempobs0 ))
1 - .ionempobs0 else {
prob0.est <- sum(Size[y == 0]) / sum(Size)
if ( .imethod == 1) {
(prob0.est - (1 - mustart)^Size) / (1 -
(1 - mustart)^Size)
} else
if ( .imethod == 2) {
prob0.est
} else {
prob0.est * 0.5
}
}
phi.init[phi.init <= -0.10] <- 0.50 # Lots of
phi.init[phi.init <= 0.01] <- 0.05 # Last resort
phi.init[phi.init >= 0.99] <- 0.95 # Last resort
if (!length(etastart)) {
etastart <-
cbind(theta2eta( mustart, .lprob , .eprob ),
theta2eta(1 - phi.init, .lonempobs0 , .eonempobs0 ))
mustart <- NULL
}
}), list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0,
.iprob = iprob, .ionempobs0 = ionempobs0,
.imethod = imethod,
.type.fitted = type.fitted ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
NOS <- ncol(eta) / c(M1 = 2)
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "onempobs0"))[1]
prob <- eta2theta(eta[, 1], .lprob , .eprob )
onempobs0 <- eta2theta(eta[, 2], .lonempobs0 , .eonempobs0 )
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
ans <- switch(type.fitted,
"mean" = onempobs0 * prob / (1 - (1 - prob)^Size),
"prob" = prob,
"pobs0" = 1 - onempobs0, # P(Y=0)
"onempobs0" = onempobs0) # P(Y>0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0 ))),
last = eval(substitute(expression({
misc$link <- c(prob = .lprob , onempobs0 = .lonempobs0 )
misc$earg <- list(prob = .eprob , onempobs0 = .eonempobs0 )
misc$imethod <- .imethod
misc$zero <- .zero
misc$expected <- TRUE
}), list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0,
.zero = zero,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
prob <- eta2theta(eta[, 1], .lprob , .eprob )
onempobs0 <- eta2theta(eta[, 2], .lonempobs0 , .eonempobs0 )
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
orig.w * dzabinom(x = round(y * Size), size = Size,
prob = prob, pobs0 = 1 - onempobs0,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0 ))),
vfamily = c("zabinomialff"),
validparams = eval(substitute(function(eta, y, extra = NULL) {
prob <- eta2theta(eta[, 1], .lprob , .eprob )
onempobs0 <- eta2theta(eta[, 2], .lonempobs0 , .eonempobs0 )
okay1 <- all(is.finite(onempobs0)) &&
all(0 < onempobs0 & onempobs0 < 1) &&
all(is.finite(prob )) &&
all(0 < prob & prob < 1)
okay1
}, list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0 ))),
deriv = eval(substitute(expression({
M1 <- 2
NOS <- if (length(extra$NOS)) extra$NOS else 1
orig.w <- if (length(extra$orig.w)) extra$orig.w else 1
new.w <- if (length(extra$new.w)) extra$new.w else 1
Size <- round(new.w / orig.w)
prob <- eta2theta(eta[, 1], .lprob , .eprob )
onempobs0 <- eta2theta(eta[, 2], .lonempobs0 , .eonempobs0 )
phi0 <- 1 - onempobs0
dprob.deta <- dtheta.deta(prob , .lprob ,
earg = .eprob )
donempobs0.deta <- dtheta.deta(onempobs0, .lonempobs0 ,
earg = .eonempobs0 )
df0.dprob <- -Size * (1 - prob)^(Size - 1)
df02.dprob2 <- Size * (Size - 1) * (1 - prob)^(Size - 2)
prob0 <- (1 - prob)^(Size)
oneminusf0 <- 1 - prob0
dl.dprob <- c(w) * (y / prob - (1 - y) / (1 - prob)) +
c(orig.w) * df0.dprob / oneminusf0
dl.donempobs0 <- +1 / (onempobs0)
dl.donempobs0[y == 0] <-
-1 / (1 - onempobs0[y == 0]) # Do it in 1 line
skip <- extra$skip.these
for (spp. in 1:NOS) {
dl.dprob[skip[, spp.], spp.] <- 0
}
ans <- cbind( dl.dprob * dprob.deta,
c(orig.w) * dl.donempobs0 * donempobs0.deta)
ans
}), list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0 ))),
weight = eval(substitute(expression({
wz <- matrix(0.0, n, M1)
usualmeanY <- prob
meanY <- (1 - phi0) * usualmeanY / oneminusf0
term1 <- c(Size) * (meanY / prob^2 -
meanY / (1 - prob)^2) +
c(Size) * (1 - phi0) / (1 - prob)^2
term2 <- -(1 - phi0) * df02.dprob2 / oneminusf0
term3 <- -(1 - phi0) * (df0.dprob / oneminusf0)^2
ned2l.dprob2 <- term1 + term2 + term3
wz[, iam(1, 1, M)] <- ned2l.dprob2 * dprob.deta^2
mu.phi0 <- phi0
tmp100 <- mu.phi0 * (1.0 - mu.phi0)
tmp200 <- if (FALSE &&
.lonempobs0 == "logitlink" &&
is.empty.list( .eonempobs0 )) {
tmp100
} else {
(donempobs0.deta^2) / tmp100
}
wz[, iam(2, 2, M)] <- tmp200
c(orig.w) * wz
}),
list( .lprob = lprob, .lonempobs0 = lonempobs0,
.eprob = eprob, .eonempobs0 = eonempobs0 ))))
}
zageometric <-
function(lpobs0 = "logitlink", lprob = "logitlink",
type.fitted = c("mean", "prob", "pobs0", "onempobs0"),
imethod = 1,
ipobs0 = NULL, iprob = NULL,
zero = NULL) {
lpobs0 <- as.list(substitute(lpobs0))
epobs0 <- link2list(lpobs0)
lpobs0 <- attr(epobs0, "function.name")
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "onempobs0"))[1]
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
if (length(iprob))
if (!is.Numeric(iprob, positive = TRUE) ||
max(iprob) >= 1)
stop("argument 'iprob' out of range")
if (length(ipobs0))
if (!is.Numeric(ipobs0, positive = TRUE) ||
max(ipobs0) >= 1)
stop("argument 'ipobs0' out of range")
new("vglmff",
blurb = c("Zero-altered geometric ",
"(Bernoulli and positive-geometric conditional ",
"model)\n\n",
"Links: ",
namesof("pobs0", lpobs0, epobs0, tag = FALSE), ", ",
namesof("prob" , lprob , eprob , tag = FALSE), "\n",
"Mean: (1 - pobs0) / prob"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = FALSE,
parameters.names = c("pobs0", "prob"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted
))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$skip.these <- skip.these <-
matrix(as.logical(y0), n, NOS)
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("pobs0", ncoly, skip1 = TRUE)
mynames2 <- param.names("prob", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lpobs0 , .epobs0 , tag = FALSE),
namesof(mynames2, .lprob , .eprob , tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
if (!length(etastart)) {
foo <- function(x) mean(as.numeric(x == 0))
phi0.init <- matrix(apply(y, 2, foo),
n, ncoly, byrow = TRUE)
if (length( .ipobs0 ))
phi0.init <- matrix( .ipobs0 , n, ncoly, byrow = TRUE)
prob.init <-
if ( .imethod == 2)
1 / (1 + y + 1/16) else
if ( .imethod == 1)
(1 - phi0.init) / (1 +
matrix(colSums(y * w) / colSums(w) + 1/16,
n, ncoly, byrow = TRUE)) else
(1 - phi0.init) / (1 +
matrix(apply(y, 2, median),
n, ncoly, byrow = TRUE) + 1/16)
if (length( .iprob ))
prob.init <- matrix( .iprob , n, ncoly, byrow = TRUE)
etastart <- cbind(theta2eta(phi0.init, .lpobs0 , .epobs0 ),
theta2eta(prob.init, .lprob , .eprob ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
}
}), list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob,
.ipobs0 = ipobs0, .iprob = iprob,
.imethod = imethod,
.type.fitted = type.fitted ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "onempobs0"))[1]
M1 <- 2
NOS <- ncol(eta) / M1
phi0 <- cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lpobs0 , earg = .epobs0 ))
prob <- cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lprob , earg = .eprob ))
ans <- switch(type.fitted,
"mean" = (1 - phi0) / prob,
"prob" = prob,
"pobs0" = phi0, # P(Y=0)
"onempobs0" = 1 - phi0) # P(Y>0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob ))),
last = eval(substitute(expression({
temp.names <- c(rep_len( .lpobs0 , NOS),
rep_len( .lprob , NOS))
temp.names <- temp.names[interleave.VGAM(M1*NOS, M1 = M1)]
misc$link <- temp.names
misc$earg <- vector("list", M1 * NOS)
names(misc$link) <-
names(misc$earg) <-
c(mynames1, mynames2)[interleave.VGAM(M1*NOS, M1 = M1)]
for (ii in 1:NOS) {
misc$earg[[M1*ii-1]] <- .epobs0
misc$earg[[M1*ii ]] <- .eprob
}
misc$expected <- TRUE
misc$imethod <- .imethod
misc$ipobs0 <- .ipobs0
misc$iprob <- .iprob
misc$multipleResponses <- TRUE
}), list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob,
.ipobs0 = ipobs0, .iprob = iprob,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
NOS <- extra$NOS
M1 <- 2
phi0 <- cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lpobs0 , earg = .epobs0 ))
prob <- cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lprob , earg = .eprob ))
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzageom(y, pobs0 = phi0, prob = prob,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob ))),
vfamily = c("zageometric"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
phi0 <- cbind(eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
.lpobs0 , earg = .epobs0 ))
prob <- cbind(eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
.lprob , earg = .eprob ))
rzageom(nsim * length(prob), prob = prob, pobs0 = phi0)
}, list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
phi0 <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
.lpobs0 , earg = .epobs0 )
prob <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
.lprob , earg = .eprob )
okay1 <- all(is.finite(phi0)) && all(0 < phi0 & phi0 < 1) &&
all(is.finite(prob)) && all(0 < prob & prob < 1)
okay1
}, list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob ))),
deriv = eval(substitute(expression({
M1 <- 2
NOS <- ncol(eta) / M1 # extra$NOS
y0 <- extra$y0
skip <- extra$skip.these
phi0 <- cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lpobs0 , earg = .epobs0 ))
prob <- cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lprob , earg = .eprob ))
dl.dprob <- 1 / prob - (y - 1) / (1 - prob)
dl.dphi0 <- -1 / (1 - phi0)
for (spp. in 1:NOS) {
dl.dphi0[skip[, spp.], spp.] <- 1 / phi0[skip[, spp.], spp.]
dl.dprob[skip[, spp.], spp.] <- 0
}
dphi0.deta <- dtheta.deta(phi0, .lpobs0 , earg = .epobs0 )
dprob.deta <- dtheta.deta(prob, .lprob , earg = .eprob )
ans <- c(w) * cbind(dl.dphi0 * dphi0.deta,
dl.dprob * dprob.deta)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
ans
}), list( .lpobs0 = lpobs0, .lprob = lprob,
.epobs0 = epobs0, .eprob = eprob ))),
weight = eval(substitute(expression({
wz <- matrix(0.0, n, M1*NOS)
ned2l.dprob2 <- (1 - phi0) / (prob^2 * (1 - prob))
wz[, NOS+(1:NOS)] <- c(w) * ned2l.dprob2 * dprob.deta^2
mu.phi0 <- phi0
tmp100 <- mu.phi0 * (1.0 - mu.phi0)
tmp200 <- if ( .lpobs0 == "logitlink" &&
is.empty.list( .epobs0 )) {
cbind(c(w) * tmp100)
} else {
cbind(c(w) * (dphi0.deta^2) / tmp100)
}
wz[, 1:NOS] <- tmp200
wz <- wz[, interleave.VGAM(ncol(wz), M1 = M1)]
wz
}), list( .lpobs0 = lpobs0,
.epobs0 = epobs0 ))))
} # End of zageometric
zageometricff <-
function(lprob = "logitlink", lonempobs0 = "logitlink",
type.fitted = c("mean", "prob", "pobs0", "onempobs0"),
imethod = 1,
iprob = NULL, ionempobs0 = NULL,
zero = "onempobs0") {
lprob <- as.list(substitute(lprob))
eprob <- link2list(lprob)
lprob <- attr(eprob, "function.name")
lonempobs0 <- as.list(substitute(lonempobs0))
eonempobs0 <- link2list(lonempobs0)
lonempobs0 <- attr(eonempobs0, "function.name")
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "onempobs0"))[1]
if (!is.Numeric(imethod, length.arg = 1,
integer.valued = TRUE, positive = TRUE) ||
imethod > 3)
stop("argument 'imethod' must be 1 or 2 or 3")
if (length(iprob))
if (!is.Numeric(iprob, positive = TRUE) ||
max(iprob) >= 1)
stop("argument 'iprob' out of range")
if (length(ionempobs0))
if (!is.Numeric(ionempobs0, positive = TRUE) ||
max(ionempobs0) >= 1)
stop("argument 'ionempobs0' out of range")
new("vglmff",
blurb = c("Zero-altered geometric ",
"(Bernoulli and positive-geometric ",
"conditional model)\n\n",
"Links: ",
namesof("prob" , lprob , eprob ,
tag = FALSE), ", ",
namesof("onempobs0", lonempobs0, eonempobs0,
tag = FALSE), "\n",
"Mean: onempobs0 / prob"),
constraints = eval(substitute(expression({
constraints <- cm.zero.VGAM(constraints, x = x, .zero ,
M = M, M1 = 2,
predictors.names = predictors.names)
}), list( .zero = zero ))),
infos = eval(substitute(function(...) {
list(M1 = 2,
Q1 = 1,
expected = TRUE,
multipleResponses = TRUE,
parameters.names = c("prob", "onempobs0"),
type.fitted = .type.fitted ,
zero = .zero )
}, list( .zero = zero,
.type.fitted = type.fitted
))),
initialize = eval(substitute(expression({
M1 <- 2
temp5 <-
w.y.check(w = w, y = y,
Is.nonnegative.y = TRUE,
Is.integer.y = TRUE,
ncol.w.max = Inf,
ncol.y.max = Inf,
out.wy = TRUE,
colsyperw = 1,
maximize = TRUE)
w <- temp5$w
y <- temp5$y
extra$y0 <- y0 <- ifelse(y == 0, 1, 0)
extra$NOS <- NOS <- ncoly <- ncol(y) # Number of species
extra$skip.these <-
skip.these <- matrix(as.logical(y0), n, NOS)
extra$type.fitted <- .type.fitted
extra$colnames.y <- colnames(y)
mynames1 <- param.names("prob", ncoly, skip1 = TRUE)
mynames2 <- param.names("onempobs0", ncoly, skip1 = TRUE)
predictors.names <-
c(namesof(mynames1, .lprob , earg = .eprob ,
tag = FALSE),
namesof(mynames2, .lonempobs0 , earg = .eonempobs0 ,
tag = FALSE))[
interleave.VGAM(M1*NOS, M1 = M1)]
if (!length(etastart)) {
foo <- function(x) mean(as.numeric(x == 0))
phi0.init <- matrix(apply(y, 2, foo),
n, ncoly, byrow = TRUE)
if (length( .ionempobs0 ))
phi0.init <- matrix( 1 - .ionempobs0 ,
n, ncoly, byrow = TRUE)
prob.init <-
if ( .imethod == 2)
1 / (1 + y + 1/16) else
if ( .imethod == 1)
(1 - phi0.init) / (1 +
matrix(colSums(y * w) / colSums(w) + 1/16,
n, ncoly, byrow = TRUE)) else
(1 - phi0.init) / (1 +
matrix(apply(y, 2, median), n, ncoly, byrow = TRUE) + 1/16)
if (length( .iprob ))
prob.init <- matrix( .iprob , n, ncoly, byrow = TRUE)
etastart <-
cbind(theta2eta( prob.init, .lprob , .eprob ),
theta2eta(1 - phi0.init, .lonempobs0 , .eonempobs0 ))
etastart <- etastart[, interleave.VGAM(ncol(etastart),
M1 = M1)]
}
}), list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob,
.ionempobs0 = ionempobs0, .iprob = iprob,
.imethod = imethod,
.type.fitted = type.fitted ))),
linkinv = eval(substitute(function(eta, extra = NULL) {
type.fitted <- if (length(extra$type.fitted))
extra$type.fitted else {
warning("cannot find 'type.fitted'. ",
"Returning the 'mean'.")
"mean"
}
type.fitted <- match.arg(type.fitted,
c("mean", "prob", "pobs0", "onempobs0"))[1]
M1 <- 2
NOS <- ncol(eta) / M1
prob <-
cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lprob , earg = .eprob ))
onempobs0 <-
cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lonempobs0 , earg = .eonempobs0 ))
ans <- switch(type.fitted,
"mean" = onempobs0 / prob,
"prob" = prob,
"pobs0" = 1 - onempobs0, # P(Y=0)
"onempobs0" = onempobs0) # P(Y>0)
label.cols.y(ans, colnames.y = extra$colnames.y, NOS = NOS)
}, list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob ))),
last = eval(substitute(expression({
temp.names <- c(rep_len( .lprob , NOS),
rep_len( .lonempobs0 , NOS))
temp.names <- temp.names[interleave.VGAM(M1*NOS, M1 = M1)]
misc$link <- temp.names
misc$earg <- vector("list", M1 * NOS)
names(misc$link) <-
names(misc$earg) <-
c(mynames1, mynames2)[interleave.VGAM(M1*NOS, M1 = M1)]
for (ii in 1:NOS) {
misc$earg[[M1*ii-1]] <- .eprob
misc$earg[[M1*ii ]] <- .eonempobs0
}
misc$expected <- TRUE
misc$imethod <- .imethod
misc$ionempobs0 <- .ionempobs0
misc$iprob <- .iprob
misc$multipleResponses <- TRUE
}), list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob,
.ionempobs0 = ionempobs0, .iprob = iprob,
.imethod = imethod ))),
loglikelihood = eval(substitute(
function(mu, y, w, residuals = FALSE, eta,
extra = NULL,
summation = TRUE) {
NOS <- extra$NOS
M1 <- 2
prob <-
cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lprob , earg = .eprob ))
onempobs0 <-
cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lonempobs0 , earg = .eonempobs0 ))
if (residuals) {
stop("loglikelihood residuals not implemented yet")
} else {
ll.elts <-
c(w) * dzageom(y, pobs0 = 1 - onempobs0, prob = prob,
log = TRUE)
if (summation) {
sum(ll.elts)
} else {
ll.elts
}
}
}, list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob ))),
vfamily = c("zageometricff"),
simslot = eval(substitute(
function(object, nsim) {
pwts <- if (length(pwts <- object@prior.weights) > 0)
pwts else weights(object, type = "prior")
if (any(pwts != 1))
warning("ignoring prior weights")
eta <- predict(object)
onempobs0 <-
cbind(eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
.lonempobs0 , earg = .eonempobs0 ))
prob <-
cbind(eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
.lprob , earg = .eprob ))
rzageom(nsim * length(prob),
pobs0 = 1 - onempobs0, prob = prob)
}, list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob ))),
validparams = eval(substitute(function(eta, y, extra = NULL) {
prob <- eta2theta(eta[, c(TRUE, FALSE), drop = FALSE],
.lprob , earg = .eprob )
onempobs0 <- eta2theta(eta[, c(FALSE, TRUE), drop = FALSE],
.lonempobs0 , earg = .eonempobs0 )
okay1 <- all(is.finite(onempobs0)) &&
all(0 < onempobs0 & onempobs0 < 1) &&
all(is.finite(prob )) &&
all(0 < prob & prob < 1)
okay1
}, list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob ))),
deriv = eval(substitute(expression({
M1 <- 2
NOS <- ncol(eta) / M1 # extra$NOS
y0 <- extra$y0
skip <- extra$skip.these
prob <-
cbind(eta2theta(eta[, M1*(1:NOS)-1, drop = FALSE],
.lprob , earg = .eprob ))
onempobs0 <-
cbind(eta2theta(eta[, M1*(1:NOS)-0, drop = FALSE],
.lonempobs0 , earg = .eonempobs0 ))
pobs0 <- 1 - onempobs0
dl.dprob <- 1 / prob - (y - 1) / (1 - prob)
dl.donempobs0 <- +1 / (onempobs0)
for (spp. in 1:NOS) {
dl.donempobs0[skip[, spp.], spp.] <-
-1 / pobs0[skip[, spp.], spp.]
dl.dprob[skip[, spp.], spp.] <- 0
}
dprob.deta <- dtheta.deta(prob,
.lprob , .eprob )
donempobs0.deta <- dtheta.deta(onempobs0,
.lonempobs0 , .eonempobs0 )
ans <- c(w) * cbind(dl.dprob * dprob.deta,
dl.donempobs0 * donempobs0.deta)
ans <- ans[, interleave.VGAM(ncol(ans), M1 = M1)]
ans
}), list( .lonempobs0 = lonempobs0, .lprob = lprob,
.eonempobs0 = eonempobs0, .eprob = eprob ))),
weight = eval(substitute(expression({
wz <- matrix(0.0, n, M1*NOS)
ned2l.dprob2 <- (1 - pobs0) / (prob^2 * (1 - prob))
wz[, (1:NOS)] <- c(w) * ned2l.dprob2 * dprob.deta^2
mu.phi0 <- pobs0 # phi0
tmp100 <- mu.phi0 * (1.0 - mu.phi0)
tmp200 <- if ( FALSE &&
.lonempobs0 == "logitlink" &&
is.empty.list( .eonempobs0 )) {
cbind(c(w) * tmp100)
} else {
cbind(c(w) * (donempobs0.deta^2) / tmp100)
}
wz[, NOS+(1:NOS)] <- tmp200
wz <- wz[, interleave.VGAM(ncol(wz), M1 = M1)]
wz
}), list( .lonempobs0 = lonempobs0,
.eonempobs0 = eonempobs0 ))))
} # End of zageometricff
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.