logbin.em <- function(mt, mf, Y, offset, mono, start, control, accelerate = c("em","squarem","pem","qn"),
control.method, warn)
{
accelerate = match.arg(accelerate)
control2 <- control
control2$trace <- (control$trace > 1)
reparam <- logbin.reparam(mt, mf, "em", mono)
if (!is.null(start)) {
start.expand <- logbin.expand(start, reparam, "em")
}
if(control$trace > 0) cat("logbin parameterisation 1/1\n")
if (length(reparam$Vmat) == 0) {
X <- model.matrix(mt, mf)
Amat <- diag(ncol(X))
}
else {
des <- logbin.design(mt, mf, "em", reparam)
X <- des$X.reparam
Amat <- des$A
}
thismodel <- nplbin(Y, X, offset, if (!is.null(start)) start.expand$coefs.exp else NULL,
Amat = Amat, control = control2, accelerate = accelerate,
control.accelerate = list(control.method))
if(control$trace > 0 & control$trace <= 1)
cat("Deviance =", thismodel$deviance, "Iterations -", thismodel$iter, "\n")
if (length(reparam$Vmat) == 0) {
np.coefs <- coefs <- coefs.boundary <- thismodel$coefficients
nn.design <- design <- X
if (control$coeftrace) coefhist <- thismodel$coefhist
} else {
np.coefs <- thismodel$coefficients
nn.design <- X
coefs <- as.vector(logbin.reduce(np.coefs, des$A))
names(coefs) <- gsub("`", "", colnames(des$X.orig))
design <- des$X.orig
coefs.boundary <- np.coefs[logbin.expand(coefs, reparam, "em")$which.boundary]
if (control$coeftrace) {
coefhist <- coefhist.reduce(thismodel$coefhist, des$A, names(coefs))
}
}
nvars <- length(coefs)
vardiff <- length(np.coefs) - nvars
aic.c <- thismodel$aic - 2 * vardiff + 2 * nvars * (nvars + 1) / (NROW(Y) - nvars - 1)
boundary <- any(coefs.boundary > -control$bound.tol)
if (warn) {
if (!thismodel$converged) {
if (identical(accelerate, "em"))
warning(gettextf("nplbin: algorithm did not converge within %d iterations -- increase 'maxit'.",
control$maxit), call. = FALSE)
else
warning(gettextf("nplbin(%s): algorithm did not converge within %d iterations -- increase 'maxit' or try with 'accelerate = \"em\"'.",
accelerate, control$maxit), call. = FALSE)
}
if(boundary) {
if(coefs.boundary[1] > -control$bound.tol)
warning("nplbin: fitted probabilities numerically 1 occurred", call. = FALSE)
else
warning("nplbin: MLE on boundary of constrained parameter space", call. = FALSE)
}
}
fit <- list(coefficients = coefs, residuals = thismodel$residuals,
fitted.values = thismodel$fitted.values,
linear.predictors = thismodel$linear.predictors, deviance = thismodel$deviance,
loglik = thismodel$loglik, aic = thismodel$aic - 2*vardiff, aic.c = aic.c,
null.deviance = thismodel$null.deviance, iter = thismodel$iter,
prior.weights = thismodel$prior.weights,
df.residual = thismodel$df.residual + vardiff, df.null = thismodel$df.null,
y = thismodel$y, x = design, converged = thismodel$converged,
boundary = boundary, np.coefficients = np.coefs,
nn.x = nn.design)
if (control$coeftrace) fit$coefhist <- coefhist
fit
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.