https://github.com/cran/emplik
Raw File
Tip revision: 2986b7285f1d1ab7b6e57d82fe3f20ec4616b3a0 authored by Mai Zhou on 08 August 1977, 00:00 UTC
version 0.8
Tip revision: 2986b72
emplik.Hs.test.R

emplik.Hs.test <- function(y, d, gmat, theta, error=1e-8, maxit=20) {


#Some simulation results to confirm that for right censored data and k (k>1)
#hazard constraints, the limiting distribution is chi-square with k df.
#
#> set.seed(123)
#> for(i in 1:1000) result[i] <- simu(200)
#> for(i in 1001:5000) result[i] <- simu(200)
#> simu
#function(num) {
#
#xx <- rexp(num)
#cc <- 2*rexp(num)
#yy <- pmin(xx, cc)
#dd <- as.numeric(cc > xx)
#
#temp0 <- summary(survfit(Surv(yy,dd),se.fit=FALSE,conf.type="none" ) )

temp0 <- DnR(y,d) 
times1 <- temp0$time
k <- length(times1)
d1 <- temp0$n.event
R1 <- temp0$n.risk
if(d1[k]==R1[k]) {
                   d1<-d1[-k]
                   R1 <- R1[-k]
                   times1 <- times1[-k]
                   }

#col1 <- 0.5*times1*as.numeric(times1 < 1 )
#col2 <- exp(-times1)
#col3 <- as.numeric(times1 < 0.9)
#theta <- c(-0.25,-1,-0.9)
#
#gmat <- cbind(col1,col2,col3)
#
#iqpHazards(times1,d1,R1,gmat,theta)$ChiSq
#
#}
#
#
#temp <- DnR()
#
#iqpHazards <- function(times1, d1, R1, gmat, theta, error=1e-9, maxit=20) {

### This function computes the emp. lik. ratio for one sample right 
### censored data---with many hazard constraints. The input should be 
### in the style of summary(survfit())$times; n.events; n.risk. 
### which can be obtained by DnR() and Wdataclean2() too.
### The gmat is the matrix of g() values on the times1. 
### theta is a vector that the constraints are equal to:
### sum g(times1) log( 1- dH(times))=int g(t) dlog(1-F(t)) = theta
### Last changed  7/2001 by  Mai Zhou mai@ms.uky.edu 
### Also, d1/R1 must be < 1. If d1(last)=R1(last) then delete this obs.
### When the constraints are true,  -2EmpLikRatio should behave as Chisq (k).

w01 <- d1/R1
u01 <- log(1-w01)
mat0a <- (exp(-u01/2)-exp(u01/2))/sqrt(d1)
mat0b <- R1 - d1 - (d1*exp(u01))/(1-exp(u01))
mat0c <- as.matrix(gmat)
dd <- u01 %*% mat0c
k <- length(dd)
mat0d <- theta - dd

# value0 <- solve.QP(diag(mat0a),mat0b,mat0c,mat0d,meq=k,factorized=TRUE)
value0 <- solve3.QP(mat0a,mat0b,mat0c,mat0d,meq=k,factorized=TRUE)
wp <- u01 + value0$solution

###### iteration

m<-1
matd <- rep(0, k)
diff <- 10

while(diff > error && m < maxit) {
    mata <- (exp(-wp/2)-exp(wp/2))/sqrt(d1)
    matb <- R1-d1-(d1*exp(wp))/(1-exp(wp))
##  value <- solve.QP(diag(mata), matb, mat0c, matd, meq=k, factorized=TRUE)
 value <- solve3.QP(mata, matb, mat0c, matd, meq=k, factorized=TRUE)
    wp <- wp+value$solution
    diff <- sum(abs(value$solution))
    m<-m+1
    }
loglik2 <- sum( d1*log(1-exp(wp))+(R1-d1)*wp)
loglik1 <- sum(d1*log(w01)+(R1-d1)*u01)
list(weight = 1-exp(wp), ChiSq=2*(loglik1-loglik2))
}

#The Q-Q plot is in pp2.ps It agrees with chi square 3 df quite well.
back to top