R/hMave.r

Defines functions hMave

Documented in hMave

#' @title Hazard Mave for Censored Survival Data
#' @name hMave
#' @description This is an almost direct R translation of Xia, Zhang & Xu's (2010) hMave Matlab code. We implemented further options for setting a different initial value. The computational algorithm does not utilize the orthogonality constrained optimization.
#' @param x A matrix for features.
#' @param y A vector of observed time.
#' @param censor A vector of censoring indicator.
#' @param m0 number of dimensions to use
#' @param B0 initial value of B. This is a feature we implemented.
#' @return A list consisting of
#' \item{B}{The estimated B matrix}
#' \item{cv}{Leave one out cross-validation error}
#' @references Xia, Y., Zhang, D., & Xu, J. (2010). Dimension reduction and semiparametric estimation of survival models. Journal of the American Statistical Association, 105(489), 278-290.
#' \url{http://dx.doi.org/10.1198/jasa.2009.tm09372}.
#' @examples
#' # generate some survival data
#' set.seed(1)
#' P = 7
#' N = 150
#' dataX = matrix(runif(N*P), N, P)
#' failEDR = as.matrix(cbind(c(1, 1.3, -1.3, 1, -0.5, 0.5, -0.5, rep(0, P-7))))
#' T = exp(dataX %*% failEDR + rnorm(N))
#' C = runif(N, 0, 15)
#' Y = pmin(T, C)
#' Censor = (T < C)
#'
#' # fit the model
#' hMave.fit = hMave(dataX, Y, Censor, 1)

