R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("emplik-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('emplik')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "BJnoint"
> 
> ### * BJnoint
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: BJnoint
> ### Title: The Buckley-James censored regression estimator
> ### Aliases: BJnoint
> ### Keywords: nonparametric htest
> 
> ### ** Examples
> 
> x <- matrix(c(rnorm(50,mean=1), rnorm(50,mean=2)), ncol=2,nrow=50)
> ## Suppose now we wish to test Ho: 2mu(1)-mu(2)=0, then
> y <- 2*x[,1]-x[,2]
> xx <- c(28,-44,29,30,26,27,22,23,33,16,24,29,24,40,21,31,34,-2,25,19)
> 
> 
> 
> cleanEx(); ..nameEx <- "bjtest"
> 
> ### * bjtest
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bjtest
> ### Title: Test the Buckley-James estimator by Empirical Likelihood
> ### Aliases: bjtest
> ### Keywords: nonparametric htest
> 
> ### ** Examples
> 
> xx <- c(28,-44,29,30,26,27,22,23,33,16,24,29,24,40,21,31,34,-2,25,19)
> 
> 
> 
> cleanEx(); ..nameEx <- "bjtest1d"
> 
> ### * bjtest1d
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bjtest1d
> ### Title: Test the Buckley-James estimator by Empirical Likelihood, 1-dim
> ###   only
> ### Aliases: bjtest1d
> ### Keywords: nonparametric htest
> 
> ### ** Examples
> 
> xx <- c(28,-44,29,30,26,27,22,23,33,16,24,29,24,40,21,31,34,-2,25,19)
> 
> 
> 
> cleanEx(); ..nameEx <- "el.cen.EM"
> 
> ### * el.cen.EM
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.cen.EM
> ### Title: Empirical likelihood ratio for mean with right, left or doubly
> ###   censored data, by EM algorithm
> ### Aliases: el.cen.EM
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> ## example with tied observations
> x <- c(1, 1.5, 2, 3, 4, 5, 6, 5, 4, 1, 2, 4.5)
> d <- c(1,   1, 0, 1, 0, 1, 1, 1, 1, 0, 0,   1)
> el.cen.EM(x,d,mu=3.5)
$loglik
[1] -17.04924

$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.1489836 0.1516380 0.1447664 0.1108301 0.1277623 0.2251085 0.0909112

$funMLE
[1] 4.065923

$"-2LLR"
[1] 1.246634

$Pval
[1] 0.2641962

> ## we should get "-2LLR" = 1.2466....
> myfun5 <- function(x, theta, eps) {
+ u <- (x-theta)*sqrt(5)/eps 
+ INDE <- (u < sqrt(5)) & (u > -sqrt(5)) 
+ u[u >= sqrt(5)] <- 0 
+ u[u <= -sqrt(5)] <- 1 
+ y <- 0.5 - (u - (u)^3/15)*3/(4*sqrt(5)) 
+ u[ INDE ] <- y[ INDE ] 
+ return(u)
+ }
> el.cen.EM(x, d, fun=myfun5, mu=0.5, theta=3.5, eps=0.1)
$loglik
[1] -17.26379

$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.12291732 0.14295095 0.23413173 0.08333333 0.10416667 0.20833333 0.10416667

$funMLE
[1] 0.2928571

$"-2LLR"
[1] 1.675727

$Pval
[1] 0.1954932

> 
> 
> 
> cleanEx(); ..nameEx <- "el.cen.EM2"
> 
> ### * el.cen.EM2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.cen.EM2
> ### Title: Empirical likelihood ratio test for a vector of means with
> ###   right, left or doubly censored data, by EM algorithm
> ### Aliases: el.cen.EM2
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> ## censored regression with one right censored observation.
> ## we check the estimation equation, with the MLE inside myfun7. 
> y <- c(3, 5.3, 6.4, 9.1, 14.1, 15.4, 18.1, 15.3, 14, 5.8, 7.3, 14.4)
> x <- c(1, 1.5, 2,   3,   4,    5,    6,    5,    4,  1,   2,   4.5)
> d <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0)
> myfun7 <- function(y, xmat) {
+ temp1 <- y - ( 1.392885 +  2.845658 * xmat)
+ return( cbind( temp1, xmat*temp1) )
+ }
> el.cen.EM2(y,d, xc=1:12, fun=myfun7, mu=c(0,0), xmat=x)
$loglik
[1] -27.56954

$times
 [1]  3.0  5.3  5.8  6.4  7.3  9.1 14.0 14.1 15.3 15.4 18.1

$prob
 [1] 0.08333336 0.08333334 0.08333330 0.08333335 0.08333333 0.08333336
 [7] 0.08333329 0.08333329 0.11111113 0.11111112 0.11111113

$iters
[1] 25

$"-2LLR"
[1] 1.037392e-12

$Pval
[1] 1

