diff --git a/R/method.R b/R/method.R index c4eb1b5..d110e5e 100644 --- a/R/method.R +++ b/R/method.R @@ -199,33 +199,35 @@ method.CC_nloglik <- function() { plogis(trimLogit(predY, trim = control$trimLogit) %*% matrix(coef)) } computeCoef = function(Z, Y, libraryNames, obsWeights, control, verbose, ...) { - cvRisk <- apply(Z, 2, function(x) { -sum(2 * obsWeights * ifelse(Y, log(x), log(1-x))) } ) + logitZ = trimLogit(Z, control$trimLogit) + cvRisk <- apply(logitZ, 2, function(x) -sum(2 * obsWeights * + ifelse(Y, plogis(x, log.p=TRUE), + plogis(x, log.p=TRUE, lower.tail=FALSE)))) names(cvRisk) <- libraryNames obj_and_grad <- function(y,x, w=NULL) { y <- y x <- x function(beta) { - xB <- x %*% cbind(beta) + xB <- x %*% cbind(beta) loglik <- y * plogis(xB, log.p=TRUE) + (1-y) * plogis(xB, log.p=TRUE, lower.tail=FALSE) if (!is.null(w)) loglik <- loglik * w obj <- -2 * sum(loglik) - p <- plogis(xB) grad <- if (is.null(w)) 2 * crossprod(x, cbind(p - y)) else 2 * crossprod(x, w*cbind(p - y)) list(objective=obj, gradient=grad) } - } + } r <- nloptr(x0=rep(1/ncol(Z), ncol(Z)), - eval_f=obj_and_grad(Y, Z), + eval_f=obj_and_grad(Y, logitZ), lb=rep(0, ncol(Z)), ub=rep(1, ncol(Z)), eval_g_eq = function(beta) (sum(beta)-1), eval_jac_g_eq = function(beta) rep(1, length(beta)), opts=list("algorithm"="NLOPT_LD_SLSQP","ftol_rel"=1.0e-8)) - if (r$status < 1 || r$status > 4) { - warning(r$messate) + if (r$status < 1 || r$status > 4) { + warning(r$message) } coef <- r$solution if (any(is.na(coef))) {