kmc.solve {kmc}

Calculate NPMLE with constriants for right censored data

Description

This function calculate the Kaplan-Meier estimator with mean constraints recursively.

El(F)=\prod_{i=1}^{n}(\Delta F(T_i))^{\delta_i}(1-F(T_i))^{1-\delta_i}

with constraints

\sum_i g(T_i)\Delta F(T_i)=0,\quad,i=1,2,...

It uses Lagrange multiplier directly.

Usage

kmc.solve(x, d, g, em.boost = T, using.num = T, using.Fortran =
                 T, using.C = F, tmp.tag = T, rtol = 1e-09, control =
                 list(nr.it = 20, nr.c = 1, em.it = 3),...)

Arguments

x

Non-negative real vector. The observed time.

d

0/1 vector. Censoring status indictator, 0: right censored; 1 uncensored

g

list of contraint functions. It should be a list of functions list(f1,f2,...)

em.boost

A logical value. It determines whether the EM algorithm is used to get the initial value, default=TRUE. See 'Details' for EM control.

using.num

A logical value. It determines whether the numeric derivatives is used in iterations, default=TRUE.

using.Fortran

A logical value. It determines whether Fortran is used in root solving, default=F.

using.C

A logical value. It determines whether to use Rcpp in iteraruib, default=T. This option will promote the computational efficiency of the KMC algorithm. Development version works on one constraint only, otherwise it will generate a Error information. It won't work on using.num=F.

tmp.tag

Development version needs it, keep it as TRUE.

rtol

Tolerance used in rootSolve(multiroot) package, see 'rootSolve::multiroot'.

control

A list. The entry nr.it controls max iterations allowed in N-R algorithm default=20; nr.c is the scaler used in N-R algorithm default=1; em.it is max iteration if use EM algorithm (em.boost) to get the initial value of lambda, default=3.

...

Unspecified yet.

Details

The function check_G_mat checks whether the solution space is null or not under the constraint. But due to the computational complexity, it will detect at most two conditions.

Value

A list with the following components:

loglik.ha

The log empirical likelihood without constraints

loglik.h0

The log empirical likelihood with constraints

"-2LLR"

The -2 Log empirical likelihood ratio

phat

\Delta F(T_i)

pvalue

The p-value of the test

df

Degree(s) of freedom. It equals the number of constraints.

lambda

The lambda is the Lagrangian multiplier described in reference.

Author(s)

Mai Zhou([e-mail]), Yifan Yang([e-mail])

References

Zhou, M. and Yang, Y. (2015). A recursive formula for the Kaplan-Meier estimator with mean constraints and its application to empirical likelihood Computational Statistics Online ISSN 1613-9658.

See Also

plotkmc2D.

Examples

# positive time
x <- c( 1, 1.5, 2, 3, 4.2, 5.0, 6.1, 5.3, 4.5, 0.9, 2.1, 4.3) 
# status censored/uncensored
d <- c( 1, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, 1)               

#################
# dim =1
#################

f <- function(x) {x-3.7}                     # \sum f(ti) wi ~ 0 
g <- list(f=f) ;                             #define constraint as a list

kmc.solve(x, d, g) ;                         #using default
kmc.solve(x, d, g, using.C=TRUE) ;           #using Rcpp

#################
# dim =2
#################

myfun5 <- function( x)  { 
 x^2-16.5
} 

g <- list( f1=f,f2=myfun5) ;                 #define constraint as a list

re0 <- kmc.solve(x,d,g);

###################################################
# Print Estimation and other information 
# with option: digits = 5
###################################################

#' Print kmc object
#' 
#' @param x kmc object
#' @param digits minimal number of significant digits, see print.default.
f_print <- function(x, digits = 5){
  cat("\n----------------------------------------------------------------\n")
  cat("A Recursive Formula for the Kaplan-Meier Estimator with Constraint\n")
  cat("Information:\n")
  cat("Number of Constraints:\t", length(x$g), "\n")
  cat("lamda(s):\t", x$lambda,'\n');
  cat("\n----------------------------------------------------------------\n")
  names <- c("Log-likelihood(Ha)", "Log-likelihood(H0)",
  "-2LLR", paste("p-Value(df=", length(x$g), ")",sep = ""))
  re <- matrix(c(x[[1]], x[[2]], x[[3]], 1 - pchisq(x[[3]],
  length(x$g))), nrow = 1)
  colnames(re) <- names
  rownames(re) <- "Est"
  print.default(format(re, digits = digits), print.gap = 2,
                quote = FALSE, df = length(x$g))
  cat("------------------------------------------------------------------\n")
}

f_print(re0)