> ## we should get, Pval = 1 , as the MLE should.
> ## for different values of (a, b) inside myfun7, you get other Pval
> rqfun1 <- function(y, xmat, beta, tau = 0.5) {
+ temp1 <- tau - (1-myfun55(y-beta*xmat))
+ return(xmat * temp1)
+ }
> myfun55 <- function(x, eps=0.001){
+ u <- x*sqrt(5)/eps
+ INDE <- (u < sqrt(5)) & (u > -sqrt(5))
+ u[u >= sqrt(5)] <- 0
+ u[u <= -sqrt(5)] <- 1
+ y <- 0.5 - (u - (u)^3/15)*3/(4*sqrt(5))
+ u[ INDE ] <- y[ INDE ]
+ return(u)
+ }
> el.cen.EM2(x=y,d=d,xc=1:12,fun=rqfun1,mu=0,xmat=x,beta=3.08,tau=0.44769875)
$loglik
[1] -27.56954

$times
 [1]  3.0  5.3  5.8  6.4  7.3  9.1 14.0 14.1 15.3 15.4 18.1

$prob
 [1] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333
 [7] 0.08333333 0.08333333 0.11111111 0.11111111 0.11111111

$iters
[1] 25

$"-2LLR"
[1] -7.105427e-15

$Pval
[1] 1

> ## default tau=0.5 
> el.cen.EM2(x=y,d=d,xc=1:12,fun=rqfun1,mu=0,xmat=x,beta=3.0799107404)
$loglik
[1] -27.56954

$times
 [1]  3.0  5.3  5.8  6.4  7.3  9.1 14.0 14.1 15.3 15.4 18.1

$prob
 [1] 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333 0.08333333
 [7] 0.08333334 0.08333334 0.11111110 0.11111112 0.11111110

$iters
[1] 25

$"-2LLR"
[1] 7.105427e-15

$Pval
[1] 1

> 
> 
> 
> cleanEx(); ..nameEx <- "el.cen.test"
> 
> ### * el.cen.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.cen.test
> ### Title: Empirical likelihood ratio for mean with right censored data, by
> ###   QP.
> ### Aliases: el.cen.test
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> el.cen.test(rexp(100), c(rep(0,25),rep(1,75)), mu=1.5)
$"-2LLR"
[1] 1.417819

$Pval
[1] 0.2337627

$weights
 [1] 0.009019268 0.009028218 0.009032587 0.009032695 0.009051070 0.009164393
 [7] 0.009486021 0.009490275 0.009490770 0.009507151 0.009515380 0.009524518
[13] 0.009526481 0.009541247 0.009657965 0.009658282 0.009663630 0.009673750
[19] 0.009805447 0.009807579 0.009832183 0.009834047 0.009858158 0.009864585
[25] 0.010003775 0.010004404 0.010054296 0.010223778 0.010226256 0.010390030
[31] 0.010557269 0.010962220 0.011001654 0.011025395 0.011423997 0.011453160
[37] 0.011462136 0.011463276 0.011467164 0.011485269 0.011854772 0.011872727
[43] 0.011875093 0.011887101 0.011902320 0.011908137 0.011915362 0.012206705
[49] 0.012552455 0.012583048 0.012585686 0.012619518 0.012662855 0.013517308
[55] 0.013987469 0.014044652 0.014045917 0.014086134 0.014123368 0.014851512
[61] 0.015067699 0.015079092 0.015098910 0.015447541 0.016801730 0.016998212
[67] 0.017168285 0.017401661 0.017466323 0.020399336 0.020488989 0.024270060
[73] 0.025697606 0.029936026 0.074326613

$xtime
 [1] 0.03726853 0.05205545 0.05926121 0.05943916 0.08967408 0.13257141
 [7] 0.20351035 0.20986658 0.21060689 0.23502745 0.24726425 0.26082824
[13] 0.26373826 0.28559098 0.30128300 0.30174093 0.30944786 0.32401015
[19] 0.34888835 0.35187050 0.38619356 0.38878677 0.42224244 0.43113213
[25] 0.44660565 0.44745189 0.51417430 0.55464140 0.55782935 0.57871246
[31] 0.59461765 0.68122911 0.72521430 0.75154269 0.77418776 0.80417091
[37] 0.81336825 0.81453581 0.81851419 0.83700649 0.97739581 0.99455579
[43] 0.99681296 1.00825646 1.02272588 1.02824690 1.03509715 1.04360851
[49] 1.07988114 1.10593627 1.10817666 1.13683142 1.17331211 1.23379151
[55] 1.25310535 1.29226165 1.29312466 1.32046793 1.34564402 1.43528534
[61] 1.56524135 1.57198685 1.58369608 1.78476540 2.00783240 2.10037723
[67] 2.17877257 2.28385347 2.31247163 2.73038931 2.75924379 2.90988727
[73] 3.21778902 3.95893285 4.83281274

$iteration
[1] 4

$error
[1] 3.601390e-12

