R/el.cen.test2.R

In emplik: Empirical Likelihood Ratio for Censored/Truncated Data

el.cen.test2 <- function(x,d,fun=function(x){cbind(x, pmin(x,3))}, mu, error = 1e-09, maxit = 15) 
{
    xvec <- as.vector(x)
    n <- length(xvec)
    kk <- length(mu) 
    if (n <= 2*kk +1) 
        stop("Need more observations")
    if (length(d) != n) 
        stop("length of x and d must agree")
    if (any((d != 0) & (d != 1))) 
        stop("d must be 0(right-censored) or 1(uncensored)")
    if (!is.numeric(xvec)) 
        stop("x must be numeric")
    if (!is.numeric(mu)) 
        stop("mu must be numeric")

       
    temp <- Wdataclean2(xvec, d)
    dd <- temp$dd
    dd[length(dd)] <- 1
    if (all(dd == 1)) 
        stop("there is no censoring, please use el.test()")
    xx <- temp$value
    n <- length(xx)
    ww <- temp$weight
    w0 <- WKM(x = xx, d = dd, w = ww)$jump
    uncenw0 <- w0[dd == 1]
    funxx <- as.matrix(fun(xx))
    nxkk <- dim(funxx)
      if( kk != nxkk[2] ) stop("mu and col(fun) not match")
      if( n != nxkk[1] ) stop("check dim of row(fun) value")

    for(i in 1:kk)  
            { if ((mu[i] > max(funxx[,i])) | (mu[i] < min(funxx[,i]))) 
        stop("check the value of mu/fun")
                     }

    xbar <- as.vector( uncenw0 %*% funxx[dd == 1,] )  # xbar is of length kk


    dvec01 <- uncenw0
    rk <- 1:n
    cenrk <- rk[dd == 0]
    mm <- length(cenrk)
    dvec02 <- rep(0, mm)
    for (j in 1:mm) dvec02[j] <- sum(w0[cenrk[j]:n])
    dvec00 <- rep(0, n)
    dvec00[dd == 1] <- dvec01
    dvec00[dd == 0] <- dvec02
    dvec0 <- ww/dvec00
    Dmat0 <- dvec00/sqrt(ww)
    mat <- matrix(rep(dd, mm), ncol = mm, nrow = n)
    for (i in 1:mm) {
        mat[1:cenrk[i], i] <- 0
        mat[cenrk[i], i] <- -1
    }

    for(i in 1:kk) { funxx[,i] <- funxx[,i]*dd }  ## need improve??

    Amat <- as.matrix(cbind(dd, funxx, mat))
    bvec0 <- c(0, as.vector(mu - xbar), rep(0, mm))
    value0 <- solve3.QP(Dmat0, dvec0, Amat, bvec0, meq =mm+1+kk, factorized = TRUE)
    w <- dvec00 + value0$solution
    
    if (any(w <= 0)) 
        stop("There is no probability satisfying the constraints")
 
    #**********end initial calculation **********************
    #**********Get ready to begin iteration ******************************
    # update vector bvec after initial calculation    
	
	bvec <- rep(0, mm+1+kk)
    diff <- 10
    m <- 0
    while ((diff > error) & (m < maxit)) {
        dvec <- ww/w
        Dmat <- w/sqrt(ww)
        value0 <- solve3.QP(Dmat, dvec, Amat, bvec, meq = mm+1+kk, factorized = TRUE)
        w <- w + value0$solution
        diff <- sum(abs(value0$solution))
        m <- m + 1
    }
	#**********end iteration *****************************
	
    lik00 <- sum(ww * log(dvec00))
    tval <- 2 * (lik00 - sum(ww * log(w)))
    list(llik00=lik00, llik11=sum(ww*log(w)), `-2LLR` = tval, 
         Pval = 1 - pchisq(tval, df = kk), prob1 = w[dd == 1], 
         xtime = xx[dd == 1], iteration = m, error = diff)
}
##### This function is very similar to the 1-dim cousin el.cen.test() which is older.