nplbin <- function(y, x, offset, start, Amat = diag(ncol(x)), control = logbin.control(),
accelerate = c("em","squarem","pem","qn"), control.accelerate = list(list())) {
control <- do.call("logbin.control", control)
accelerate <- match.arg(accelerate)
x <- as.matrix(x)
xnames <- dimnames(x)[[2L]]
ynames <- if (is.matrix(y))
rownames(y)
else names(y)
if (any(x < 0)) stop("x must be non-negative")
if (any(apply(x, 2, function(col) all(col==0)))) stop("x contains column with all 0")
nvars <- ncol(x)
nobs <- NROW(y)
n <- weights <- rep(1, nobs)
if (is.null(offset)) offset <- rep.int(0, nobs)
if (any(offset > 0))
stop("offset must be non-positive")
fam <- binomial(link = log)
eval(fam$initialize)
if (max(y) == 0) stop("cannot fit a log-binomial model if all y==0")
mu.eta <- fam$mu.eta
linkinv <- fam$linkinv
dev.resids <- fam$dev.resids
aic <- fam$aic
y1 <- round(n*y)
y2 <- round(n*(1-y))
x.mins <- apply(x, 2, function(t) min(t[t > 0]))
x.scale <- 1/x.mins
converged <- FALSE
coefold <- if (!is.null(start)) {
if (length(start) != nvars)
stop(gettextf("length of 'start' should equal %d and correspond to initial coefs for %s",
nvars, paste(deparse(xnames), collapse = ", ")), domain = NA)
else if (any(start >= -control$bound.tol))
stop("'start' is on or outside the boundary of the parameter space (consider 'bound.tol')", domain = NA)
else start
} else {
simple <- log(mean(y)) / colMeans(x) - 2 * control$bound.tol
logy <- log(y)
logy[is.infinite(logy)] <- min(logy[!is.infinite(logy)]) - 2 * control$bound.tol
trymat <- tryCatch(as.vector(solve(t(x) %*% x) %*% t(x) %*% (logy)) + 2*control$bound.tol,
error = function(e) NULL)
if (is.null(trymat)) simple
else if (any(trymat >= -control$bound.tol)) simple
else trymat
}
fixptfn <- function(p, y1, y2, n, x, x.s, o, nobs, nvars, fam, bound.tol) {
x1 <- sweep(x, 2, FUN = "*", x.s)
p.scale <- p / x.s
eta <- drop(x1 %*% p.scale) + o
estep <- y1 + y2 * ((matrix(fam$linkinv(p.scale), nobs, nvars, byrow = TRUE) - fam$linkinv(eta))/(1 - fam$linkinv(eta)))
pnew.scale <- log(colSums(estep * x1) / colSums(n * x1))
pnew <- pnew.scale * x.s
pnew[pnew >= 0] <- -bound.tol / 2
# Sometimes estimates go to negative infinity (so exp(p) = 0). Just reduce them by a bit.
pnew[is.infinite(pnew)] <- pmin(p[is.infinite(pnew)] - .Machine$double.neg.eps,
log(.Machine$double.neg.eps))
return(pnew)
}
objfn <- function(p, y1, y2, n, x, x.s, o, nobs, nvars, fam, bound.tol) {
eta <- drop(x %*% p) + o
mu <- n * fam$linkinv(eta)
negll <- -sum(dbinom(y1, size = n, prob = mu / n, log = TRUE))
return(negll)
}
validparams <- function(p) return(all(p <= 0))
projfn <- function(p) {
pnew <- p
pnew[p > 0] <- -.Machine$double.eps
}
conv.user <- conv.test(Amat, control$epsilon)
emargs <- list(par = coefold, fixptfn = fixptfn, objfn = objfn, method = accelerate,
pconstr = validparams, project = projfn, y1 = y1, y2 = y2, n = n, x = x, x.s = x.scale,
o = offset, nobs = nobs, nvars = nvars, fam = fam, bound.tol = control$bound.tol,
control.run = list(convtype = "parameter", tol = control$epsilon,
stoptype = "maxiter", maxiter = control$maxit,
convfn.user = conv.user, trace = control$trace,
keep.paramval = control$coeftrace),
control.method = control.accelerate)
res <- do.call(turboEM::turboem, emargs)
if (res$fail[1]) stop(res$errors[1])
iter <- res$itr
if (control$coeftrace) {
coefhist <- res$trace.paramval[[1]]
colnames(coefhist) <- xnames
}
coefnew <- res$pars[1,]
names(coefnew) <- xnames
eta <- drop(x %*% coefnew) + offset
mu <- n * linkinv(eta)
residuals <- (y - (mu / n)) / mu.eta(eta)
names(y) <- names(mu) <- names(eta) <- names(residuals) <- ynames
dev.new <- sum(dev.resids(y, mu / n, n))
aic.model <- aic(y, n, mu / n, weights, dev.new) + 2 * nvars
aic.c.model <- aic.model + 2 * nvars * (nvars + 1) / (nobs - nvars - 1)
wtdmu <- sum(n * y) / sum(n)
nulldev <- sum(dev.resids(y, wtdmu, n))
nulldf <- nobs - 1
resdf <- nobs - nvars
boundary <- any(coefnew > -control$bound.tol)
result <- list(coefficients = coefnew, residuals = residuals, fitted.values = mu / n,
rank = nvars, family = fam, linear.predictors = eta, deviance = dev.new,
aic = aic.model, aic.c = aic.c.model, null.deviance = nulldev, iter = iter,
prior.weights = n, df.residual = resdf, df.null = nulldf, y = y,
converged = res$convergence[1], boundary = boundary, loglik = -res$value.objfn[1],
nn.design = x)
if (control$coeftrace) result$coefhist <- coefhist
result
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.