> ## second example with tied observations
> x <- c(1, 1.5, 2, 3, 4, 5, 6, 5, 4, 1, 2, 4.5)
> d <- c(1,   1, 0, 1, 0, 1, 1, 1, 1, 0, 0,   1)
> el.cen.test(x,d,mu=3.5)
$"-2LLR"
[1] 1.246634

$Pval
[1] 0.2641962

$weights
[1] 0.1489836 0.1516380 0.1447664 0.1108301 0.1277623 0.2251085 0.0909112

$xtime
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$iteration
[1] 5

$error
[1] 9.16558e-11

> # we should get  "-2LLR" = 1.246634  etc. 
> 
> 
> 
> cleanEx(); ..nameEx <- "el.ltrc.EM"
> 
> ### * el.ltrc.EM
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.ltrc.EM
> ### Title: Empirical likelihood ratio for mean with left truncated and
> ###   right censored data, by EM algorithm
> ### Aliases: el.ltrc.EM
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> ## example with tied observations
> y <- c(0, 0, 0.5, 0, 1, 2, 2, 0, 0, 0, 0, 0 )
> x <- c(1, 1.5, 2, 3, 4, 5, 6, 5, 4, 1, 2, 4.5)
> d <- c(1,   1, 0, 1, 0, 1, 1, 1, 1, 0, 0,   1)
> el.ltrc.EM(y,x,d,mu=3.5)
$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.1667054 0.1515407 0.1190516 0.1041636 0.1246927 0.2320187 0.1018273

$"-2LLR"
[1] 0.4993391

$Pval
[1] 0.4797907

> ypsy <- c(51, 58, 55, 28, 25, 48, 47, 25, 31, 30, 33, 43, 45, 35, 36)
> xpsy <- c(52, 59, 57, 50, 57, 59, 61, 61, 62, 67, 68, 69, 69, 65, 76)
> dpsy <- c(1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1 )
> el.ltrc.EM(ypsy,xpsy,dpsy,mu=64)
$times
[1] 50 52 57 59 61 69 76

$prob
[1] 0.07070221 0.06705267 0.13065270 0.12218681 0.12451658 0.30359323 0.18129580

$"-2LLR"
[1] 0.1256373

$Pval
[1] 0.722999

> 
> 
> 
> cleanEx(); ..nameEx <- "el.test"
> 
> ### * el.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.test
> ### Title: Empirical likelihood ratio test for the means, uncensored data
> ### Aliases: el.test
> ### Keywords: nonparametric htest
> 
> ### ** Examples
> 
> x <- matrix(c(rnorm(50,mean=1), rnorm(50,mean=2)), ncol=2,nrow=50)
> el.test(x, mu=c(1,2))
$"-2LLR"
[1] 1.545630

$Pval
[1] 0.4617116

$lambda
[1] 0.1400026 0.1403307

$grad
[1] 1.734669e-14 1.989300e-15

$hess
          [,1]      [,2]
[1,] 44.718279 -1.234859
[2,] -1.234859 43.388479

$wts
 [1] 1.0328857 1.0640285 1.0742530 0.9390914 0.8017773 0.8598133 0.9835657
 [8] 1.0451033 0.8616530 1.0657657 0.6457137 0.9532230 0.9902791 1.4412007
[15] 0.9494970 0.9801975 1.3432876 0.7474949 0.8799105 0.7204446 0.8365486
[22] 0.9902227 0.9122898 1.6937581 1.0978709 0.9680161 1.0917260 1.2590558
[29] 1.0598945 1.0248159 0.9005626 1.0345116 0.8199431 1.2842510 1.1228927
[36] 1.0115081 0.9140902 1.0537301 0.8292346 0.8738748 1.1100824 0.8818112
[43] 0.7933884 0.8501986 0.8878993 1.0211150 1.1469240 0.9735678 1.2308897
[50] 0.9461412

$nits
[1] 5

> ## Suppose now we wish to test Ho: 2mu(1)-mu(2)=0, then
> y <- 2*x[,1]-x[,2]
> el.test(y, mu=0)
$"-2LLR"
[1] 0.09229383

$Pval
[1] 0.761281

$lambda
[1] 0.02199372

$grad
[1] 4.20111e-13

$hess
         [,1]
[1,] 192.8757

$wts
 [1] 1.0376802 0.9789154 1.0463093 0.9132327 1.0173181 1.0865392 0.9713290
 [8] 0.9474709 0.9873648 1.0105735 0.9865057 0.9823074 1.0443828 1.1086903
[15] 0.9382358 1.0061666 0.9624773 0.9908015 0.9682859 1.0221391 0.9709066
[22] 0.9523641 1.0102563 1.0717679 0.9480141 1.0089585 0.9971118 1.0691973
[29] 1.0231934 0.9696030 0.9325992 1.0015507 1.0089370 0.9697984 1.0794897
[36] 1.0262487 1.0424543 0.9959355 0.9613081 0.9730507 0.9953269 1.0391880
[43] 0.9948902 0.9909959 1.0697441 1.0453735 0.9577498 0.9556447 0.9784813
[50] 0.9531348

