# Example 5.2. Estimating a survival function with interval censored data.

#install.packages("Icens")

library(Icens)
data(cosmesis)
cosmesis.radonly<-subset(cosmesis, subset=(Trt==0), select=c(L,R))
cosmesis.radchemo<-subset(cosmesis, subset=(Trt==1), select=c(L,R))
NPMLE.radonly<-EM(cosmesis.radonly)
#in addition to EM, may also use (faster but less stable?)
# EMICM(), VEM(), ISDM(), PGM().

NPMLE.radchemo<-EM(cosmesis.radchemo)
#a warning is given on nonconvergence with the default maxiter=500;
# try more interations:
NPMLE.radchemo2<-EM(cosmesis.radchemo, maxiter=2000)
# it took NPMLE.radchemo2$numiter=714 to converge;
# but the two estimates are almost identical.

pdf("example5.2.pdf")
plot(NPMLE.radonly, surv=T, type="lc")
plot(NPMLE.radchemo, surv=T, type="lc", lty=2, new=F)
#plot(NPMLE.radchemo2, surv=T, type="lc", lty=4, new=F)
dev.off()

#see results with diff # of iterations:
> NPMLE.radonly
$pf
 [1] 4.634677e-02 3.336337e-02 8.866737e-02 7.075292e-02 0.000000e+00
 [6] 0.000000e+00 9.264584e-02 0.000000e+00 8.178576e-02 0.000000e+00
[11] 8.671167e-13 1.208798e-01 6.597546e-12 4.655581e-01

$numiter
[1] 370

$converge
[1] TRUE

$intmap
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,]    4    6    7   11   15   17   24   25   33    34    36    38    40    46
[2,]    5    7    8   12   16   18   25   26   34    35    37    40    44    48

attr(,"class")
[1] "icsurv"

NPMLE.radonly1<-EM(cosmesis.radonly, maxiter=1)
#Warning message:
#EM may have failed to converge in: EM(cosmesis.radonly, maxiter = 1)

### default tol=1e-12 is used to judge convergence!

NPMLE.radonly1$pf
 [1] 0.04066496 0.05570201 0.07259438 0.05244333 0.00986622 0.02454595
 [7] 0.05599433 0.02838403 0.05541771 0.03654528 0.05663381 0.09981587
[13] 0.10187317 0.30951895


NPMLE.radonly2<-EM(cosmesis.radonly, maxiter=2)
NPMLE.radonly2$pf
 [1] 0.041376793 0.052475833 0.076333016 0.058741849 0.004334793 0.020773815
 [7] 0.061246513 0.019210289 0.053469554 0.022239159 0.038752519 0.083777150
[13] 0.081972007 0.385296710