kmc.solve: Calculate NPMLE with constriants for right censored data

In yfyang86/kmc: Kaplan-Meier Estimator with Constraints for Right Censored Data -- a Recursive Computational Algorithm

Calculate NPMLE with constriants for right censored data

Description

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

El(F)=∏_{i=1}^{n}(Δ F(T_i))^{δ_i}(1-F(T_i))^{1-δ_i}

with constraints

∑_i g(T_i)Δ 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	Δ 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)

yfyang86/kmc documentation built on Nov. 29, 2022, 1:27 p.m.