$nits
[1] 4

> xx <- c(28,-44,29,30,26,27,22,23,33,16,24,29,24,40,21,31,34,-2,25,19)
> el.test(xx, mu=15)  #### -2LLR = 1.805702
$"-2LLR"
[1] 1.805702

$Pval
[1] 0.1790247

$lambda
[1] 0.01078493

$grad
[1] 3.256062e-08

$hess
         [,1]
[1,] 28962.89

$wts
 [1] 0.8770359 2.7496028 0.8688180 0.8607526 0.8939472 0.8854108 0.9298048
 [8] 0.9205734 0.8374306 0.9893301 0.9115235 0.8688180 0.9115235 0.7876352
[15] 0.9392232 0.8528356 0.8299349 1.2245056 0.9026498 0.9586443

$nits
[1] 7

> 
> 
> 
> cleanEx(); ..nameEx <- "el.test.wt"
> 
> ### * el.test.wt
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.test.wt
> ### Title: Weighted Empirical Likelihood ratio for mean, uncensored data
> ### Aliases: el.test.wt
> ### Keywords: nonparametric htest
> 
> ### ** Examples
> 
> ## example with tied observations
> x <- c(1, 1.5, 2, 3, 4, 5, 6, 5, 4, 1, 2, 4.5)
> d <- c(1,   1, 0, 1, 0, 1, 1, 1, 1, 0, 0,   1)
> el.cen.EM(x,d,mu=3.5)
$loglik
[1] -17.04924

$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.1489836 0.1516380 0.1447664 0.1108301 0.1277623 0.2251085 0.0909112

$funMLE
[1] 4.065923

$"-2LLR"
[1] 1.246634

$Pval
[1] 0.2641962

> ## we should get "-2LLR" = 1.2466....
> myfun5 <- function(x, theta, eps) {
+ u <- (x-theta)*sqrt(5)/eps 
+ INDE <- (u < sqrt(5)) & (u > -sqrt(5)) 
+ u[u >= sqrt(5)] <- 0 
+ u[u <= -sqrt(5)] <- 1 
+ y <- 0.5 - (u - (u)^3/15)*3/(4*sqrt(5)) 
+ u[ INDE ] <- y[ INDE ] 
+ return(u)
+ }
> el.cen.EM(x, d, fun=myfun5, mu=0.5, theta=3.5, eps=0.1)
$loglik
[1] -17.26379

$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.12291732 0.14295095 0.23413173 0.08333333 0.10416667 0.20833333 0.10416667

$funMLE
[1] 0.2928571

$"-2LLR"
[1] 1.675727

$Pval
[1] 0.1954932

> 
> 
> 
> cleanEx(); ..nameEx <- "el.test.wt2"
> 
> ### * el.test.wt2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.test.wt2
> ### Title: Weighted Empirical Likelihood ratio for mean(s), uncensored data
> ### Aliases: el.test.wt2
> ### Keywords: nonparametric htest
> 
> ### ** Examples
> 
> ## example with tied observations
> x <- c(1, 1.5, 2, 3, 4, 5, 6, 5, 4, 1, 2, 4.5)
> d <- c(1,   1, 0, 1, 0, 1, 1, 1, 1, 0, 0,   1)
> el.cen.EM(x,d,mu=3.5)
$loglik
[1] -17.04924

$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.1489836 0.1516380 0.1447664 0.1108301 0.1277623 0.2251085 0.0909112

$funMLE
[1] 4.065923

$"-2LLR"
[1] 1.246634

$Pval
[1] 0.2641962

> ## we should get "-2LLR" = 1.2466....
> myfun5 <- function(x, theta, eps) {
+ u <- (x-theta)*sqrt(5)/eps 
+ INDE <- (u < sqrt(5)) & (u > -sqrt(5)) 
+ u[u >= sqrt(5)] <- 0 
+ u[u <= -sqrt(5)] <- 1 
+ y <- 0.5 - (u - (u)^3/15)*3/(4*sqrt(5)) 
+ u[ INDE ] <- y[ INDE ] 
+ return(u)
+ }
> el.cen.EM(x, d, fun=myfun5, mu=0.5, theta=3.5, eps=0.1)
$loglik
[1] -17.26379

$times
[1] 1.0 1.5 3.0 4.0 4.5 5.0 6.0

$prob
[1] 0.12291732 0.14295095 0.23413173 0.08333333 0.10416667 0.20833333 0.10416667

$funMLE
[1] 0.2928571

$"-2LLR"
[1] 1.675727

$Pval
[1] 0.1954932

