https://github.com/cran/emplik
Raw File
Tip revision: 31d32093eee1aa45b156c2a215f600da3894af7f authored by Mai Zhou on 26 December 2006, 00:00:00 UTC
version 0.9-3-1
Tip revision: 31d3209
el.cen.EM2.R
el.cen.EM2 <- function(x,d,xc=1:length(x),fun,mu,maxit=25,error=1e-9,...){
####
#### xc is collaps control: if for index i and j, (x[] d[]) are identical
####     should they be merged into one obs. with weight 2?
####     if xc[i] != xc[j] then they do not merge. 
####     For regression applications, you should not merge -- use
####     xc = 1:length(x) .
####     The order of censoring (re-distribute weights)
####     is based on x, not fun(x).
####
   xvec <- as.vector(x) 
   d <- as.vector(d) 
   mu <- as.vector(mu)
   xc <- as.vector(xc)
   n <- length(d) 

   if (length(xvec)!=n) stop ("length of d and x must agree")
   if (length(xc)!=n) stop ("length of xc and d must agree")
   if(n <= 2*length(mu)+1) stop ("Need more observations")
   if(any((d!=0)&(d!=1)&(d!=2)))
     stop("d must be 0(right-censored) or 1(uncensored) or 2(left-censored)")
   if(!is.numeric(xvec)) stop("x must be numeric")
   if(!is.numeric(mu)) stop("mu must be numeric")

   funx <- as.matrix(fun(xvec, ...))   ##get the matrix before sorting xvec
   pp <- ncol(funx)
   if(length(mu)!=pp) stop("length of mu and ncol of fun(x) must agree")

   temp <- Wdataclean5(z=xvec,d,zc=xc, xmat=funx)  ## collaps control zc
   x <- temp$value
   d <- temp$dd
   w <- temp$weight
   funx <- temp$xxmat     ## make sure funx and x are in same order

   ###### change the last obs. among d=1,0, so that d=1
   ###### change the first obs. among d=1,2 so that d=1
   ###### this ensures we got a proper CDF for NPMLE. (no mass escapes)
   INDEX10 <- which(d != 2)
   d[ INDEX10[length(INDEX10)] ] <- 1
   INDEX12 <- which(d != 0)
   d[ INDEX12[1] ] <- 1

   xd1 <- x[d==1]
    if(length(xd1) <= 1) stop("need more distinct uncensored obs.")

   funxd1 <- funx[d==1,] 
   xd0 <- x[d==0]
   xd2 <- x[d==2]
   wd1 <- w[d==1]
   wd0 <- w[d==0]
   wd2 <- w[d==2]
   m <- length(xd0)
   mleft <- length(xd2)

##############################################
#### do the computation in 4 different cases.#
##############################################
 if( (m>0) && (mleft>0) ) { 
   pnew <- el.test.wt2(x=funxd1, wt=wd1, mu=mu)$prob
   n <- length(pnew)
   k <- rep(NA, m)
   for(i in 1:m) { k[i] <- 1+n - sum( xd1 > xd0[i] ) }
   kk <- rep(NA, mleft)
   for(j in 1:mleft) { kk[j] <- sum( xd1 < xd2[j] ) }
   num <- 1
   while(num < maxit) {
     wd1new <- wd1
     sur <- rev(cumsum(rev(pnew)))
     cdf <- 1 - c(sur[-1],0)
     for(i in 1:m)
        {wd1new[k[i]:n] <- wd1new[k[i]:n] + wd0[i]*pnew[k[i]:n]/sur[k[i]]}
     for(j in 1:mleft) {
     wd1new[1:kk[j]] <- wd1new[1:kk[j]] + wd2[j]*pnew[1:kk[j]]/cdf[kk[j]]}
     pnew <- el.test.wt2(x=funxd1, wt=wd1new, mu=mu)$prob
     num <- num +1
     }
   logel <- sum(wd1*log(pnew)) + sum(wd0*log(sur[k])) + sum(wd2*log(cdf[kk]))
   logel00 <- WCY(x=x,d=d,wt=w)$logEL   ####  should I add zc?
   }
  if( (m>0) && (mleft==0) ) {
   pnew <- el.test.wt2(x=funxd1, wt=wd1, mu=mu)$prob
   n <- length(pnew)
   k <- rep(NA,m)
   for(i in 1:m) { k[i] <- 1 + n - sum( xd1 > xd0[i] ) }
   num <- 1
   while(num < maxit) {
     wd1new <- wd1
     sur <- rev(cumsum(rev(pnew)))
     for(i in 1:m)
        {wd1new[k[i]:n] <- wd1new[k[i]:n] + wd0[i]*pnew[k[i]:n]/sur[k[i]]}
     pnew <- el.test.wt2(x=funxd1, wt=wd1new, mu=mu)$prob
     num <- num+1
     }
   sur <- rev(cumsum(rev(pnew)))
   logel <- sum( wd1*log(pnew)) + sum( wd0*log(sur[k]) )
   logel00 <- WKM(x=x,d=d,zc=xc,w=w)$logel
   }
  if( (m==0) && (mleft>0) ) {
   kk <- rep(NA, mleft)
   for(j in 1:mleft) { kk[j] <- sum( xd1 < xd2[j] ) }
   pnew <- el.test.wt2(x=funxd1, wt=wd1, mu=mu)$prob
   n <- length(pnew)
   num <- 1
   while(num < maxit) {
     wd1new <- wd1
     cdf <- cumsum(pnew) 
     for(j in 1:mleft)
       {wd1new[1:kk[j]] <- wd1new[1:kk[j]] + wd2[j]*pnew[1:kk[j]]/cdf[kk[j]]}
     pnew <- el.test.wt2(x=funxd1, wt=wd1new, mu=mu)$prob
     num <- num+1
     }
   logel <- sum( wd1*log(pnew)) + sum( wd2*log( cdf[kk] ) )
   dleft <- d 
   dleft[dleft==2] <- 0
   templeft <- WKM(x=rev(x),d=rev(dleft),zc=rev(xc),w=rev(w))
   logel00 <- templeft$logel   ### ?? do I need a left WKM ??
   }
  if( (m==0) && (mleft==0) ) { 
    num <- 0
    pnew <- el.test.wt2(x=funxd1, wt=wd1, mu)$prob
    logel <- sum( wd1*log(pnew) ) 
    logel00 <- sum( wd1*log( wd1/sum(wd1) ) )
  }
  tval <- 2*(logel00 - logel)
  list(loglik=logel, times=xd1, prob=pnew, iters=num, 
             "-2LLR"=tval, Pval= 1-pchisq(tval, df=length(mu)) )
}
back to top