R/kl.compreg2.R

Defines functions kl.compreg2

Documented in kl.compreg2

kl.compreg2 <- function(y, x, xnew = NULL, tol = 1e-07, maxiters = 50) {
  X <- model.matrix(y~., data.frame(x))
  p <- dim(X)[2]
  Y <- y[, -1]
  dm <- dim(Y)
  n <- dm[1]
  d <- dm[2]
  m <- Rfast::colmeans(Y)
  b0 <- Rfast::Log(m / (1 - m) )
  b1 <- matrix( c(b0, numeric(p * d - d) ), nrow = p, ncol = d, byrow = TRUE)
  e <- Y - rep(m, rep(n, d) )
  id <- matrix(1:c(p * d), ncol = d)
  der <- numeric(d * p)
  der2 <- matrix(0, p * d, p * d)
  for (i in 1:d) {
    der[id[, i]] <- Rfast::colsums( e[, i] * X )
    for (j in i:d) {
      if (i != j) {
        der2[id[, i], id[, j]] <- der2[id[, j], id[, i]] <-  - crossprod(m[i] * m[j] * X, X)
      } else  der2[id[, i], id[, i]] <- crossprod(m[i] * (1 - m[i]) * X, X)
    }
  }
  b2 <- b1 + solve(der2, der)
  k <- 2
  res <- try(
  while ( sum( abs(b2 - b1) ) > tol  &  k < maxiters) {
    k <- k + 1
    b1 <- b2
    m1 <- exp(X %*% b1)
    m <- m1 / (Rfast::rowsums(m1) + 1)
    e <- Y - m
    for (i in 1:d) {
      der[id[, i]] <- Rfast::colsums( e[, i] * X )
      for (j in i:d) {
        if (i != j) {
          der2[id[, i], id[, j]] <- der2[id[, j], id[, i]] <-   - crossprod(m[, i] * m[, j] * X, X)
        } else  der2[id[, i], id[, i]] <- crossprod(m[, i] * (1 - m[, i]) * X, X)
      }
    }
    b2 <- b1 + solve(der2, der)
  },
  silent = TRUE)
  if ( class(res) == "try-error" )   b2 <- b1
  est <- NULL
  if ( !is.null(xnew) ) {
    xnew <- model.matrix(~., data.frame(xnew))
    mu <- cbind(1, exp(xnew %*% b2))
    est <- mu/Rfast::rowsums(mu)
  }
  colnames(b2) <- paste("Y", 1:d, sep = "")
  rownames(b2) <- colnames(X)
  m <- cbind(1, m1)
  m <- m / Rfast::rowsums(m)
  loglik <-  - sum(y * log(y/m), na.rm = TRUE )
  list(iters = k, loglik = loglik, be = b2, est = est)
}

Try the Compositional package in your browser

Any scripts or data that you put into this service are public.

Compositional documentation built on Jan. 14, 2018, 5:03 p.m.