> 
> 
> 
> cleanEx(); ..nameEx <- "el.trun.test"
> 
> ### * el.trun.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: el.trun.test
> ### Title: Empirical likelihood ratio for mean with left truncated data
> ### Aliases: el.trun.test
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> ## example with tied observations
> vet <- c(30, 384, 4, 54, 13, 123, 97, 153, 59, 117, 16, 151, 22, 56, 21, 18,
+          139, 20, 31, 52, 287, 18, 51, 122, 27, 54, 7, 63, 392, 10)
> vetstart <- c(0,60,0,0,0,33,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
> el.trun.test(vetstart, vet, mu=80, maxit=15)
$NPMLE
 [1] 0.03571429 0.03571429 0.03571429 0.03571429 0.03571429 0.07142857
 [7] 0.03571429 0.03571429 0.03571429 0.03571429 0.03571429 0.03571429
[13] 0.03348214 0.03348214 0.06696429 0.03348214 0.03348214 0.03043831
[19] 0.03043831 0.03043831 0.03043831 0.03043831 0.03043831 0.03043831
[25] 0.03043831 0.03043831 0.03043831 0.03043831

$NPMLEmu
 [1] 0.03606427 0.03605032 0.03603639 0.03602246 0.03600855 0.07199857
 [7] 0.03599002 0.03598539 0.03598076 0.03595763 0.03594377 0.03593915
[13] 0.03358417 0.03358014 0.06714415 0.03356402 0.03355194 0.03044728
[19] 0.03033500 0.03026935 0.03025297 0.03024970 0.03019745 0.03015838
[25] 0.03015187 0.02972252 0.02941927 0.02939453

$"-2LLR"
[1] 0.003933062

> 
> 
> 
> cleanEx(); ..nameEx <- "emplikH.disc"
> 
> ### * emplikH.disc
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: emplikH.disc
> ### Title: Empirical likelihood ratio for discrete hazard with right
> ###   censored, left truncated data
> ### Aliases: emplikH.disc
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> fun4 <- function(x, theta) { as.numeric(x <= theta) }
> x <- c(1, 2, 3, 4, 5, 6, 5, 4, 3, 4, 1, 2.4, 4.5)
> d <- c(1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1)
> # test if -H(4) = -0.7 
> emplikH.disc(x=x,d=d,K=-0.7,fun=fun4,theta=4)
$"discrete.-2LLR"
[1] 0.1446316

$lambda
[1] -0.09625977

$times
[1] 1.0 2.4 3.0 4.0 4.5 5.0

$jumps
[1] 0.08511636 0.11430370 0.25811038 0.17395470 0.25000000 0.66666667

> # we should get "discrete.-2logemlikRatio" 0.1446316  etc....
> y <- c(-2,-2, -2, 1.5, -1)
> emplikH.disc(x=x,d=d,y=y,K=-0.7,fun=fun4,theta=4)
$"discrete.-2LLR"
[1] 0.1258724

$lambda
[1] -0.0898586

$times
[1] 1.0 2.4 3.0 4.0 4.5 5.0

$jumps
[1] 0.09232043 0.11322671 0.25536789 0.17147252 0.25000000 0.66666667

> 
> 
> 
> cleanEx(); ..nameEx <- "emplikH.disc2"
> 
> ### * emplikH.disc2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: emplikH.disc2
> ### Title: Two sample empirical likelihood ratio for discrete hazards with
> ###   right censored, left truncated data, one parameter.
> ### Aliases: emplikH.disc2
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> library(boot)
> data(channing)
> ymale <- channing[1:97,2]
> dmale <- channing[1:97,5]
> xmale <- channing[1:97,3]
> yfemale <- channing[98:462,2]
> dfemale <- channing[98:462,5]
> xfemale <- channing[98:462,3]
> fun1 <- function(x) { as.numeric(x <= 960) }
> emplikH.disc2(x1=xfemale, d1=dfemale, y1=yfemale, x2=xmale,d2=dmale, y2=ymale, theta=0.2, fun1=fun1, fun2=fun1, maxi=4, mini=-10)
[1]  0.5925540 -0.2118239  0.5925535
$"-2LLR"
[1] 1.511239

$lambda
[1] -4.060575

$times1
  [1]  804  822  830  840  845  861  868  873  883  885  895  897  901  905  908
 [16]  911  912  915  919  923  926  928  930  931  932  934  936  940  941  944
 [31]  948  954  959  963  966  969  970  975  976  978  982  983  986  989  990
 [46]  991  992  994  995  996  998  999 1000 1001 1003 1004 1005 1006 1010 1011
 [61] 1012 1013 1014 1015 1017 1018 1019 1020 1021 1023 1024 1027 1029 1030 1033
 [76] 1040 1041 1043 1044 1047 1056 1063 1064 1068 1070 1072 1073 1074 1083 1084
 [91] 1085 1086 1089 1093 1097 1115 1122 1131 1132 1142 1152 1172 1192 1200

$times2
 [1]  777  869  872  876  893  894  898  907  909  911  927  932  945  948  957
[16]  966  969  971  983  985  989  993  998 1009 1012 1022 1025 1029 1031 1033
[31] 1036 1043 1044 1053 1055 1059 1060 1080 1085 1094 1128 1139

$wts1
  [1] 0.059033881 0.031309267 0.023843913 0.018539315 0.016144806 0.011773096
  [7] 0.010533032 0.009906932 0.009013928 0.008776593 0.007579236 0.007410733
 [13] 0.007145949 0.013894734 0.013991941 0.006852158 0.006852158 0.006947367
 [19] 0.006899434 0.006538536 0.006454135 0.006581570 0.006581570 0.006581570
 [25] 0.006581570 0.006669360 0.006291705 0.006371885 0.006454135 0.012427036
 [31] 0.006454135 0.006496062 0.006412746 0.006369427 0.006451613 0.012345679
 [37] 0.006250000 0.006329114 0.006451613 0.006535948 0.013071895 0.006622517
 [43] 0.006849315 0.006944444 0.028368794 0.007299270 0.007352941 0.014925373
 [49] 0.015151515 0.023255814 0.008000000 0.008130081 0.008196721 0.008333333
 [55] 0.008474576 0.008547009 0.008547009 0.017391304 0.009090909 0.009345794
 [61] 0.019230769 0.010101010 0.010000000 0.010416667 0.010526316 0.021276596
 [67] 0.021739130 0.011627907 0.012048193 0.024691358 0.012987013 0.012987013
 [73] 0.012820513 0.013157895 0.014084507 0.059701493 0.047619048 0.016666667
 [79] 0.017241379 0.017857143 0.041666667 0.023809524 0.025000000 0.052631579
 [85] 0.027777778 0.028571429 0.030303030 0.031250000 0.033333333 0.034482759
 [91] 0.071428571 0.040000000 0.043478261 0.047619048 0.047619048 0.066666667
 [97] 0.076923077 0.083333333 0.090909091 0.100000000 0.125000000 0.142857143
[103] 0.250000000 0.666666667

$wts2
 [1] 0.16500085 0.03563719 0.03441088 0.03441088 0.02698285 0.02698285
 [7] 0.02773112 0.02627391 0.02698285 0.02698285 0.02435426 0.02496220
[13] 0.02496220 0.02560126 0.02496220 0.02702703 0.02631579 0.02631579
[19] 0.02564103 0.02564103 0.02564103 0.05263158 0.02857143 0.03225806
[25] 0.06250000 0.03571429 0.03846154 0.04000000 0.04166667 0.04761905
[31] 0.05000000 0.04545455 0.05000000 0.05555556 0.05882353 0.06666667
[37] 0.07142857 0.09090909 0.10000000 0.25000000 0.33333333 0.50000000

> ######################################################
> ### You should get "-2LLR" = 1.511239 and a lot more other outputs.
> ########################################################
> emplikH.disc2(x1=xfemale, d1=dfemale, y1=yfemale, x2=xmale,d2=dmale, y2=ymale, theta=0.25, fun1=fun1, fun2=fun1, maxi=4, mini=-5)
[1]  0.54255402 -0.09407195  0.54255350
$"-2LLR"
[1] 1.150098

$lambda
[1] -3.195481

$times1
  [1]  804  822  830  840  845  861  868  873  883  885  895  897  901  905  908
 [16]  911  912  915  919  923  926  928  930  931  932  934  936  940  941  944
 [31]  948  954  959  963  966  969  970  975  976  978  982  983  986  989  990
 [46]  991  992  994  995  996  998  999 1000 1001 1003 1004 1005 1006 1010 1011
 [61] 1012 1013 1014 1015 1017 1018 1019 1020 1021 1023 1024 1027 1029 1030 1033
 [76] 1040 1041 1043 1044 1047 1056 1063 1064 1068 1070 1072 1073 1074 1083 1084
 [91] 1085 1086 1089 1093 1097 1115 1122 1131 1132 1142 1152 1172 1192 1200

$times2
 [1]  777  869  872  876  893  894  898  907  909  911  927  932  945  948  957
[16]  966  969  971  983  985  989  993  998 1009 1012 1022 1025 1029 1031 1033
[31] 1036 1043 1044 1053 1055 1059 1060 1080 1085 1094 1128 1139

$wts1
  [1] 0.056165516 0.030483605 0.023362019 0.018246670 0.015922421 0.011654398
  [7] 0.010437921 0.009822747 0.008944182 0.008710459 0.007529864 0.007363525
 [13] 0.007102045 0.013811724 0.013907769 0.006811779 0.006811779 0.006905862
 [19] 0.006858498 0.006501760 0.006418299 0.006544309 0.006544309 0.006544309
 [25] 0.006544309 0.006631101 0.006257645 0.006336954 0.006418299 0.012360594
 [31] 0.006418299 0.006459760 0.006377367 0.006369427 0.006451613 0.012345679
 [37] 0.006250000 0.006329114 0.006451613 0.006535948 0.013071895 0.006622517
 [43] 0.006849315 0.006944444 0.028368794 0.007299270 0.007352941 0.014925373
 [49] 0.015151515 0.023255814 0.008000000 0.008130081 0.008196721 0.008333333
 [55] 0.008474576 0.008547009 0.008547009 0.017391304 0.009090909 0.009345794
 [61] 0.019230769 0.010101010 0.010000000 0.010416667 0.010526316 0.021276596
 [67] 0.021739130 0.011627907 0.012048193 0.024691358 0.012987013 0.012987013
 [73] 0.012820513 0.013157895 0.014084507 0.059701493 0.047619048 0.016666667
 [79] 0.017241379 0.017857143 0.041666667 0.023809524 0.025000000 0.052631579
 [85] 0.027777778 0.028571429 0.030303030 0.031250000 0.033333333 0.034482759
 [91] 0.071428571 0.040000000 0.043478261 0.047619048 0.047619048 0.066666667
 [97] 0.076923077 0.083333333 0.090909091 0.100000000 0.125000000 0.142857143
[103] 0.250000000 0.666666667

$wts2
 [1] 0.19247497 0.03677082 0.03546668 0.03546668 0.02762776 0.02762776
 [7] 0.02841274 0.02688499 0.02762776 0.02762776 0.02487842 0.02551315
[13] 0.02551315 0.02618111 0.02551315 0.02702703 0.02631579 0.02631579
[19] 0.02564103 0.02564103 0.02564103 0.05263158 0.02857143 0.03225806
[25] 0.06250000 0.03571429 0.03846154 0.04000000 0.04166667 0.04761905
[31] 0.05000000 0.04545455 0.05000000 0.05555556 0.05882353 0.06666667
[37] 0.07142857 0.09090909 0.10000000 0.25000000 0.33333333 0.50000000

> ########################################################
> ### This time you get "-2LLR" = 1.150098 etc. etc.
> ##############################################################
> 
> 
> 
> cleanEx(); ..nameEx <- "emplikH1.test"
> 
> ### * emplikH1.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: emplikH1.test
> ### Title: Empirical likelihood for hazard with right censored, left
> ###   truncated data
> ### Aliases: emplikH1.test
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> fun <- function(x) { as.numeric(x <= 6.5) }
> emplikH1.test( x=c(1,2,3,4,5), d=c(1,1,0,1,1), theta=2, fun=fun) 
$"-2LLR"
[1] 0.006804456

$lambda
[1] -0.02672995

$times
[1] 1 2 4 5

$wts
[1] 0.2054928 0.2586419 0.5358051 1.0000000

$nits
NULL

$message
NULL

> fun2 <- function(x) {exp(-x)}  
> emplikH1.test( x=c(1,2,3,4,5), d=c(1,1,0,1,1), theta=0.2, fun=fun2) 
$"-2LLR"
[1] 0.5918754

$lambda
[1] -1.292062

$times
[1] 1 2 4 5

$wts
[1] 0.3811868 0.3199291 0.5314413 1.0000000

$nits
NULL

$message
NULL

> 
> 
> 
> cleanEx(); ..nameEx <- "emplikH2.test"
> 
> ### * emplikH2.test
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: emplikH2.test
> ### Title: Empirical likelihood for weighted hazard with right censored,
> ###   left truncated data
> ### Aliases: emplikH2.test
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> z1<-c(1,2,3,4,5)
> d1<-c(1,1,0,1,1)
> fun4 <- function(x, theta) { as.numeric(x <= theta) }
> emplikH2.test(x=z1,d=d1, K=0.5, fun=fun4, theta=3.5)
$"-2LLR"
[1] 0.02270185

$lambda
[1] -0.087689

$times
[1] 1 2 4 5

$wts
[1] 0.2192235 0.2807762 0.5000000 1.0000000

$nits
NULL

$message
NULL

> #Next, test if H(3.5) = log(2) .
> emplikH2.test(x=z1,d=d1, K=log(2), fun=fun4, theta=3.5)
$"-2LLR"
[1] 0.4263614

$lambda
[1] -0.3060978

$times
[1] 1 2 4 5

$wts
[1] 0.2882251 0.4049385 0.5000000 1.0000000

$nits
NULL

$message
NULL

> #Next, try one sample log rank test
> indi <- function(x,y){ as.numeric(x >= y) }
> fun3 <- function(t,z){rowsum(outer(z,t,FUN="indi"),group=rep(1,length(z)))} 
> emplikH2.test(x=z1, d=d1, K=sum(0.25* z1), fun=fun3, z=z1) 
$"-2LLR"
[1] 0.02206569

$lambda
[1] 0.01818070

$times
[1] 1 2 4 5

$wts
       [,1]      [,2]      [,3] [,4]
1 0.1833343 0.2291678 0.4583357    1

$nits
NULL

$message
NULL

> 
> 
> 
> cleanEx(); ..nameEx <- "emplikHs.disc2"
> 
> ### * emplikHs.disc2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: emplikHs.disc2
> ### Title: Two sample empirical likelihood ratio for discrete hazards with
> ###   right censored, left truncated data. Many constraints.
> ### Aliases: emplikHs.disc2
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> library(boot)
> data(channing)
> ymale <- channing[1:97,2]
> dmale <- channing[1:97,5]
> xmale <- channing[1:97,3]
> yfemale <- channing[98:462,2]
> dfemale <- channing[98:462,5]
> xfemale <- channing[98:462,3]
> fun1 <- function(x) { as.numeric(x <= 960) }
> ########################################################
> emplikHs.disc2(x1=xfemale, d1=dfemale, y1=yfemale, x2=xmale,d2=dmale, y2=ymale, theta=0.25, fun1=fun1, fun2=fun1)
         [,1]
[1,] 0.792554
$"-2LLR"
[1] 1.150098

$lambda
[1] -3.195481

$"-2LLR(sample1)"
[1] 0.06772553

$times1
  [1]  804  822  830  840  845  861  868  873  883  885  895  897  901  905  908
 [16]  911  912  915  919  923  926  928  930  931  932  934  936  940  941  944
 [31]  948  954  959  963  966  969  970  975  976  978  982  983  986  989  990
 [46]  991  992  994  995  996  998  999 1000 1001 1003 1004 1005 1006 1010 1011
 [61] 1012 1013 1014 1015 1017 1018 1019 1020 1021 1023 1024 1027 1029 1030 1033
 [76] 1040 1041 1043 1044 1047 1056 1063 1064 1068 1070 1072 1073 1074 1083 1084
 [91] 1085 1086 1089 1093 1097 1115 1122 1131 1132 1142 1152 1172 1192 1200

$times2
 [1]  777  869  872  876  893  894  898  907  909  911  927  932  945  948  957
[16]  966  969  971  983  985  989  993  998 1009 1012 1022 1025 1029 1031 1033
[31] 1036 1043 1044 1053 1055 1059 1060 1080 1085 1094 1128 1139

> ########################################################
> ### This time you get "-2LLR" = 1.150098 etc. etc.
> ##############################################################
> fun2 <- function(x){ cbind(as.numeric(x <= 960), as.numeric(x <= 860))}
> ############ fun2 has matrix output ###############
> emplikHs.disc2(x1=xfemale, d1=dfemale, y1=yfemale, x2=xmale,d2=dmale, y2=ymale, theta=c(0.25,0), fun1=fun2, fun2=fun2)
         [,1]     [,2]
[1,] 0.792554 0.561548
$"-2LLR"
[1] 1.554386

$lambda
[1]  1.119534 -6.007162

$"-2LLR(sample1)"
[1] 0.1342396

$times1
  [1]  804  822  830  840  845  861  868  873  883  885  895  897  901  905  908
 [16]  911  912  915  919  923  926  928  930  931  932  934  936  940  941  944
 [31]  948  954  959  963  966  969  970  975  976  978  982  983  986  989  990
 [46]  991  992  994  995  996  998  999 1000 1001 1003 1004 1005 1006 1010 1011
 [61] 1012 1013 1014 1015 1017 1018 1019 1020 1021 1023 1024 1027 1029 1030 1033
 [76] 1040 1041 1043 1044 1047 1056 1063 1064 1068 1070 1072 1073 1074 1083 1084
 [91] 1085 1086 1089 1093 1097 1115 1122 1131 1132 1142 1152 1172 1192 1200

$times2
 [1]  777  869  872  876  893  894  898  907  909  911  927  932  945  948  957
[16]  966  969  971  983  985  989  993  998 1009 1012 1022 1025 1029 1031 1033
[31] 1036 1043 1044 1053 1055 1059 1060 1080 1085 1094 1128 1139

> ################# you get "-2LLR" = 1.554386, etc ###########
> 
> 
> 
> cleanEx(); ..nameEx <- "emplikHs.test2"
> 
> ### * emplikHs.test2
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: emplikHs.test2
> ### Title: Two sample empirical likelihood ratio for hazards with right
> ###   censored, left truncated data
> ### Aliases: emplikHs.test2
> ### Keywords: nonparametric survival htest
> 
> ### ** Examples
> 
> library(boot)
> data(channing)
> ymale <- channing[1:97,2]
> dmale <- channing[1:97,5]
> xmale <- channing[1:97,3]
> yfemale <- channing[98:462,2]
> dfemale <- channing[98:462,5]
> xfemale <- channing[98:462,3]
> fun1 <- function(x) { as.numeric(x <= 960) }
> ########################################################
> fun2 <- function(x){ cbind(as.numeric(x <= 960), as.numeric(x <= 860))}
> ############ fun2 has matrix output ###############
> emplikHs.test2(x1=xfemale, d1=dfemale, y1=yfemale, x2=xmale,d2=dmale, y2=ymale, theta=c(0,0), fun1=fun2, fun2=fun2)
          [,1]      [,2]
[1,] -1.595082 -1.370471
$"-2LLR"
[1] 38.66794

$lambda
[1] -22.007684   2.121989

$"-2LLR(sample1)"
[1] 32.39357

> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  1.12 0.08 1.34 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')