hMave <- function(x, y, censor, m0, B0 = NULL)
{
  # this function returns B and cv

  # Input : x (n x p); y (n x 1); m0 (integer)
  # Output: B (p x m0)
  # m0: dimension e.g. it is 1 for single-index model

  if (!is.matrix(x)) stop("x must be a matrix")
  if (!is.matrix(y)) stop("y must be a matrix")
  if (!is.matrix(censor)) stop("censor must be a matrix")

  censor = as.matrix(as.numeric(censor))
  if (!all(censor %in% c(0,1)))
    stop("censor can only have two levels 0/1")

  n = nrow(x)
  p = ncol(x)

  onen = matrix(1, n, 1)
  x = scale(x, scale = FALSE)

  xeigen = eigen(t(x) %*% x /n)

  ss = xeigen$vectors[, rev(1:p), drop = FALSE]
  D = diag(rev(xeigen$values))

  ss = (ss/repmat(matrix(sqrt(diag(D)), 1, p), p, 1))

  x = x %*% ss

  # if initial value is given
  if (!is.null(B0))
  {
    if (!is.matrix(B0)) stop("B0 must be a matrix")
    if (nrow(B0) != ncol(x)) stop("Dimension of B0 is not correct")
    B0 = solve(ss) %*% B0
    B = B0
    m = ncol(B)
  }else{
    m = p
    B = diag(p)
  }

  Ba = B
  BI = Ba
  B = B[, 1:m, drop = FALSE]

  noip0 = 0
  noip1 = 1
  iterstop = 0
  Btmp = B

  rige = sd(y)*mean(apply(x, 2, sd))
  y = scale(y)

  #x = x./repmat(std(x,[],1),n,1);

  I = order(y)
  y = y[I, , drop = FALSE]
  x = x[I, , drop = FALSE]
  censor = censor[I, ,drop = FALSE]
  censor0 = censor

  niter = floor(p*3/2)
  #ch = (sqrt(p)/n^(1/(p+4))*n^(2/(m0+4)))^(1/niter);
  yc = y

  ky2 = repmat(yc, 1, n)
  ky2 = (ky2-t(ky2))^2
  censor = repmat(censor, 1, n)


  for (iter in 1:niter)
  {
    adj = p^(1/iter)*n^(1/(m0+4)-1/(m+4))
    hy = sd(yc)/n^(1/5)*p^(1/(iter+1))
    xB = x %*% Ba
    h = p^(1/(iter+1))*mean(apply(xB, 2, sd))/n^(1/(m+4))
    h2 = 2*h*h*adj

    H = matrix(1, 1, n)

    for (i in 1:n)
    {
      xBi = (xB-repmat(xB[i, ,drop = FALSE], n, 1))
      k = exp(-rowSums(xBi^2)/h2)
      kH = pnorm( (y[i] - y)/hy )
      H[, i] = t(k) %*% kH/sum(k)
    }

    ky = censor*exp(-ky2/(2*hy*hy))/hy/repmat(H,n,1)
    n1 = ncol(ky)

    ABI = matrix(0, m, n*n)

    for (iiter in 1:max(1, (m < p)*floor(m/2)))
    {
      dd = matrix(0, m*p, m*p)
      dc = matrix(0, m*p, 1)

      for (j in 1:n)
      {
        xij = x - repmat(x[j, , drop = FALSE], n, 1)
        sxij = rowSums((xij %*% Ba)^2) + rowSums(xij^2)/1.5^iter
        ker = exp(-sxij/h2)
        rker = repmat(as.matrix(ker), 1, p+1)
        onexi = cbind(xij %*% B, onen)
        xk = t(onexi*rker[, 1:(m+1)])
        abi = solve(xk %*% onexi + diag(1, m+1)/n) %*% (xk %*% ky)

        kxij = t(xij*rker[ , 1:p, drop = FALSE])
        kxijy = kxij %*% (ky - repmat(abi[m+1, , drop = FALSE],n,1))
        ddx = kxij %*% xij

        for (k1 in 1:m)
        {
          ka = ((k1-1)*p+1) : (k1*p)
          dc[ka,] = dc[ka,] + kxijy %*% t(abi[k1,,drop = FALSE])
        }

        tmp = abi[1:m, , drop = FALSE] %*% t(abi[1:m, , drop = FALSE])

        dd = dd + kronecker(tmp, ddx)

        ABI[, ((j-1)*n1+1):(j*n1)] = abi[1:m, , drop = FALSE]
      }

      B = solve(dd + rige* diag(1, length(dc))/n ) %*% dc

      B0 = matrix(B, p, m) %*% ABI

      B = eigen(B0 %*% t(B0))$vectors[, rev(1:p), drop = FALSE]

      B = B[, (p-m+1):p, drop = FALSE]
      Ba = B

      if (max(svd(B %*% t(B) - BI %*% t(BI))$d) < 0.001)
        break

      BI = B
    }

    mb = m
    ma = max(m-1,m0)
    m = ma
    B = B[ , (mb-m+1):mb, drop = FALSE]
    Ba = B

    if (max(svd(B %*% t(B) - BI %*% t(BI))$d) < 0.001 &  iter > p+3 )
      break

    BI = Ba

    if (iter == 1)
      B00 = B[, (ncol(B)-m0+1):ncol(B), drop = FALSE]

  }

  h2 = 2*n^(-2/(m0+3))

  x = x[censor0 == 1, , drop = FALSE]
  y = y[censor0 == 1, , drop = FALSE]
  n = nrow(y)
  ky = scale(y)
  ky = exp(-(repmat(ky,1,n) - repmat(t(ky),n,1))^2/h2)

  kye = ky

  cv = 0

  for (j in 1:n)
  {
    xij = x - repmat(x[j, , drop = FALSE], n, 1)
    sxij = rowSums((xij %*% Ba)^2)
    ker = as.matrix(exp(-sxij/h2))
    ker[j] = 0
    kye[j, ] = t(ker) %*% ky / sum(ker)
  }

  cv = sum((ky-kye)^2)
  B = ss %*% B

  apply(B, 2, function(x) x / sqrt(sum(x^2)))

  return(list("B" = B, "cv" = cv))
}

Try the orthoDr package in your browser

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

orthoDr documentation built on Sept. 5, 2019, 5:03 p.m.