Nothing
flexrsurv.ll.fromto.brass0.wceadd.fitCoptim<-function (Y, X0, X, Z, W,
BX0,
Id, FirstId, LastId,
isEnter,
isEnd,
expected_rate,
expected_logit_end,
expected_logit_enter,
weights=NULL,
Ycontrol=NULL, BX0control,
weightscontrol,
Idcontrol, FirstIdcontrol, LastIdcontrol,
isEntercontrol,
isEndcontrol,
expected_ratecontrol,
expected_logit_endcontrol,
expected_logit_entercontrol,
Spline_t =BSplineBasis(knots=NULL, degree=3, keep.duplicates=TRUE), Intercept_t_NPH=TRUE,
ISpline_W =MSplineBasis(knots=NULL, degree=3, keep.duplicates=TRUE), Intercept_W=TRUE,
Spline_B =LEBSplineBasis(knots=NULL, degree=3, keep.duplicates=TRUE), Intercept_B=TRUE,
init=list(eta0=NULL, alpha0=NULL, beta0=NULL, alpha=NULL, beta=NULL,brass0=NULL, balpha0=NULL),
optim.control=list(trace=100, REPORT=1, fnscale=-1, maxit=25),
Coptim.control=list(),
method=list(int_meth="CAV_SIM", step=diff(range(Y[,ncol(Y)-1]))/100, optim_meth="BFGS",
constOptim_meth="BFGS", lower=-Inf, upper = Inf),
vartype = c("oim", "opg", "none"),
varmethod = c("optim", "numDeriv.hessian", "numDeriv.jacobian"),
numDeriv.method.args=list(eps=5e-7, d=0.001, zero.tol=sqrt(.Machine$double.eps/7e-4), r=4, v=2),
namebrass="CorTable",
debug=FALSE, debug.ll=FALSE, debug.gr=FALSE,
...
) {
# flexible relative survival model using full likelihood and
# one WCE additif proportional with NLL, NPH, NPHNLL effects
# no baseline hazard, the baseline is the WCE effet
# the WCE link is identity (ie WCE(t) = sum b_j B_j(t-Tk) X_k
# corection of lifetable according to generalized brass method
# Cohort-independent generalized Brass model in an age-cohort table
# stratified brass model eccording to fixed effects BX0 (one brass function per combination)
# rate = brass0(expected-rate, expected_logit)*exp(BX0 balpha0) + WCE(W,t)*exp{ time-independent-effect(LL + NLL)(X0) + NPH(X) + NPHNLL(Z)}
#
# input :
#
#
# Y : object of class Surv but the matrix has 4 columns : beginning(1) and end(2) of interval, status(3) and end of followup(4)
# end of followup is assumed constant by Id
#################################################################################################################
# Y : object of class Surv but the matrix has 4 columns :
# Y[,1] beginning(1) , fromT, start
# Y[,2] end(2), toT, time
# Y[,3] status(3) fail, status
# Y[,4] end of followup(4) FinalT
# end of followup is assumed constant by Id
#################################################################################################################
# X0 : non-time dependante variable (may contain spline bases expended for non-loglinear terms)
# X : log lineair but time dependante variable
# Z : object of class time d?pendent variables (spline basis expended)
# W : Exposure variable used in Weighted Cumulative Exposure Models
# BX0 : non-time dependante variable (may contain spline bases expended for non-loglinear terms) for the Brass model
# Id : varibale indicating individuals Id, lines with the same Id are considered to be from the same individual
# FirstId : all lines in FirstId[iT]:iT in the data comes from the same individual
# LastId : all lines in FirstId[iT]:LastId[iT] in the data comes from the same individual Id[iT]
# expected_rate : expected rate at event time T
# expected_logit : expected logit at event time T (used in the Brass model
# weights : vector of weights : LL = sum_i w_i ll_i
# Spline_t, spline object for time dependant effects, with evaluate() method
# Intercept_t=FALSE, option for evaluate, = TRUE all the basis, =FALSE all but first basis
# ISpline_W, list of one integrated spline objects for the WCE effect, with evaluate() method
# Intercept_W=TRUE, option for evaluate the ith wce effect, = TRUE all the basis, =FALSE all but first basis
# Spline_B, spline objects for cohort-independent Brass function, with evaluate() method
# Intercept_B=FALSE, option for evaluate the brass function = TRUE all the basis, =FALSE all but first basis
# init : list of initial values
# optime.control : control parameters/options for optim()
# method : optimisation method (optim_meth) for optim(), numerical integration method (int_meth),
# vartype : type of variance matrix : observed inf. mat (oim inv(-H)), robust/sandwich (robust H inv(S'S) H ),
# outer product of the gradients (opg inv(S'S)), wher where S is the matrix of scores
# namebrass="CorrectionTable" : used to build the names of the parameters of the brass function
#
# output : coef(with name, structure as attributes), var, conv & LL
# coef=c(eta0, alpha0, beta0, beta, alpha, brass0, balpha0 )
# with
# eta0 : vector of all the coef for the WCE effects
# alpha0 ; vector of all coefs for non time dependant variables (may contain non-loglinear terms such as spline)
# beta0 ; matrix of all coefs for log-linear but time dependant variables X%*%beta0(t)
# beta : matrix of coefs for beta(t) nTbasis * nTDvars for NLG and NPH
# alpha : vector of coef for alpha(z) for NLG and NPH
# brass0 : vector of coef for the cohort independent Brass function
# balpha0 ; vector of all coefs for non time dependant variables (may contain non-loglinear terms such as spline) of the brass model
# for correcting life table
#
vartype <- match.arg(vartype) # type of var
varmethod <- match.arg(varmethod) # function for computing oim
if( vartype == "opg"){
stop("vartype == 'opg' not implemented for WCE effect", call.=TRUE)
}
dohessian <- ( vartype=="oim") & (varmethod == "optim")
# if (ncol(Y)==4) {
# start <- Y[,1]
# time <- Y[,2]
# status <- Y[,3]
# finaltime <- Y[,4]
# }
# else {
# stop("wrong dim of Y", call.=TRUE)
# }
start <- Y[,1]
time <- Y[,2]
status <- Y[,3]
n <- nrow(Y)
# structure of parameter vector
# WCE effects
if(!is.null(W)){
is.WCE <- TRUE
if(is.matrix(W)) {
nW <- dim(W)[2]
if(nW>1){
stop("wrong dim of W, only one WCE effect allowed", call.=TRUE)
}
W <- as.vector(W)
} else if(is.vector(W)) {
} else {
stop("wrong dim of W", call.=TRUE)
}
nW <- 1L
nWbasis <- getNBases(ISpline_W) - 1 + Intercept_W
neta0 <- sum(nWbasis)
ieta0 <- 1:neta0
Ieta0 <- 1:neta0
First.eta0 <- 1
first.eta0 <- 1
iWend <- cumsum(nWbasis)
iWbeg <- c(1, iWend[-nW]-1)
} else {
# no WCE, constant effect WCE(t) = 1, excess model = exp(LIN+NPH+NPHNLL)
is.WCE <- FALSE
nWbasis <- 0L
nW <- 0L
neta0 <- 0L
eta0 <- NULL
ieta0 <- NULL
Ieta0 <- NULL
First.eta0 <- NULL
first.eta0 <- NULL
iWend <- NULL
iWbeg <- NULL
}
# PHLIN and PHNLIN effects
if(!is.null(X0)){
is.PH <- TRUE
if(is.matrix(X0)) {
nX0 <- dim(X0)[2]
} else if(is.vector(X0)) {
nX0 <- 1L
} else {
stop("wrong type of X0", call.=TRUE)
}
nalpha0<-nX0
ialpha0<-1:nX0 + neta0
Ialpha0<-1:nX0
First.alpha0<-neta0+1
first.alpha0<-1
} else {
is.PH <- FALSE
nX0 <- 0L
nalpha0 <- 0L
alpha0 <- NULL
ialpha0 <- NULL
Ialpha0 <- NULL
First.alpha0 <- NULL
first.alpha0 <- NULL
}
# NPHLIN effects
if(!is.null(X)){
is.NPHLIN <- TRUE
nTbasis <- getNBases(Spline_t)
nTbasis_NPH <- getNBases(Spline_t) - 1 + Intercept_t_NPH
if(is.matrix(X)) {
nX <- dim(X)[2]
} else if(is.vector(X)) {
nX <- 1L
} else {
stop("wrong type of X", call.=TRUE)
}
nbeta0 <- sum(nTbasis_NPH)
ibeta0 <- 1:nbeta0 + neta0 + nalpha0
Ibeta0 <- 1:nbeta0
First.beta0 <- neta0 + nalpha0 + 1
first.beta0 <- 1
} else {
is.NPHLIN <- FALSE
nTbasis <- 0L
nTbasis_NPH <- 0L
nX <- 0L
nbeta0 <- 0L
beta0 <- NULL
ibeta0 <- NULL
Ibeta0 <- NULL
First.beta0 <- NULL
first.beta0 <- NULL
}
# NPHNLIN effects
if(!is.null(Z)) {
if( getNvar(Z)>0){
is.NPHNLIN <- TRUE
nTbasis <- getNBases(Spline_t)
nZ<-getNvar(Z)
nTbasis_NPHNLL <- rep(getNBases(Spline_t), nZ)
nalpha <- getNparam(Z)
ialpha <- 1:nalpha + neta0 + nalpha0 + nbeta0
Ialpha <- 1:nalpha + nalpha0
First.alpha <- neta0 + nalpha0 + nbeta0 + 1
first.alpha <- nalpha0 + 1
# as first beta is constraints, nTbasis -1 beta per Z variable
nbeta <- nZ * (nTbasis-1)
ibeta <- 1:nbeta + neta0 + nalpha0 + nbeta0 + nalpha
Ibeta <- 1:nbeta + nbeta0
First.beta <- neta0 + nalpha0 + nbeta0 + nalpha + 1
first.beta <- 1 + nbeta0
}
} else {
is.NPHNLIN <- FALSE
# nTbasis <- 0 already done with beta0
nZ <- 0L
nalpha <- 0L
alpha <- NULL
Ialpha <- NULL
ialpha <- NULL
First.alpha <- NULL
first.alpha <- NULL
nbeta <- 0L
beta <- NULL
Ibeta <- NULL
ibeta <- NULL
First.beta <- NULL
first.beta <- NULL
}
# Brass model
# Brass function
if(!is.null(Spline_B)){
nBbasis <- dim(evaluate(Spline_B, expected_logit_end, intercept=Intercept_B))[2]
# first basis is the linear basis, with coef ==1
nbrass0 <- nBbasis -1
ibrass0<-1:nbrass0 + neta0 + nalpha0 + nbeta0 + nalpha + nbeta
Ibrass0<-1:nbrass0
First.brass0<-neta0 + nalpha0 + nbeta0 + nalpha + nbeta + 1
first.brass0<-1
}
else {
# no Brass function
nBbasis <- 0
nbrass0 <- 0
ibrass0<-NULL
Ibrass0<-NULL
First.brass0<-NULL
first.brass0<-NULL
}
# Proportional corrections
if(!is.null(BX0)){
is.BPH <- TRUE
if(is.matrix(BX0)) {
nBX0 <- dim(BX0)[2]
} else if(is.vector(BX0)) {
nBX0 <- 1L
} else {
stop("wrong type of X0", call.=TRUE)
}
nbalpha0<-nBX0
ibalpha0<-1:nBX0 + neta0 + nalpha0 + nbeta0 + nalpha + nbeta + nbrass0
Ibalpha0<-1:nBX0
First.balpha0<-neta0 + + nalpha0 + nbeta0 + nalpha + nbeta + nbrass0 + 1
first.balpha0<-1
} else {
is.PH <- FALSE
nBX0 <- 0L
nbalpha0 <- 0L
balpha0 <- NULL
ibalpha0 <- NULL
Ibalpha0 <- NULL
First.balpha0 <- NULL
first.balpha0 <- NULL
}
# number of variables
nvar <- nX0 + nX + nZ +nW
# nb of parameters
nparam = neta0 + nalpha0 + nbeta0 + nalpha + nbeta + nbrass0 + nbalpha0
LastId <- getLastId(FirstId)
if(!is.null(Ycontrol)) {
LastIdcontrol <- getLastId(FirstIdcontrol)
}
# default init values for eta0 (same value at each basis)
# contribution to the cumulatice WCE of each line (if eta==1)
# it is NPHterm of the ll_function when I_Spline in unscaled and no time dependent terms
wce2 <- predictwce(object=integrate(ISpline_W), t=Y[,2], Increment=W, fromT=Y[,1], tId=(1:dim(Y)[1]),
FirstId=FirstId, LastId=LastId, intercept=Intercept_W, outer.ok=TRUE)
wce1 <- predictwce(object=integrate(ISpline_W), t=Y[,1], Increment=W, fromT=Y[,1], tId=(1:dim(Y)[1]),
FirstId=FirstId, LastId=LastId, intercept=Intercept_W, outer.ok=TRUE)
cumwcebasis <- wce2 - wce1
cumwcebasis2 <- predictwce(object=integrate(ISpline_W), t=Y[Y[,3]==1,2], Increment=W, fromT=Y[,1], tId=(1:dim(Y)[1])[Y[,3]==1],
FirstId=FirstId, LastId=LastId, intercept=Intercept_W, outer.ok=TRUE)
initeta0 <- rep(sum(status)/sum(cumwcebasis), neta0)
if (!missing(expected_rate) && !is.null(expected_rate)) {
WCEevent0 <- predictwce(object=ISpline_W, t=Y[,2], Increment=W, fromT=Y[,1], tId=(1:dim(Y)[1]),
FirstId=FirstId, LastId=LastId, intercept=Intercept_W, outer.ok=TRUE)
# keep valaues only at event
rateevent <-expected_rate[Y[,3]==1]
WCEevent0 <- WCEevent0[Y[,3]==1]
# initeta0 <- initeta0 - mean((rateevent/WCEevent0)[WCEevent0>0])
initeta0 <- initeta0 - mean(rateevent)
}
if (missing(expected_rate) || is.null(expected_rate))
expected_rate <- rep(0, n)
if (missing(weights) || is.null(weights)) {
weights <- NULL
} else {
if (any(weights <= 0))
stop("Invalid weights, must be >0", call.=TRUE)
storage.mode(weights) <- "double"
}
# control of init values
if (!missing(init) && !is.null(init)) {
if (nW ){
if (!is.null(init$eta0)) {
if (length(init$eta0) != neta0){
stop("Wrong length for initial values for eta0", call.=TRUE)
}
else {
initeta0 <- init$eta0
}
}
else {
# initeta0 is the default value
}
}
else {
initeta0 <- NULL
}
if ( nX0 ) {
if(!is.null(init$alpha0)) {
if (length(init$alpha0) != nX0 ) {
stop("Wrong length for initial values for alpha0", call.=TRUE)
} else {
initalpha0 <- init$alpha0
}
} else {
initalpha0 <- rep(0, nX0)
}
} else {
initalpha0 <- NULL
}
if (nX ){
if (!is.null(init$beta0)) {
if (length(init$beta0) != nbeta0){
stop("Wrong length for initial values for beta0", call.=TRUE)
} else {
initbeta0 <- init$beta0
}
} else {
initbeta0 <- rep(0, nbeta0)
}
} else {
initbeta0 <- NULL
}
if (nZ ){
if (!is.null(init$alpha)) {
if (length(init$alpha) != nalpha ){
stop("Wrong length for initial values for alpha", call.=TRUE)
} else {
initalpha <- init$alpha
}
} else {
initalpha <- rep(0, nalpha)
}
if (!is.null(init$beta)) {
if (length(init$beta) != nbeta) {
stop("Wrong length for initial values for beta", call.=TRUE)
} else {
initbeta <- init$beta
}
} else {
# else initbeta=0
initbeta <- rep(0, nbeta)
}
}
else {
initalpha <- NULL
initbeta <- NULL
}
if (nbrass0) {
if (!is.null(init$brass0)) {
if (length(init$brass0) != nbrass0){
stop("Wrong length for initial values for brass0", call.=TRUE)
}
else {
initbrass0 <- init$brass0
}
}
else {
initbrass0 <- rep(0, nbrass0)
}
} else {
initbrass0 <- NULL
}
if (nBX0) {
if (!is.null(init$balpha0)) {
if (length(init$balpha0) != nbalpha0){
stop("Wrong length for initial values for balpha0", call.=TRUE)
}
else {
initbalpha0 <- init$balpha0
}
}
else {
initbalpha0 <- rep(0, nbalpha0)
}
} else {
initbalpha0 <- NULL
}
# end if (!missing(init) && !is.null(init))
} else {
# user did not provide init values
cat("default init values used\n")
# no init values
if (nW) {
# initeta0 is the default value (see above)
} else {
initeta0 <- NULL
}
if ( nX0) {
initalpha0 <- rep(0, nX0)
} else {
initalpha0 <- NULL
}
if (nX) {
initbeta0 <- rep(0, nbeta0)
} else {
initbeta0 <- NULL
}
if (nZ) {
initalpha <- rep(0, nalpha)
initbeta <- rep(0 , nbeta)
} else {
initalpha <- NULL
initbeta <- NULL
}
if (nbrass0) {
initbrass0 <- rep(0, nbrass0)
} else {
initbrass0 <- NULL
}
if (nBX0) {
initbalpha0 <- rep(0, nbalpha0)
} else {
initbalpha0 <- NULL
}
}# end else if (!missing(init) && !is.null(init))
storage.mode(initalpha0) <- "double"
storage.mode(initbeta0) <- "double"
storage.mode(initalpha) <- "double"
storage.mode(initbeta) <- "double"
storage.mode(initeta0) <- "double"
storage.mode(initbrass0) <- "double"
storage.mode(initbalpha0) <- "double"
# numerical integration method
# computes steps for time integtration
if(method$int_meth == "CAV_SIM"){
int_meth <- "NC"
intTD <- intTDft_NC
intTD_debug<- intTDft_NC_debug
intTD_base<- intTDft_base_NC
intTD_WCEbase<- intTDft_WCEbase_NC
intweightsfunc <-intweights_CAV_SIM
step <-method$step
mult <- 2
} else if(method$int_meth == "SIM_3_8"){
int_meth <- "NC"
intTD <- intTDft_NC
intTD_base<- intTDft_base_NC
intTD_WCEbase<- intTDft_WCEbase_NC
intweightsfunc <-intweights_SIM_3_8
step <-method$step
mult <- 3
} else if(method$int_meth == "BOOLE"){
int_meth <- "NC"
intTD <- intTDft_NC
intTD_base<- intTDft_base_NC
intTD_WCEbase<- intTDft_WCEbase_NC
intweightsfunc <-intweights_BOOLE
step <-method$step
mult <- 4
} else if(method$int_meth == "Gauss-Legendre"){
int_meth <- "GL"
intTD <- intTDft_GL
intTD_base<- intTDft_base_GL
intTD_WCEbase<- intTDft_WCEbase_GL
intweightsfunc <-NULL
gq <- gauss.quad(method$npoints, kind="legendre")
step <- gq$nodes
Nstep <- gq$weights
} else if(method$int_meth == "GLM"){
int_meth <- "GLM"
intTD <- intTDft_GLM
intTD_base<- fastintTDft_base_GLM
intTD_WCEbase<- fastintTDft_WCEbase_GLM
intweightsfunc <- NULL
step <-GLMStepParam(cuts=method$bands)
Firststep <- WhichBandInf(start, step) + 1L
Laststep <- WhichBand(time, step)- 1L
Nstep <- cbind(Firststep, Laststep)
} else {
stop("unknown method$int_meth", call.=TRUE)
}
if( int_meth == "NC"){
STEPS<-cutTfromto(start, time, step=method$step, mult=mult)
Nstep<-STEPS$Nstep
step<-STEPS$step
}
LL <- +Inf
# objective, gradiant functions
ll_function <- ll_flexrsurv_fromto_1WCEaddBr0Control
gr_function <- gr_ll_flexrsurv_fromto_1WCEaddBr0Control
# gr_function <- NULL
opgFunction <- opg_flexrsurv_fromto_1WCEaddBr0Control
computeLinearPredictor <- .computeLinearPredictor_fromto_1wceadd
computeCumulativeHazard <- .computeCumulativeHazard_fromto_1wceadd
alpha0 <- initalpha0
alpha <- initalpha
beta0 <- initbeta0
beta <- initbeta
eta0 <- initeta0
brass0 <- initbrass0
balpha0 <- initbalpha0
initGA0B0ABE0Br0 <- c(eta0, alpha0, beta0, alpha, beta, brass0, balpha0)
if(debug) cat("********************get init LL values \n")
LL <- ll_function(allparam=initGA0B0ABE0Br0,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol = LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0=ibeta0, nX=nX,
ialpha=ialpha,
ibeta=ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
debug=debug.ll
)
if (debug) {
print("*************** init values")
print(initGA0B0ABE0Br0)
cat("LL at init", LL, "\n")
}
LLold<- LL
conv <- TRUE
iter <- -1
algo <- "optim"
if(!is.null(gr_function)){
GRInit <- gr_function(allparam=initGA0B0ABE0Br0,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0=ibeta0, nX=nX,
ialpha=ialpha,
ibeta=ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug.gr=debug.ll
)
# method="L-BFGS-B",
# method="BFGS",
# lower=lowerGA0B0ABE0Br0,
# upper=upperGA0B0ABE0Br0,
# method="Nelder-Mead",
if (debug) {
cat("GR at init", "\n")
print(cbind(initGA0B0ABE0Br0, GRInit))
cat("GR at init", "\n")
}
}
else {
GRInit <- NULL
}
if(!is.null(Spline_B)) {
# constraint optimisation
# for brass model : deriv(Spline_B) > 0
# first basis is contraints to one
# constraits is betaclt %*% evaluate(deriv(Spline_B), PivotConstraintSplinceCLT.R2bBSplineBasis(Spline_B)[-1]) > 0
# constraints on brass parameters : brass0[1] unconstraints (la constante)
# brass0[-1] >= 0
# other parameters unconstraines
# dim(ui) = c( nbrass0-1, nparam)
# ui[, 1: First.brass0-1] == 0
# ui[, First.brass0 -1+(1:nbrass0)] == evaluate(deriv(Spline_B), pivot)[-1, -1]
# ci = evaluate(deriv(Spline_B), pivot)[-1, 1]
# ui[, (First.brass0+nbrass0):nparam] = 0
pivot <- PivotConstraintSplinceCLT.R2bBSplineBasis(Spline_B)
matder <- evaluate(deriv(Spline_B), pivot)
ui <- matder[-1, -1]
ci <- - matder[-1, 1]
if(First.brass0>1){
ui <- cbind(matrix(0, ncol = First.brass0-1, nrow=nbrass0), ui)
}
if(!is.null(BX0)){
ui <- cbind(ui, matrix(0, ncol = nBX0 , nrow=nbrass0))
}
# ui <- diag(nbrass0-1)
# ci = rep(0, nbrass0-1)
# if(First.brass0>1){
# ui <- cbind(matrix(0, ncol = First.brass0, nrow=nbrass0-1), ui)
# }
# if(!is.null(BX0)){
# ui <- cbind(ui, matrix(0, ncol = nBX0 , nrow=nbrass0-1))
# }
# unconstraint optimisation with optim
# fit<-optim(initGA0B0ABE0Br0,
# fn=ll_function,
# gr = gr_function,
# method=method$optim_meth,
# constraint optimisation with constrOptim
# fit<-constrOptim(initGA0B0ABE0Br0,
# f=ll_function,
# grad = gr_function,
# ui=ui,
# ci=rep(0, nbrass0-1),
# method="L-BFGS-B",
## method="Nelder-Mead",
# constraints=NULL,
# lower=lowerGA0B0ABE0Br0,
# upper=upperGA0B0ABE0Br0,
Con <- list(outer.iterations = 100, outer.eps = 1e-05, mu = 1e-04)
nCon <- names(Con)
Con[(nCoptim.control <- names(Coptim.control))] <- Coptim.control
Con <- Con[nCon]
#print("************************************************** fit")
#print(initGA0B0ABE0Br0)
#print(ui)
#print("************************************************** fit")
#print(optim.control)
options(show.error.messages = FALSE)
fit<-try(constrOptim(initGA0B0ABE0Br0,
f=ll_function,
grad = gr_function,
# grad = NULL,
ui=ui,
ci=rep(0, nbrass0-1),
method=method$constOptim_meth,
lower=method$lower,
upper=method$upper,
constraints=NULL,
outer.iterations = Con$outer.iterations, outer.eps = Con$outer.eps, mu = Con$mu,
control = optim.control,
hessian = dohessian,
# ll_flexrsurv_fromto_1WCEaddBr0Control args
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0= ibeta0, nX=nX,
ialpha=ialpha,
ibeta= ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug=debug.ll,
debug.gr=debug.gr),
silent=TRUE)
options(show.error.messages = TRUE)
#print("fin fit")
}
else {
#print("####################################### # no Brass correction, unconstrait optimisation")
# no Brass correction, unconstrait optimisation
options(show.error.messages = FALSE)
fit<-try(optim(initGA0B0ABE0Br0,
fn=ll_function,
gr = gr_function,
# gr = NULL,
method=method$optim_meth,
lower=method$lower,
upper=method$upper,
constraints=NULL,
# ll_flexrsurv_fromto_GA0B0ABE0Br0 args
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0= ibeta0, nX=nX,
ialpha=ialpha,
ibeta= ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
control = optim.control,
debug=debug.ll,
debug.gr=debug.gr,
hessian = dohessian),
silent=TRUE)
options(show.error.messages = TRUE)
}
if( inherits(fit, "try-error")){
msg <- fit
cat(msg)
res <- list(coefficients = rep(NA, ),
list_coefficients = NULL,
linear.predictors=NULL,
fitted.values=NULL,
cumulative.hazard=NULL,
var = numeric(0),
informationMatrix= numeric(0),
loglik = NA,
weights = weights,
optimfunction=algo,
conv=msg,
method = "flexrsurv.ll.fromto.brass.wce.fit",
fit=fit,
start=list(eta0=initeta0, alpha0=initalpha0, beta0=initbeta0, alpha=initalpha, beta=initbeta,
brass0=initbrass0, balpha0 = initbalpha0))
res
} else {
#print("LLfinal")
LLfinal <- ll_function(allparam=fit$par,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0=ibeta0, nX=nX,
ialpha=ialpha,
ibeta=ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug=debug.ll
)
#print("fin LLfinal")
if(!is.null(gr_function)){
#print("GRfinal")
GRfinal <- gr_function(allparam=fit$par,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0=ibeta0, nX=nX,
ialpha=ialpha,
ibeta=ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug.gr=debug.ll
)
#print("finGRfinal")
}
else {
GRfinal <- NULL
}
print("********************************* fin")
print(c(LLfinal, fit$value, LLfinal - fit$value))
print(cbind(fit$par, GRfinal, GRInit))
print("********************************* fin")
conv <- converged(fit, step="maximisation eta0 alpha0 beta0 alpha and beta")
#print("message conv")
#print(conv)
if(nW){
eta0<-fit$par[ieta0]
}
if (nX0){
alpha0 <- fit$par[ialpha0]
}
if(nX){
beta0<-fit$par[ibeta0]
}
if (nZ){
alpha <- fit$par[ialpha]
beta <- fit$par[ibeta]
}
if (nbrass0){
brass0 <- fit$par[ibrass0]
}
if (nBX0){
balpha0 <- fit$par[ibalpha0]
}
LL <- fit$value
iter <- NA
# prepare returned value
if(nW){
nameseta0 <- NULL
for(iw in 1:nW) {
nameseta0 <- c(nameseta0,
paste(paste("WCE(",
dimnames(W)[[2]][iw],
", ",
dimnames(Y)[[2]][ncol(Y)],
"):",
sep=""),
(1-Intercept_W[iw])+1:nWbasis[iw],
sep="")
)
}
names(eta0) <- nameseta0
}
if (nX0){
names(alpha0) <- dimnames(X0)[[2]]
}
if(nX){
namesbeta0 <- NULL
for(ix in 1:nX) {
namesbeta0 <- c(namesbeta0,
paste(paste("NPH(",
dimnames(X)[[2]][ix],
", ",
dimnames(Y)[[2]][ncol(Y)-1],
"):",
sep=""),
(2-Intercept_t_NPH[ix]):nTbasis,
sep="")
)
}
names(beta0) <- namesbeta0
}
if (nZ){
names(alpha) <-paste("NPHNLL:", dimnames(getDesignMatrix(Z))[[2]], sep="")
namesbeta <- NULL
for(iz in 1:nZ) {
namesbeta <- c(namesbeta,
paste(paste("NPHNLL(", getNames(Z), sep=""),
paste(dimnames(Y)[[2]][ncol(Y)-1],
"):",
2:nTbasis,
sep=""),
sep="_"))
names(beta) <- namesbeta
}
}
if (nbrass0){
namesbrass0 <- paste(namebrass, ":",
(1-Intercept_B)+1:nbrass0,
sep="")
names(brass0) <- namesbrass0
}
if (nBX0){
names(balpha0) <- paste(namebrass, ":", dimnames(BX0)[[2]], sep="")
}
lcoef <- list(eta0=eta0, alpha0=alpha0, beta0=beta0 , alpha=alpha, beta=beta, brass0=brass0, balpha0=balpha0)
coef <- c(lcoef$eta0, lcoef$alpha0, lcoef$beta0 , lcoef$alpha, lcoef$beta, lcoef$brass0, lcoef$balpha0)
# variance computation
#print("# variance computation")
var <- NULL
informationMatrix <- NULL
cholinformationMatrix <- NULL
if( vartype == "oim"){
# the variance matrix is obtained from the observed information matrix
# numericNHessian in package maxLik
var <- NULL
if( varmethod == "optim" ){
#print("varmethod = optim")
informationMatrix <- - fit$hessian
#print("varmethod = optim 2")
} else if( varmethod == "numDeriv.hessian" ){
#print("varmethod == numDeriv.hessian")
informationMatrix <- - hessian(ll_function, fit$par,
method.args=numDeriv.method.args,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0=ibeta0, nX=nX,
ialpha=ialpha,
ibeta=ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug=debug.ll)
#print("varmethod == numDeriv.hessian2")
} else if( varmethod == "numDeriv.jacobian" ){
#print("varmethod == numDeriv.jacobian")
informationMatrix <- - jacobian(gr_function, fit$par,
method.args=numDeriv.method.args,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0=ibeta0, nX=nX,
ialpha=ialpha,
ibeta=ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug=debug.ll)
#print("varmethod == numDeriv.jacobian2")
informationMatrix <- (informationMatrix + t(informationMatrix))/2
#print("varmethod == numDeriv.jacobian3")
}
} else if( vartype == "opg"){
# the variance matrix is obtained from the outer product of the gradient
# in the case of Type 1 censoring
informationMatrix <- opgFunction(allparam=fit$par,
Y=Y, X0=X0, X=X, Z=Z, W=W,
BX0=BX0,
Id=Id, FirstId=FirstId, LastId=LastId,
isEnter=isEnter,
isEnd=isEnd,
expected_rate=expected_rate,
expected_logit_enter=expected_logit_enter,
expected_logit_end=expected_logit_end,
weights = weights,
Ycontrol = Ycontrol, BX0control = BX0control,
weightscontrol = weightscontrol,
Idcontrol = Idcontrol, FirstIdcontrol = FirstIdcontrol, LastIdcontrol=LastIdcontrol,
isEntercontrol = isEntercontrol,
isEndcontrol = isEndcontrol,
expected_ratecontrol = expected_ratecontrol,
expected_logit_endcontrol = expected_logit_endcontrol,
expected_logit_entercontrol = expected_logit_entercontrol,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
intTD_base=intTD_base,
intTD_WCEbase=intTD_WCEbase,
ialpha0=ialpha0, nX0=nX0,
ibeta0= ibeta0, nX=nX,
ialpha=ialpha,
ibeta= ibeta,
nTbasis=nTbasis,
ieta0=ieta0, iWbeg=iWbeg, iWend=iWend, nW=nW,
ibrass0=ibrass0, nbrass0=nbrass0,
ibalpha0=ibalpha0, nBX0=nBX0,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
nBbasis=nBbasis,
Spline_B=Spline_B, Intercept_B=Intercept_B,
debug.gr=debug.gr
)
}
var <- NULL
if(vartype != "none"){
# first try inverting from its Choleski decomposition.
#print(" # first try inverting from its Choleski decomposition.")
cholinformationMatrix <- NULL
options(show.error.messages = FALSE)
cholinformationMatrix <- try(chol(informationMatrix), silent=TRUE)
#print(" cholinformationMatrix <- try(chol(informationMatrix), silent=TRUE)")
options(show.error.messages = TRUE)
if( !inherits(cholinformationMatrix, "try-error")){
options(show.error.messages = FALSE)
#print(" var <- try( chol2inv(cholinformationMatrix) , silent=TRUE) 1")
var <- try( chol2inv(cholinformationMatrix) , silent=TRUE)
#print(" var <- try( chol2inv(cholinformationMatrix) , silent=TRUE) 2")
options(show.error.messages = TRUE)
if( inherits(var, "try-error")){
var <- NULL
cat(geterrmessage())
}
else{
attr(var, "type") <- vartype
attr(var, "solve") <- "chol"
}
}
if( is.null(var)){
# second try inverting with solve()
#print(" # second try inverting with solve()")
options(show.error.messages = FALSE)
var <- try( solve(informationMatrix) , silent=TRUE)
#print(" # second try inverting with solve() 2")
options(show.error.messages = TRUE)
if( inherits(var, "try-error")){
var <- numeric(0)
cat(geterrmessage())
}
else{
attr(var, "type") <- vartype
attr(var, "solve") <- "solve"
}
}
}
attr(informationMatrix, "type") <- vartype
attr(var, "type") <- vartype
# computes the linear predictor relative to the excess hazard objfit$linear.predictor
# Remarque, make sure that .computeLinearPredictor_GA0B0ABE0 and .computeCumulativeHazard_fromto_GA0B0ABE0
# is able to get the wright parameters from ialpha0, ibeta0, ialpah and ibeta and iteat0,
# even if ther is parameters for brass model in fit$par
#print("linpred")
linearPredictors <- computeLinearPredictor(allparam=fit$par,
Y=Y, X0=X0, X=X, Z=Z, W=W,
Id=Id, FirstId=FirstId, LastId=LastId,
ialpha0=ialpha0, nX0=nX0,
ibeta0= ibeta0, nX=nX,
ialpha=ialpha,
ibeta= ibeta,
ieta0=ieta0,
nTbasis=nTbasis,
Spline_t = Spline_t,
Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
wcelink="identity",
debug=debug)
# computes the fited rate
#print(" # computes the fited rate")
fittedValues <- exp(linearPredictors[,1])*linearPredictors[,2]
# # computes the cumulative exxcess hazard ??
#print("# # computes the cumulative exxcess hazard ?? ")
cumulativeHazard <- computeCumulativeHazard(allparam=fit$par,
Y=Y, X0=X0, X=X, Z=Z, W=W,
Id=Id, FirstId=FirstId, LastId=LastId,
step=step, Nstep=Nstep,
intTD=intTD, intweightsfunc=intweightsfunc,
ialpha0=ialpha0, nX0=nX0,
ibeta0= ibeta0, nX=nX,
ialpha=ialpha,
ibeta= ibeta,
ieta0=ieta0,
nTbasis=nTbasis,
Spline_t = Spline_t, Intercept_t_NPH=Intercept_t_NPH,
ISpline_W = ISpline_W,
Intercept_W=Intercept_W,
debug=debug)
# convergence assesment
#print(" # convergence assesment")
if (conv != TRUE) {
warning("Ran out of iterations and did not converge")
}
res <- list(coefficients = coef,
list_coefficients = lcoef,
linear.predictors=linearPredictors,
fitted.values=fittedValues,
cumulative.hazard=cumulativeHazard,
var = var,
informationMatrix=informationMatrix,
loglik = LL,
weights = weights,
optimfunction=algo,
conv=conv,
method = "flexrsurv.ll.fromto.brass.wce.fit",
numerical_integration_method = method,
fit=fit,
start=list(eta0= initeta0, alpha0=initalpha0, beta0=initbeta0, alpha=initalpha, beta=initbeta,
brass0=initbrass0, balpha0 = initbalpha0))
#print("fin res")
res
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.