R/bcontSurvGunivMIXED_ExcessHazard_LeftTruncation.R
In KironmoyDas/KD-STAT0035-GMupdate: Generalised Joint Regression Modelling
Defines functions
bcontSurvGunivMIXED_ExcessHazard_LeftTruncation
Documented in
bcontSurvGunivMIXED_ExcessHazard_LeftTruncation
bcontSurvGunivMIXED_ExcessHazard_LeftTruncation <- function(params, respvec, VC, ps, AT = FALSE){
c1 <- VC$hrate
c2 <- exp(-VC$d.lchrate)
c3 <- exp(-VC$d.rchrate)
c2.td <- exp(-VC$d.lchrate.td)
c3.td <- exp(-VC$d.rchrate.td)
# C.i = VC$indvU*(-VC$d.lchrate) + VC$indvR*(-VC$d.lchrate)
# C.i.td = VC$indvUT*(-VC$d.chrate.td) + VC$indvRT*(-VC$d.lchrate.td)
#################################################################################
monP <- monP1 <- monP2 <- k <- 0; V <- list()
etad <- etas1 <- etas2 <- l.ln <- NULL
params1 <- params[1:VC$X1.d2]
params1[VC$mono.sm.pos] <- exp( params1[VC$mono.sm.pos] )
eta1 <- VC$X1%*%params1
Xd1P <- VC$Xd1%*%params1
indN <- as.numeric(Xd1P < 0)
#if(!is.null(VC$indexT)) print(table(indN))
Xd1P <- ifelse(Xd1P < VC$min.dn, VC$min.dn, Xd1P )
if( any(indN == TRUE) && !is.null(VC$indexT) ){
monP22 <- matrix(0, length(params),length(params))
for(i in 1:length(VC$pos.pb)){
V[[i]] <- as.numeric(diff(params1[ VC$pos.pb[[i]] ]) < 0)
monP22[ VC$pos.pb[[i]], VC$pos.pb[[i]] ] <- t(VC$D[[i]]*V[[i]])%*%VC$D[[i]]
}
k <- VC$my.env$k
monP2 <- k*monP22
monP <- k/2*crossprod(params, monP22)%*%params
monP1 <- k*(monP22%*%params)
VC$my.env$k <- k*2
}
##################
der.par1 <- der2.par1 <- params1
der.par1[-c( VC$mono.sm.pos )] <- 1
der2.par1[-c( VC$mono.sm.pos )] <- 0
##################
pd1 <- probmS(eta1, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p1 <- pd1$pr
dS1eta1 <- pd1$dS
d2S1eta1 <- pd1$d2S
d3S1eta1 <- pd1$d3S
der2eta1dery1b1 <- t(t(VC$Xd1)*der.par1)
dereta1derb1 <- t(t(VC$X1)*der.par1)
if( any(unique(VC$cens) == "I") ){
eta2 <- VC$X2%*%params1
pd2 <- probmS(eta2, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p2 <- pd2$pr
dS2eta2 <- pd2$dS
d2S2eta2 <- pd2$d2S
d3S2eta2 <- pd2$d3S
dereta2derb1 <- t(t(VC$X2)*der.par1)
}else{
eta2 <- rep(0.1, length(eta1)) # just to get the code running, a bit inefficient but clear what is below
pd2 <- probmS(eta2, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p2 <- pd2$pr
dS2eta2 <- pd2$dS
d2S2eta2 <- pd2$d2S
d3S2eta2 <- pd2$d3S
dereta2derb1 <- t(t(VC$X2)*der.par1)
}
# THIRD PREDICTOR due to the LEFT-TRUNCATED observations
if( any(unique(VC$cens) %in% c('UT', 'LT', 'RT', 'IT') ) ){
eta3 <- VC$X3%*%params1
pd3 <- probmS(eta3, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p3 <- pd3$pr
dS3eta3 <- pd3$dS
d2S3eta3 <- pd3$d2S
d3S3eta3 <- pd3$d3S
dereta3derb1 <- t(t(VC$X3)*der.par1)
} else {
eta3 <- rep(0.1, length(eta1)) # just to get the code running
pd3 <- probmS(eta3, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p3 <- pd3$pr
dS3eta3 <- pd3$dS
d2S3eta3 <- pd3$d2S
d3S3eta3 <- pd3$d3S
dereta3derb1 <- t(t(VC$X3)*der.par1)
}
####################################################################################################################
l.par <- VC$weights*(
VC$indvU*log( c1*p1 - dS1eta1*Xd1P ) +
VC$indvR*log( p1 ) +
VC$indvL*log( 1 - p1*c2 ) +
VC$indvI*log( mm(p1*c2 - p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr) ) +
VC$indvUT*log( c1*p3^-1*p1 - p3^-1*dS1eta1*Xd1P ) +
VC$indvRT*log( p3^-1*p1 ) +
VC$indvLT*log( mm( 1 - p3^-1*p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr) ) +
VC$indvIT*log( mm(p3^-1*p1*c2.td - p3^-1*p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr) ) # +
# C.i + C.i.td
)
res <- -sum(l.par)
########################################################################################################################
dl.dbe1 <- -VC$weights*(
VC$indvU*( c((c1*p1 - dS1eta1*Xd1P)^-1)*( c(c1*dS1eta1)*dereta1derb1 - c(d2S1eta1*Xd1P)*dereta1derb1 - c(dS1eta1)*der2eta1dery1b1 ) ) +
VC$indvR*c(p1^-1*dS1eta1)*dereta1derb1 +
VC$indvL*-c( ( 1 - p1*c2 )^-1*dS1eta1*c2 )*dereta1derb1 + # minus here in front from L censoring
VC$indvI*c( mm( p1*c2 - p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr ) )^-1*(c(dS1eta1*c2)*dereta1derb1 - c(dS2eta2*c3)*dereta2derb1) +
VC$indvUT*( c((c1*p1 - dS1eta1*Xd1P)^-1)*( c(c1*dS1eta1)*dereta1derb1 - c(d2S1eta1*Xd1P)*dereta1derb1 - c(dS1eta1)*der2eta1dery1b1 ) - c(p3^-1*dS3eta3)*dereta3derb1 ) +
VC$indvRT*( c(p1^-1*dS1eta1)*dereta1derb1 - c(p3^-1*dS3eta3)*dereta3derb1 ) +
VC$indvLT*( -c( mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1*dS1eta1*c2.td )*dereta1derb1 + c( mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1*dS3eta3 - p3^-1*dS3eta3 )*dereta3derb1 ) +
VC$indvIT*( c( mm( p1*c2.td - p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr ) )^-1*( c(dS1eta1*c2.td)*dereta1derb1 - c(dS2eta2*c3.td)*dereta2derb1 ) - c(p3^-1*dS3eta3)*dereta3derb1 )
)
G <- colSums(dl.dbe1)
########################################################################################################################################################################
H <- -(
# UNCENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvU*( ( c1*p1 - dS1eta1*Xd1P )^-1*(c1*d2S1eta1) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-1*(d3S1eta1*Xd1P) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*dereta1derb1, der2eta1dery1b1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*der2eta1dery1b1, dereta1derb1 ) +
diag( colSums( t( t(c(VC$weights*VC$indvU*( c1*p1 - dS1eta1*Xd1P )^-1*(c1*dS1eta1) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvU*-( c1*p1 - dS1eta1*Xd1P )^-1*(d2S1eta1*Xd1P) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvU*-( c1*p1 - dS1eta1*Xd1P )^-1*dS1eta1 )*VC$Xd1)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*(d2S1eta1*Xd1P)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*dS1eta1^2 ))*der2eta1dery1b1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvU*(2*( c1*p1 - dS1eta1*Xd1P )^-2*( (c1*dS1eta1)*(d2S1eta1*Xd1P) ) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvU*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvU*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*der2eta1dery1b1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*der2eta1dery1b1, dereta1derb1) +
# RIGHT CENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvR*(-p1^-2*dS1eta1^2+p1^-1*d2S1eta1))*dereta1derb1, dereta1derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvR*p1^-1*dS1eta1)*VC$X1)*der2.par1 ) ) ) -
# LEFT CENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvL*(( 1 - p1*c2 )^-2*(dS1eta1*c2)^2 +
( 1 - p1*c2 )^-1*(d2S1eta1*c2) ))*dereta1derb1, dereta1derb1) -
diag( colSums( t( t(c(VC$weights*VC$indvL*( 1 - p1*c2 )^-1*(dS1eta1*c2) )*VC$X1)*der2.par1 ) ) ) +
# INTERVAL CENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvI*(mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S1eta1*c2) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvI*(-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S2eta2*c3) ))*dereta2derb1, dereta2derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvI*mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS1eta1*c2) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvI*-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS2eta2*c3) )*VC$X2)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvI*(-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2)^2))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvI*(-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS2eta2*c3)^2))*dereta2derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvI*(mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2)*(dS2eta2*c3) ))*dereta1derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvI*(mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2)*(dS2eta2*c3) ))*dereta2derb1, dereta1derb1) +
# UNCENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvUT*( ( c1*p1 - dS1eta1*Xd1P )^-1*(c1*d2S1eta1) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-1*(d3S1eta1*Xd1P) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*dereta1derb1, der2eta1dery1b1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*der2eta1dery1b1, dereta1derb1 ) +
diag( colSums( t( t(c(VC$weights*VC$indvUT*( c1*p1 - dS1eta1*Xd1P )^-1*(c1*dS1eta1) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvUT*-( c1*p1 - dS1eta1*Xd1P )^-1*(d2S1eta1*Xd1P) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvUT*-( c1*p1 - dS1eta1*Xd1P )^-1*dS1eta1 )*VC$Xd1)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*(d2S1eta1*Xd1P)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*dS1eta1^2 ))*der2eta1dery1b1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvUT*(2*( c1*p1 - dS1eta1*Xd1P )^-2*( (c1*dS1eta1)*(d2S1eta1*Xd1P) ) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvUT*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*der2eta1dery1b1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*der2eta1dery1b1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*(-p3^-1*d2S3eta3 ))*dereta3derb1, dereta3derb1 ) +
diag( colSums( t( t(-c(VC$weights*VC$indvUT*p3^-1*dS3eta3 )*VC$X3)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvUT*( p3^-2*dS3eta3^2 ))*dereta3derb1, dereta3derb1 ) +
# RIGHT CENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvRT*(-p1^-2*dS1eta1^2 + p1^-1*d2S1eta1))*dereta1derb1, dereta1derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvRT*p1^-1*dS1eta1)*VC$X1)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvRT*(p3^-2*dS3eta3^2 - p3^-1*d2S3eta3))*dereta3derb1, dereta3derb1) +
diag( colSums( t( t(-c(VC$weights*VC$indvRT*p3^-1*dS3eta3)*VC$X3)*der2.par1 ) ) ) +
# LEFT CENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvLT*( -mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)^2 - mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*d2S1eta1*c2.td ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvLT*( (p3^-2 - mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2)*dS3eta3^2 - (p3^-1 - mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1)*d2S3eta3 ))*dereta3derb1, dereta3derb1) +
crossprod(c(VC$weights*VC$indvLT*( mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*dS3eta3*dS1eta1*c2.td ))*dereta1derb1, dereta3derb1) +
crossprod(c(VC$weights*VC$indvLT*( mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*dS1eta1*dS3eta3*c2.td ))*dereta3derb1, dereta1derb1) +
diag( colSums( t( t(-c(VC$weights*VC$indvLT*mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1*(dS1eta1*c2.td) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(-c(VC$weights*VC$indvLT*( p3^-1 - mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1 )*dS3eta3 )*VC$X3)*der2.par1 ) ) ) +
# INTERVAL CENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvIT*(mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S1eta1*c2.td) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvIT*(-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S2eta2*c3.td) ))*dereta2derb1, dereta2derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvIT*mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS1eta1*c2.td) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvIT*-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS2eta2*c3.td) )*VC$X2)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvIT*(-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)^2))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvIT*(-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS2eta2*c3.td)^2))*dereta2derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvIT*(mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)*(dS2eta2*c3.td) ))*dereta1derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvIT*(mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)*(dS2eta2*c3.td) ))*dereta2derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvIT*( p3^-2*dS3eta3^2 - p3^-1*d2S3eta3 ))*dereta3derb1, dereta3derb1 ) +
diag( colSums( t( t(-c(VC$weights*VC$indvIT*p3^-1*dS3eta3 )*VC$X3)*der2.par1 ) ) )
)
############################################################################################################################################################################################
if(VC$extra.regI == "pC") H <- regH(H, type = 1)
S.h <- ps$S.h + monP2 # hess
S.h1 <- 0.5*crossprod(params, ps$S.h)%*%params + monP # lik
S.h2 <- S.h%*%params + monP1 # grad
S.res <- res
res <- S.res + S.h1
G <- G + S.h2
H <- H + S.h
if(VC$extra.regI == "sED") H <- regH(H, type = 2)
list(value = res, gradient = G, hessian = H, S.h = S.h, S.h1 = S.h1, S.h2 = S.h2, indN = indN, V = V,
l = S.res, l.ln = l.ln, l.par = l.par, ps = ps, k = VC$my.env$k, monP2 = monP2, params1 = params1,
eta1 = eta1, p1 = p1, p2 = p2, p3 = p3,
pdf1 = -dS1eta1, pdf2 = -dS2eta2, pdf3 = -dS3eta3,
dl.dbe1 = NULL, dl.dbe2 = NULL, dl.dteta.st = NULL,
hrate = VC$hrate,
d.lchrate = VC$d.lchrate,
d.rchrate = VC$d.rchrate,
d.lchrate.td = VC$d.lchrate.td,
d.rchrate.td = VC$d.rchrate.td)
}
KironmoyDas/KD-STAT0035-GMupdate documentation built on Feb. 15, 2021, 12:17 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions bcontSurvGunivMIXED_ExcessHazard_LeftTruncation
Documented in bcontSurvGunivMIXED_ExcessHazard_LeftTruncation
bcontSurvGunivMIXED_ExcessHazard_LeftTruncation <- function(params, respvec, VC, ps, AT = FALSE){
c1 <- VC$hrate
c2 <- exp(-VC$d.lchrate)
c3 <- exp(-VC$d.rchrate)
c2.td <- exp(-VC$d.lchrate.td)
c3.td <- exp(-VC$d.rchrate.td)
# C.i = VC$indvU*(-VC$d.lchrate) + VC$indvR*(-VC$d.lchrate)
# C.i.td = VC$indvUT*(-VC$d.chrate.td) + VC$indvRT*(-VC$d.lchrate.td)
#################################################################################
monP <- monP1 <- monP2 <- k <- 0; V <- list()
etad <- etas1 <- etas2 <- l.ln <- NULL
params1 <- params[1:VC$X1.d2]
params1[VC$mono.sm.pos] <- exp( params1[VC$mono.sm.pos] )
eta1 <- VC$X1%*%params1
Xd1P <- VC$Xd1%*%params1
indN <- as.numeric(Xd1P < 0)
#if(!is.null(VC$indexT)) print(table(indN))
Xd1P <- ifelse(Xd1P < VC$min.dn, VC$min.dn, Xd1P )
if( any(indN == TRUE) && !is.null(VC$indexT) ){
monP22 <- matrix(0, length(params),length(params))
for(i in 1:length(VC$pos.pb)){
V[[i]] <- as.numeric(diff(params1[ VC$pos.pb[[i]] ]) < 0)
monP22[ VC$pos.pb[[i]], VC$pos.pb[[i]] ] <- t(VC$D[[i]]*V[[i]])%*%VC$D[[i]]
}
k <- VC$my.env$k
monP2 <- k*monP22
monP <- k/2*crossprod(params, monP22)%*%params
monP1 <- k*(monP22%*%params)
VC$my.env$k <- k*2
}
##################
der.par1 <- der2.par1 <- params1
der.par1[-c( VC$mono.sm.pos )] <- 1
der2.par1[-c( VC$mono.sm.pos )] <- 0
##################
pd1 <- probmS(eta1, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p1 <- pd1$pr
dS1eta1 <- pd1$dS
d2S1eta1 <- pd1$d2S
d3S1eta1 <- pd1$d3S
der2eta1dery1b1 <- t(t(VC$Xd1)*der.par1)
dereta1derb1 <- t(t(VC$X1)*der.par1)
if( any(unique(VC$cens) == "I") ){
eta2 <- VC$X2%*%params1
pd2 <- probmS(eta2, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p2 <- pd2$pr
dS2eta2 <- pd2$dS
d2S2eta2 <- pd2$d2S
d3S2eta2 <- pd2$d3S
dereta2derb1 <- t(t(VC$X2)*der.par1)
}else{
eta2 <- rep(0.1, length(eta1)) # just to get the code running, a bit inefficient but clear what is below
pd2 <- probmS(eta2, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p2 <- pd2$pr
dS2eta2 <- pd2$dS
d2S2eta2 <- pd2$d2S
d3S2eta2 <- pd2$d3S
dereta2derb1 <- t(t(VC$X2)*der.par1)
}
# THIRD PREDICTOR due to the LEFT-TRUNCATED observations
if( any(unique(VC$cens) %in% c('UT', 'LT', 'RT', 'IT') ) ){
eta3 <- VC$X3%*%params1
pd3 <- probmS(eta3, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p3 <- pd3$pr
dS3eta3 <- pd3$dS
d2S3eta3 <- pd3$d2S
d3S3eta3 <- pd3$d3S
dereta3derb1 <- t(t(VC$X3)*der.par1)
} else {
eta3 <- rep(0.1, length(eta1)) # just to get the code running
pd3 <- probmS(eta3, VC$margins[1], min.dn = VC$min.dn, min.pr = VC$min.pr, max.pr = VC$max.pr)
p3 <- pd3$pr
dS3eta3 <- pd3$dS
d2S3eta3 <- pd3$d2S
d3S3eta3 <- pd3$d3S
dereta3derb1 <- t(t(VC$X3)*der.par1)
}
####################################################################################################################
l.par <- VC$weights*(
VC$indvU*log( c1*p1 - dS1eta1*Xd1P ) +
VC$indvR*log( p1 ) +
VC$indvL*log( 1 - p1*c2 ) +
VC$indvI*log( mm(p1*c2 - p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr) ) +
VC$indvUT*log( c1*p3^-1*p1 - p3^-1*dS1eta1*Xd1P ) +
VC$indvRT*log( p3^-1*p1 ) +
VC$indvLT*log( mm( 1 - p3^-1*p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr) ) +
VC$indvIT*log( mm(p3^-1*p1*c2.td - p3^-1*p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr) ) # +
# C.i + C.i.td
)
res <- -sum(l.par)
########################################################################################################################
dl.dbe1 <- -VC$weights*(
VC$indvU*( c((c1*p1 - dS1eta1*Xd1P)^-1)*( c(c1*dS1eta1)*dereta1derb1 - c(d2S1eta1*Xd1P)*dereta1derb1 - c(dS1eta1)*der2eta1dery1b1 ) ) +
VC$indvR*c(p1^-1*dS1eta1)*dereta1derb1 +
VC$indvL*-c( ( 1 - p1*c2 )^-1*dS1eta1*c2 )*dereta1derb1 + # minus here in front from L censoring
VC$indvI*c( mm( p1*c2 - p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr ) )^-1*(c(dS1eta1*c2)*dereta1derb1 - c(dS2eta2*c3)*dereta2derb1) +
VC$indvUT*( c((c1*p1 - dS1eta1*Xd1P)^-1)*( c(c1*dS1eta1)*dereta1derb1 - c(d2S1eta1*Xd1P)*dereta1derb1 - c(dS1eta1)*der2eta1dery1b1 ) - c(p3^-1*dS3eta3)*dereta3derb1 ) +
VC$indvRT*( c(p1^-1*dS1eta1)*dereta1derb1 - c(p3^-1*dS3eta3)*dereta3derb1 ) +
VC$indvLT*( -c( mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1*dS1eta1*c2.td )*dereta1derb1 + c( mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1*dS3eta3 - p3^-1*dS3eta3 )*dereta3derb1 ) +
VC$indvIT*( c( mm( p1*c2.td - p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr ) )^-1*( c(dS1eta1*c2.td)*dereta1derb1 - c(dS2eta2*c3.td)*dereta2derb1 ) - c(p3^-1*dS3eta3)*dereta3derb1 )
)
G <- colSums(dl.dbe1)
########################################################################################################################################################################
H <- -(
# UNCENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvU*( ( c1*p1 - dS1eta1*Xd1P )^-1*(c1*d2S1eta1) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-1*(d3S1eta1*Xd1P) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*dereta1derb1, der2eta1dery1b1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*der2eta1dery1b1, dereta1derb1 ) +
diag( colSums( t( t(c(VC$weights*VC$indvU*( c1*p1 - dS1eta1*Xd1P )^-1*(c1*dS1eta1) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvU*-( c1*p1 - dS1eta1*Xd1P )^-1*(d2S1eta1*Xd1P) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvU*-( c1*p1 - dS1eta1*Xd1P )^-1*dS1eta1 )*VC$Xd1)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*(d2S1eta1*Xd1P)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*dS1eta1^2 ))*der2eta1dery1b1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvU*(2*( c1*p1 - dS1eta1*Xd1P )^-2*( (c1*dS1eta1)*(d2S1eta1*Xd1P) ) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvU*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvU*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*der2eta1dery1b1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvU*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*der2eta1dery1b1, dereta1derb1) +
# RIGHT CENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvR*(-p1^-2*dS1eta1^2+p1^-1*d2S1eta1))*dereta1derb1, dereta1derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvR*p1^-1*dS1eta1)*VC$X1)*der2.par1 ) ) ) -
# LEFT CENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvL*(( 1 - p1*c2 )^-2*(dS1eta1*c2)^2 +
( 1 - p1*c2 )^-1*(d2S1eta1*c2) ))*dereta1derb1, dereta1derb1) -
diag( colSums( t( t(c(VC$weights*VC$indvL*( 1 - p1*c2 )^-1*(dS1eta1*c2) )*VC$X1)*der2.par1 ) ) ) +
# INTERVAL CENSORED CONTRIBUTION
crossprod(c(VC$weights*VC$indvI*(mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S1eta1*c2) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvI*(-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S2eta2*c3) ))*dereta2derb1, dereta2derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvI*mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS1eta1*c2) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvI*-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS2eta2*c3) )*VC$X2)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvI*(-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2)^2))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvI*(-mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS2eta2*c3)^2))*dereta2derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvI*(mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2)*(dS2eta2*c3) ))*dereta1derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvI*(mm(p1*c2-p2*c3, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2)*(dS2eta2*c3) ))*dereta2derb1, dereta1derb1) +
# UNCENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvUT*( ( c1*p1 - dS1eta1*Xd1P )^-1*(c1*d2S1eta1) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-1*(d3S1eta1*Xd1P) ))*dereta1derb1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*dereta1derb1, der2eta1dery1b1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-1*d2S1eta1 ))*der2eta1dery1b1, dereta1derb1 ) +
diag( colSums( t( t(c(VC$weights*VC$indvUT*( c1*p1 - dS1eta1*Xd1P )^-1*(c1*dS1eta1) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvUT*-( c1*p1 - dS1eta1*Xd1P )^-1*(d2S1eta1*Xd1P) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvUT*-( c1*p1 - dS1eta1*Xd1P )^-1*dS1eta1 )*VC$Xd1)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*(d2S1eta1*Xd1P)^2 ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*dS1eta1^2 ))*der2eta1dery1b1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvUT*(2*( c1*p1 - dS1eta1*Xd1P )^-2*( (c1*dS1eta1)*(d2S1eta1*Xd1P) ) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvUT*( ( c1*p1 - dS1eta1*Xd1P )^-2*(c1*dS1eta1^2) ))*der2eta1dery1b1, dereta1derb1 ) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*dereta1derb1, der2eta1dery1b1) +
crossprod(c(VC$weights*VC$indvUT*(-( c1*p1 - dS1eta1*Xd1P )^-2*( (d2S1eta1*Xd1P)*dS1eta1 ) ))*der2eta1dery1b1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvUT*(-p3^-1*d2S3eta3 ))*dereta3derb1, dereta3derb1 ) +
diag( colSums( t( t(-c(VC$weights*VC$indvUT*p3^-1*dS3eta3 )*VC$X3)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvUT*( p3^-2*dS3eta3^2 ))*dereta3derb1, dereta3derb1 ) +
# RIGHT CENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvRT*(-p1^-2*dS1eta1^2 + p1^-1*d2S1eta1))*dereta1derb1, dereta1derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvRT*p1^-1*dS1eta1)*VC$X1)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvRT*(p3^-2*dS3eta3^2 - p3^-1*d2S3eta3))*dereta3derb1, dereta3derb1) +
diag( colSums( t( t(-c(VC$weights*VC$indvRT*p3^-1*dS3eta3)*VC$X3)*der2.par1 ) ) ) +
# LEFT CENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvLT*( -mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)^2 - mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*d2S1eta1*c2.td ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvLT*( (p3^-2 - mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2)*dS3eta3^2 - (p3^-1 - mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1)*d2S3eta3 ))*dereta3derb1, dereta3derb1) +
crossprod(c(VC$weights*VC$indvLT*( mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*dS3eta3*dS1eta1*c2.td ))*dereta1derb1, dereta3derb1) +
crossprod(c(VC$weights*VC$indvLT*( mm(p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*dS1eta1*dS3eta3*c2.td ))*dereta3derb1, dereta1derb1) +
diag( colSums( t( t(-c(VC$weights*VC$indvLT*mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1*(dS1eta1*c2.td) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(-c(VC$weights*VC$indvLT*( p3^-1 - mm( p3 - p1*c2.td, min.pr = VC$min.pr, max.pr = VC$max.pr )^-1 )*dS3eta3 )*VC$X3)*der2.par1 ) ) ) +
# INTERVAL CENSORED LEFT-TRUNCATED CONTRIBUTION
crossprod(c(VC$weights*VC$indvIT*(mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S1eta1*c2.td) ))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvIT*(-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(d2S2eta2*c3.td) ))*dereta2derb1, dereta2derb1) +
diag( colSums( t( t(c(VC$weights*VC$indvIT*mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS1eta1*c2.td) )*VC$X1)*der2.par1 ) ) ) +
diag( colSums( t( t(c(VC$weights*VC$indvIT*-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-1*(dS2eta2*c3.td) )*VC$X2)*der2.par1 ) ) ) +
crossprod(c(VC$weights*VC$indvIT*(-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)^2))*dereta1derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvIT*(-mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS2eta2*c3.td)^2))*dereta2derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvIT*(mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)*(dS2eta2*c3.td) ))*dereta1derb1, dereta2derb1) +
crossprod(c(VC$weights*VC$indvIT*(mm(p1*c2.td-p2*c3.td, min.pr = VC$min.pr, max.pr = VC$max.pr)^-2*(dS1eta1*c2.td)*(dS2eta2*c3.td) ))*dereta2derb1, dereta1derb1) +
crossprod(c(VC$weights*VC$indvIT*( p3^-2*dS3eta3^2 - p3^-1*d2S3eta3 ))*dereta3derb1, dereta3derb1 ) +
diag( colSums( t( t(-c(VC$weights*VC$indvIT*p3^-1*dS3eta3 )*VC$X3)*der2.par1 ) ) )
)
############################################################################################################################################################################################
if(VC$extra.regI == "pC") H <- regH(H, type = 1)
S.h <- ps$S.h + monP2 # hess
S.h1 <- 0.5*crossprod(params, ps$S.h)%*%params + monP # lik
S.h2 <- S.h%*%params + monP1 # grad
S.res <- res
res <- S.res + S.h1
G <- G + S.h2
H <- H + S.h
if(VC$extra.regI == "sED") H <- regH(H, type = 2)
list(value = res, gradient = G, hessian = H, S.h = S.h, S.h1 = S.h1, S.h2 = S.h2, indN = indN, V = V,
l = S.res, l.ln = l.ln, l.par = l.par, ps = ps, k = VC$my.env$k, monP2 = monP2, params1 = params1,
eta1 = eta1, p1 = p1, p2 = p2, p3 = p3,
pdf1 = -dS1eta1, pdf2 = -dS2eta2, pdf3 = -dS3eta3,
dl.dbe1 = NULL, dl.dbe2 = NULL, dl.dteta.st = NULL,
hrate = VC$hrate,
d.lchrate = VC$d.lchrate,
d.rchrate = VC$d.rchrate,
d.lchrate.td = VC$d.lchrate.td,
d.rchrate.td = VC$d.rchrate.td)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.