yfyang86 / kmc

Kaplan Meier estimator with Constraints
2 stars 3 forks source link

New rootSolve causes failure #5

Closed yfyang86 closed 2 years ago

yfyang86 commented 3 years ago

Verison: 0.2+ BUG: root solving procedure causes -2LLR to be negative.

The potential issue is that rootSolve package changes in the past.

emplik results:

```R library(KMsurv) library(emplik) data(kidney) RMSTfun <- function(x) {pmin(x,20) - 15} el.cen.EM2( x = kidney$time, d = kidney$delta, fun = RMSTfun, mu = 0, maxit = 50) ``` reports ```R $loglik [1] -131.105 $times [1] 0.5 0.5 0.5 0.5 0.5 0.5 1.5 2.5 2.5 3.5 3.5 4.5 4.5 5.5 6.5 [16] 8.5 8.5 9.5 10.5 11.5 15.5 15.5 16.5 18.5 23.5 26.5 28.5 $prob [1] 0.01122308 0.01122308 0.01122308 0.01122308 0.01122308 0.01122308 [7] 0.01243132 0.01279473 0.01279473 0.01381316 0.01381316 0.01533046 [13] 0.01533046 0.01599773 0.01746515 0.02043914 0.02043914 0.02186346 [19] 0.02285866 0.02485554 0.03538362 0.03538362 0.03297763 0.03220609 [25] 0.06200016 0.10114908 0.39333447 $lam [1] 2.061926 $iters [1] 50 $`-2LLR` [1] 2.51874 $Pval [1] 0.1125004 ```

kmc results

library(kmc)
kmc.solve(kidney$time, kidney$delta, list(RMSTfun))
myfun9 <- function(theta, x, d){
    ff <- function(t) pmin(t, 20) - theta
kmc.solve(x, d, g = list(ff)) 
}
findUL(step = 1, fun = myfun9, MLE = 16, x = kidney$time, d = kidney$delta)

$Low
[1] 12.97085

$Up
[1] 19.56962

$FstepL
[1] 6.103516e-05

$FstepU
[1] 8.69197e-05

$Lvalue
[1] 3.840001

$Uvalue
[1] -198.9877

Warning messages:
1: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!

2: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!

3: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!

4: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!

5: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!

6: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!

7: In kmc.solve(x, d, g = list(ff)) :
The results may be not feasible!
yfyang86 commented 3 years ago

After revision, in this case, lambda should meet the following requirements:

  1. sum(prob) = 1

therefore, for the rootSolve(.) function. we need to check the root that:

  1. in the branch that contains 0
  2. satisfies sum(RMSTfun(time) * prob) == 1

An older version checks this, but the newer one fails.

The next version will fix this.

library(KMsurv)
library(emplik)
library(kmc)
library(rootSolve)
data(kidney)
kmc.clean <- function(kmc.time,delta){
  n <- length(kmc.time)
  dataOrder <- order(kmc.time, -delta)
  kmc.time <- kmc.time[dataOrder]
  delta <- delta[dataOrder] 
  FirstUnCenLocation<-which(delta==1)[1];
  if (FirstUnCenLocation==n) {stop('Only one uncensored point.');}
  if (FirstUnCenLocation!=1){
    delta=delta[FirstUnCenLocation:n];
    kmc.time=kmc.time[FirstUnCenLocation:n];
  }
  delta[length(kmc.time)]=1;
  return (list(kmc.time=kmc.time,delta=delta));
}

RMSTfun <- function(x) {pmin(x,20) - 15}

x = kidney$time
d = kidney$delta
inputData = kmc.clean(x, d)
Gmat = RMSTfun(inputData$kmc.time)

##  `kmc_native` is the KMC Solver (1-dim) which returns `prob`
fun.C <- Vectorize(
  function(lam){
    return(sum(Gmat*kmc_native(inputData$delta, lam*Gmat)));
  }
)

All = uniroot.all(fun.C, c(-3, 2), tol = 1e-8)
s = el.cen.EM2( x  = kidney$time, d = kidney$delta, fun = RMSTfun, mu = 0, maxit = 50)

curve(fun.C(x), -3, 5,main = "uniroot.all", ylim=c(-10,10))
points(All, y = rep(0, length(All)), pch = 16, cex = 2, col='steelblue')
points(-s$lam, y = 0, pch = 16, cex = 2, col='red')
yfyang86 commented 2 years ago

@yfyang86 26abbb3

Fixed