R/ll_flexrsurv_beta0beta_bh.R

In flexrsurv: Flexible Relative Survival Analysis

ll_flexrsurv_beta0beta_bh<-function(beta0beta, alpha, gamma0, alpha0,
		Y, X0, X, Z, 
		expected_rate,
		weights=NULL,
		step, Nstep, 
		intTD=intTD_NC, intweightsfunc=intweights_CAV_SIM,
		nT0basis,
		Spline_t0=BSplineBasis(knots=NULL,  degree=3,   keep.duplicates=TRUE), Intercept_t0=TRUE, 
		nX0,
		ibeta0, nX, 
		ibeta, nTbasis,
		Spline_t =BSplineBasis(knots=NULL,  degree=3,   keep.duplicates=TRUE),
		Intercept_t_NPH=FALSE,
		debug=FALSE, ...){
	# compute log likelihood of the relative survival  model
	# rate = f(t)%*%gamma * exp(  X0%*%alpha0 + X%*%beta0(t) + sum( alphai(zi)betai(t) ))
	# case where nX > 0 ie, there are NPH but LIN effects 
	# gamma : vector of coef for baseline hazard
	# alpha0 ; vector of all coefs for non time dependant variables (may contain non-loglinear terms such as spline)
	# beta0beta ; matrix of all coefs for time dependant variables  X%*%beta0(t) and alpha(Z)%*%beta(t)
	#          beta0 = beta0beta[ibeta0] and beta = beta0beta[ibeta]
	# beta does not contains coef for the first t-basis
	#################################################################################################################
	# alpha : vector of coef for alpha(z) for NLG and NPH
	# Y : object of class Surv
	# 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 "DesignMatrixNPHNLL" time dependent variables (spline basis expended)
	# expected_rate : expected rate at event time T
	# weights : vector of weights  : LL = sum_i w_i ll_i
	# step : lag of subinterval for numerical integration fr each observation
	# Nstep : number of lag for each observation
	# intTD : function to perform numerical integration 
	# intweightfunc : function to compute weightsfor numerical integration
	# nT0basis : number of spline basis 
	#  Spline_t0, spline object for baseline hazard, with evaluate() method
	#  Intercept_t0=FALSE, option for evaluate, = TRUE all the basis, =FALSE all but first basis
	
	# nTbasis : number of time spline basis for NPH or NLL effects
	# nX0   : nb of PH variables dim(X0)=c(nobs, nX0)
	# nX    : nb of NPHLIN variables dim(X)=c(nobs, nX)
	#  Spline_t, spline object for time dependant effects,  with evaluate() method
	# Intercept_t_NPH vector of intercept option for NPH spline (=FALSE when X is NLL too, ie in case of remontet additif NLLNPH)
	#  ... not used args
	# returned value : the log liikelihood of the model
	
	if(is.null(Z)){
		nZ <- 0
	} else
	{
		nZ <- Z@nZ
	}
	
	
	if(Intercept_t0){
		tmpgamma0 <- gamma0
	}
	else {
		tmpgamma0 <- c(0, gamma0)
	}
	
	# baseline hazard at the end of the interval
	
	YT0Gamma0 <- predictSpline(Spline_t0*tmpgamma0, Y[,1], intercept=Intercept_t0)
	
	
	
	# contribution of non time dependant variables
	if(nX0) PHterm <-exp(X0 %*% alpha0)
	else PHterm <- 1
	
	# contribution of time d?pendant effect
	# parenthesis are important for efficiency
	if(nZ) {
		Zalphabeta <- Z@DM %*% ( ( diag(alpha)%*% Z@signature ) %*% t(ExpandAllCoefBasis(beta0beta[ibeta], ncol=nZ,  value=1)))
		if(nX) {
			Zalphabeta <- Zalphabeta + X %*% t(ExpandCoefBasis(beta0beta[ibeta0],
							ncol=nX,
							splinebasis=Spline_t,
							expand=!Intercept_t_NPH,
							value=0))
		}
	}
	else {
		if(nX) {
			Zalphabeta <- X %*% t(ExpandCoefBasis(beta0beta[ibeta0],
							ncol=nX,
							splinebasis=Spline_t,
							expand=!Intercept_t_NPH,
							value=0))
		}
		else {
			Zalphabeta <- NULL
		}
	}
	
	
	if( nX+nZ ){
		NPHterm <- intTD(rateTD_bh_alphabeta, intTo=Y[,1], intToStatus=Y[,2],
				step, Nstep,
				intweightsfunc=intweightsfunc, 
				gamma0=gamma0, Zalphabeta=Zalphabeta, 
				Spline_t0=Spline_t0*tmpgamma0, Intercept_t0=Intercept_t0,
				Spline_t = Spline_t, Intercept_t=TRUE)
	}
	else {
#    NPHterm <- intTD(rateTD_gamma0_bh, intTo=Y[,1], intToStatus=Y[,2],
#                     step=step, Nstep=Nstep,
#                     intweightsfunc=intweightsfunc, 
#                     gamma0=gamma0,
#                     Spline_t0=Spline_t0*tmpgamma0, Intercept_t0=Intercept_t0)
		NPHterm <- predict(integrate(Spline_t0*tmpgamma0), Y[,1], intercep=Intercept_t0)
	}
	# spline bases for baseline hazard
	YT0 <- evaluate(Spline_t0, Y[,1], intercept=Intercept_t0)
	if(nX0){
		if(nX + nZ){
			# spline bases for each TD effect
			YT <- evaluate(Spline_t, Y[,1], intercept=TRUE)
			eventterm <- ifelse(Y[,2] ,
					log( PHterm * (YT0Gamma0) * exp( apply(YT * Zalphabeta, 1, sum)) + expected_rate ), 
					0)
		} 
		else {
			eventterm <- ifelse(Y[,2] , 
					log( PHterm * (YT0Gamma0) + expected_rate ), 
					0)
		}
	} 
	else {
		if(nX + nZ){
			# spline bases for each TD effect
			YT <- evaluate(Spline_t, Y[,1], intercept=TRUE)
			eventterm <- ifelse(Y[,2] ,
					log( (YT0Gamma0) * exp(apply(YT * Zalphabeta, 1, sum)) + expected_rate ), 
					0)
		} else {
			eventterm <- ifelse(Y[,2] , 
					log( (YT0Gamma0) + expected_rate ), 
					0)
		}
	}
	
	
	if (!is.null(weights)) {
		if( nX0){
			ret <- crossprod(eventterm - PHterm * NPHterm , weights)
		} else {
			ret <- crossprod(eventterm - NPHterm , weights)
		}
	}
	else {
		if( nX0){
			ret <- sum( eventterm - PHterm * NPHterm )
		} else {
			ret <- sum( eventterm - NPHterm )
		}
	}
	
	ret
	
}