R/js.compreg.R

Defines functions js.compreg

Documented in js.compreg

################################
#### Jensen-Shannon divergence based regression for compositional data
#### Tsagris Michail 5/2015
#### [email protected]
#### References: Michail Tsagris (2015)
#### A novel, divergence based, regression for compositional data
#### Proceedings of the 28th Panhellenic Statistics Conference
################################
js.compreg <- function(y, x, B = 1, ncores = 1, xnew = NULL) {
  ## y is dependent variable, the compositional data
  ## x is the independent variable(s)
  ## B is the number of bootstrap samples used to obtain
  ## standard errors for the betas
  ## if B==1 no bootstrap is performed and no standard errors are reported
  ## if ncores=1, then 1 processor is used, otherwise
  ## more are used (parallel computing)
  n <- dim(y)[1]  ## sample size
  x <- model.matrix(y ~ ., data.frame(x) )
  d <- dim(y)[2] - 1  ## dimensionality of the simplex

  jsreg <- function(para, y, x, d){
    be <- matrix(para, byrow = TRUE, ncol = d)
    mu1 <- cbind( 1, exp(x %*% be) )
    mu <- mu1 / rowSums(mu1)
    M <- ( mu + y ) / 2
    sum( - y * log1p(mu / y) + mu * log(mu / M), na.rm = TRUE )
  }

  ## the next lines minimize the kl.compreg function and obtain the estimated betas
  runtime <- proc.time()
  ini <- as.vector( t( kl.compreg(y, x[, -1, drop = FALSE])$be ) )
  options (warn = -1)
  qa <- nlm(jsreg, ini, y = y, x = x, d = d)
  qa <- nlm(jsreg, qa$estimate, y = y, x = x, d = d)
  qa <- nlm(jsreg, qa$estimate, y = y, x = x, d = d)
  be <- matrix(qa$estimate, byrow = TRUE, ncol = d)
  seb <- NULL
  runtime <- proc.time() - runtime

  if (B > 1) {
  betaboot <- matrix( nrow = B, ncol = length(ini) )
  nc <- ncores
    if (nc == 1) {
      runtime <- proc.time()
      for (i in 1:B) {
        ida <- sample( 1:n, n, replace = TRUE )
        yb <- y[ida, ]
        xb <- x[ida, ]
        ini <- as.vector( t( kl.compreg(yb, xb[, -1, drop = FALSE])$be ) ) ## initial values
        qa <- nlm(jsreg, ini, y = yb, x = xb, d = d)
        qa <- nlm(jsreg, qa$estimate, y = yb, x = xb, d = d)
        qa <- nlm(jsreg, qa$estimate, y = yb, x = xb, d = d)
        betaboot[i, ] <- qa$estimate
      }
      s <- Rfast::colVars(betaboot, std = TRUE)
      seb <- matrix(s, byrow = TRUE, ncol = d)
      runtime <- proc.time() - runtime

    } else {
      runtime <- proc.time()
      cl <- makePSOCKcluster(ncores)
      registerDoParallel(cl)
      ww <- foreach::foreach(i = 1:B, .combine = rbind, .export=c("jsreg", "kl.compreg") ) %dopar% {
        ida <- sample(1:n, n, replace = TRUE)
        yb <- y[ida, ]
        xb <- x[ida, ]
        ini <- as.vector( t( kl.compreg(yb, xb[, -1, drop = FALSE])$be ) ) ## initial values
        qa <- nlm(jsreg, ini, y = yb, x = xb, d = d)
        qa <- nlm(jsreg, qa$estimate, y = yb, x = xb, d = d)
        qa <- nlm(jsreg, qa$estimate, y = yb, x = xb, d = d)
        betaboot <- qa$estimate
      }
      stopCluster(cl)
      s <- Rfast::colVars(ww, std = TRUE)
      seb <- matrix(s, byrow = TRUE, ncol = d)
      runtime <- proc.time() - runtime
    }
  }

  if ( !is.null(xnew) ) {
    xnew <- model.matrix(~., data.frame(xnew) )
    mu <- cbind( 1, exp(xnew %*% be) )
    est <- mu / Rfast::rowsums(mu)
  }  else  est <- NULL

  rownames(be)  <- colnames(x)
  if  ( !is.null(seb) ) rownames(seb) <- colnames(x)
  list(runtime = runtime, be = be, seb = seb, 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 June 4, 2018, 5:04 p.m.