https://github.com/cran/emplik
Tip revision: 31d32093eee1aa45b156c2a215f600da3894af7f authored by Mai Zhou on 26 December 2006, 00:00:00 UTC
version 0.9-3-1
version 0.9-3-1
Tip revision: 31d3209
emplikH.disc2.R
#####################################################
#### 2 sample, discrete hazard type constraint. #####
#### Only handles one constraint. #####
#####################################################
emplikH.disc2 <- function(x1, d1, y1 = -Inf, x2, d2, y2 = -Inf,
theta, fun1, fun2, tola = 1e-6, maxi, mini) {
########Sample One########
n1 <- length(x1)
if (n1 <= 2)
stop("Need more observations")
if (length(d1) != n1)
stop("length of x1 and d1 must agree")
if (any((d1 != 0) & (d1 != 1)))
stop("d1 must be 0/1's for censor/not-censor")
if (!is.numeric(x1))
stop("x1 must be numeric values --- observed times")
newdata1 <- Wdataclean2(z=x1, d=d1)
temp1 <- DnR(newdata1$value, newdata1$dd, newdata1$weight, y=y1)
jump1 <- (temp1$n.event)/temp1$n.risk
k1 <- temp1$n.event - temp1$n.risk
index1 <- (jump1 < 1)
k1 <- k1[index1]
eve1 <- temp1$n.event[index1]
tm1 <- temp1$times[index1]
rsk1 <- temp1$n.risk[index1]
jmp1 <- jump1[index1]
funtime1 <- fun1(tm1)
funh1 <- funtime1/rsk1
########Sample two########
n2 <- length(x2)
if (n2 <= 2)
stop("Need more observations for sample 2")
if (length(d2) != n2)
stop("length of x2 and d2 must agree")
if (any((d2 != 0) & (d2 != 1)))
stop("d2 must be 0/1's for censor/not-censor")
if (!is.numeric(x2))
stop("x2 must be numeric values -- observed times")
newdata2 <- Wdataclean2(z=x2, d=d2)
temp2 <- DnR(newdata2$value, newdata2$dd, newdata2$weight, y=y2)
jump2 <- (temp2$n.event)/temp2$n.risk
k2 <- temp2$n.event - temp2$n.risk
index2 <- (jump2 < 1)
k2 <- k2[index2]
eve2 <- temp2$n.event[index2]
tm2 <- temp2$times[index2]
rsk2 <- temp2$n.risk[index2]
jmp2 <- jump2[index2]
funtime2 <- fun2(tm2)
funh2 <- funtime2/rsk2
######## constrain function ########
inthaz <- function(x, funt1, evt1, rsk1, funt2, evt2, rsk2, tht) {
sum(funt1*log(1-(evt1/(rsk1+x*funt1)))) -
sum(funt2*log(1-(evt2/(rsk2-x*funt2)))) - tht }
m0 <- inthaz(0, funtime1, eve1, rsk1, funtime2, eve2, rsk2, theta)
if (m0 >0) maxi <- 0
if (m0 <0) mini <- 0
m1 <- inthaz(mini+tola, funtime1, eve1, rsk1, funtime2, eve2, rsk2, theta)
m2 <- inthaz(maxi-tola, funtime1, eve1, rsk1, funtime2, eve2, rsk2, theta)
print(c(m0, m1, m2))
# print(m1)
# print(m2)
temp <- uniroot(inthaz, c(mini, maxi), tol=tola, funt1=funtime1,
evt1=eve1, rsk1=rsk1, funt2=funtime2, evt2=eve2, rsk2=rsk2, tht=theta)
lam <- temp$root
onePlam <- 1+lam*funh1
weights1 <- jmp1/onePlam
oneMlam <- 1-lam*funh2
weights2 <- jmp2/oneMlam
loglik1 <- sum(eve1*log(onePlam))+sum((-k1)*log((1-jmp1)/(1-weights1)))
loglik2 <- sum(eve2*log(oneMlam))+sum((-k2)*log((1-jmp2)/(1-weights2)))
loglikR <- 2*(loglik1+loglik2)
#MZ <- inthaz(lam, funtime1, eve1, rsk1, funtime2, eve2, rsk2, theta)
#print(MZ)
list("-2LLR" = loglikR, lambda = lam, times1 = tm1, times2=tm2,
wts1 = weights1, wts2 = weights2)
}