Nothing
# nobs - number of recurrent event
# nsujet - number of individuals
# nsujety - number of repeated measurements
#' Estimators of the Expected Prognostic Observed Cross-Entropy (EPOCE) for
#' evaluating predictive accuracy of joint models.
#'
#'
#' This function computes estimators of the Expected Prognostic Observed
#' Cross-Entropy (EPOCE) for evaluating the predictive accuracy of joint models
#' using \code{frailtyPenal}, \code{longiPenal}, \code{trivPenal} or
#' \code{trivPenalNL}. On the same data as used for estimation of the joint
#' model, this function computes both the Mean Prognosis Observed Loss (MPOL)
#' and the Cross-Validated Prognosis Observed Loss (CVPOL), two estimators of
#' EPOCE. The latter corrects the MPOL estimate for over-optimism by
#' approximated cross-validation. On external, this function only computes
#' MPOL.
#'
#'
#' @usage epoce(fit, pred.times, newdata = NULL, newdata.Longi = NULL)
#' @param fit A jointPenal, longiPenal, trivPenal or trivPenalNL object.
#' @param pred.times Time or vector of times to compute epoce.
#' @param newdata Optional. In case of joint models obtained with
#' \code{frailtyPenal}, \code{trivPenal} or \code{trivPenalNL}. For models
#' inheriting from \code{trivPenal} or \code{trivPenalNL} class, if
#' \code{newdata} is given, \code{newdata.Longi} must be given as well. When
#' missing, the data used for estimating the fit are used, and CVPOL and MPOL
#' are computed (internal validation). When \code{newdata} is specified, only
#' MPOL is computed on this new dataset (external validation). The new dataset
#' and the dataset used in the estimation must have the same covariates with
#' the same coding without missing data.
#' @param newdata.Longi Optional. In case of joint models obtained with
#' \code{longiPenal}, \code{trivPenal} or \code{trivPenalNL}. For models
#' inheriting from \code{longiPenal}, if the \code{newdata.Longi} is given,
#' \code{newdata} must be \code{NULL}, but for models from \code{trivPenal} or
#' \code{trivPenalNL} class, if \code{newdata.Longi} is given, \code{newdata}
#' must be provided as well. The two datasets newdata and newdata.Longi must
#' include the information concerning the same patients with the same
#' characteristics and the appropriate data on follow up (recurrences for
#' \code{newdata} and longitudinal measurements for \code{newdata.Longi}).
#' @return \item{data}{name of the data used to compute epoce}
#' \item{new.data}{a boolean which is FALSE if computation is done on the same
#' data as for estimation, and TRUE otherwise} \item{pred.times}{time or vector
#' of times used in the function} \item{mpol}{values of MPOL for each
#' pred.times} \item{cvpol}{values of CVPOL for each pred.times}
#' \item{IndivContrib}{all the contributions to the log-likelihood for each
#' pred.times} \item{AtRisk}{number of subject still at risk for each
#' pred.times}
#' @references D. Commenges, B. Liquet, C. Proust-Lima (2012). Choice of
#' prognostic estimators in joint models by estimating differences of expected
#' conditional Kullback-Leibler risks. \emph{Biometrics}, \bold{68(2)},
#' 380-387.
#' @keywords misc
#' @export
#' @examples
#'
#'
#' \dontrun{
#'
#' ########################################
#' #### EPOCE on a joint frailty model ####
#' ########################################
#'
#' data(readmission)
#'
#' modJoint.gap <- frailtyPenal(Surv(t.start,t.stop,event)~ cluster(id) +
#' dukes + charlson + sex + chemo + terminal(death),
#' formula.terminalEvent = ~ dukes + charlson + sex + chemo ,
#' data = readmission, n.knots = 8, kappa =c(2.11e+08,9.53e+11),
#' recurrentAG=TRUE)
#'
#' # computation on the same dataset
#' temps <- c(200,500,800,1100)
#' epoce <- epoce(modJoint.gap,temps)
#'
#' print(epoce)
#' plot(epoce,type = "cvpol")
#'
#' # computation on a new dataset
#' # here a sample of readmission with the first 50 subjects
#' s <- readmission[1:100,]
#' epoce <- epoce(modJoint.gap,temps,newdata=s)
#'
#' print(epoce)
#' plot(epoce,type = "mpol")
#'
#' #################################################
#' #### EPOCE on a joint model for a biomarker ####
#' ######### and a terminal event ###############
#' #################################################
#'
#' data(colorectal)
#' data(colorectalLongi)
#'
#' # Survival data preparation - only terminal events
#' colorectalSurv <- subset(colorectal, new.lesions == 0)
#'
#' modLongi <- longiPenal(Surv(time0, time1, state) ~ age +
#' treatment + who.PS, tumor.size ~ year*treatment + age +
#' who.PS, colorectalSurv, data.Longi =colorectalLongi,
#' random = c("1", "year"), id = "id", link = "Random-effects",
#' left.censoring = -3.33, hazard = "Weibull",
#' method.GH = "Pseudo-adaptive")
#'
#' # computation on the same dataset
#' time <- c(1, 1.5, 2, 2.5)
#' epoce <- epoce(modLongi,time)
#'
#' print(epoce)
#' plot(epoce, type = "cvpol")
#'
#' # computation on a new dataset
#' # here a sample of colorectal data with the first 50 subjects
#' s <- subset(colorectal, new.lesions == 0 & id%in%1:50)
#' s.Longi <- subset(colorectalLongi, id%in%1:50)
#' epoce <- epoce(modLongi, time, newdata = s, newdata.Longi = s.Longi)
#'
#' print(epoce)
#' plot(epoce, type = "mpol")
#'
#'
#' ###################################################
#' #### EPOCE on a joint model for a biomarker, ######
#' #### recurrent events and a terminal event ######
#' ###################################################
#'
#' data(colorectal)
#' data(colorectalLongi)
#'
#' # Linear model for the biomarker
#' # (computation takes around 30 minutes)
#' model.trivPenalNL <-trivPenal(Surv(gap.time, new.lesions) ~ cluster(id)
#' + age + treatment + who.PS + prev.resection + terminal(state),
#' formula.terminalEvent =~ age + treatment + who.PS + prev.resection,
#' tumor.size ~ year * treatment + age + who.PS, data = colorectal,
#' data.Longi = colorectalLongi, random = c("1", "year"), id = "id",
#' link = "Random-effects", left.censoring = -3.33, recurrentAG = FALSE,
#' hazard = "Weibull", method.GH="Pseudo-adaptive", n.nodes=7)
#'
#' # computation on the same dataset
#' time <- c(1, 1.5, 2, 2.5)
#'
#' # (computation takes around 10 minutes)
#' epoce <- epoce(model.trivPenalNL,time)
#' print(epoce)
#' plot(epoce, type = "cvpol")
#'
#' # computation on a new dataset
#' # here a sample of colorectal data with the first 100 subjects
#' s <- subset(colorectal, id%in%1:100)
#' s.Longi <- subset(colorectalLongi, id%in%1:100)
#' # (computation takes around 10 minutes)
#' epoce <- epoce(model.trivPenalNL, time, newdata = s, newdata.Longi = s.Longi)
#'
#' print(epoce)
#' plot(epoce, type = "mpol")
#'
#'
#'
#' # Non-linear model for the biomarker
#'
#' # No information on dose - creation of a dummy variable
#' colorectalLongi$dose <- 1
#'
#' # (computation can take around 40 minutes)
#' model.trivPenalNL <- trivPenalNL(Surv(time0, time1, new.lesions) ~ cluster(id) + age + treatment
#' + terminal(state), formula.terminalEvent =~ age + treatment, biomarker = "tumor.size",
#' formula.KG ~ 1, formula.KD ~ treatment, dose = "dose", time.biomarker = "year",
#' data = colorectal, data.Longi =colorectalLongi, random = c("y0", "KG"), id = "id",
#' init.B = c(-0.22, -0.16, -0.35, -0.19, 0.04, -0.41, 0.23), init.Alpha = 1.86,
#' init.Eta = c(0.5, 0.57, 0.5, 2.34), init.Biomarker = c(1.24, 0.81, 1.07, -1.53),
#' recurrentAG = TRUE, n.knots = 5, kappa = c(0.01, 2), method.GH = "Pseudo-adaptive")
#'
#' # computation on the same dataset
#' time <- c(1, 1.5, 2, 2.5)
#'
#' epoce <- epoce(model.trivPenalNL, time)
#'
#'
#' }
#'
#'
epoce <- function(fit, pred.times, newdata = NULL, newdata.Longi = NULL){
if (missing(fit)) stop("The argument fit must be specified")
if (!inherits(fit, "jointPenal") & !inherits(fit, "longiPenal") & !inherits(fit, "trivPenal") & !inherits(fit, "trivPenalNL")) stop("The argument fit must be a class 'jointPenal', 'longiPenal' or 'trivPenal'")
if (missing(pred.times)) stop("The argument pred.times must be specified")
if (!inherits(pred.times, "numeric")) stop("pred.times must contain numerical values")
if(!missing(newdata) & (!inherits(newdata, "data.frame"))) stop("The argument newdata must be a 'data.frame'")
if(!missing(newdata.Longi) & (!inherits(newdata.Longi, "data.frame"))) stop("The argument newdata must be a 'data.frame'")
if(inherits(fit, "jointPenal") & !missing(newdata.Longi))warning("The argument newdata.Longi is not required and thus ignored")
if(inherits(fit, "longiPenal") & !missing(newdata.Longi) & missing(newdata))warning("For an object of class 'longiPenal' both datasets should be given")
if(inherits(fit, "trivPenal") & !missing(newdata.Longi) & missing(newdata))warning("For an object of class 'trivPenal' both datasets should be given")
if(inherits(fit, "trivPenalNL") & missing(newdata.Longi) & !missing(newdata))warning("For an object of class 'trivPenalNL' both datasets should be given")
nt <- length(pred.times)
vopt <- fit$varHtotal
b <- fit$b
np <- length(fit$b)
typeof <- fit$typeof
nva <- fit$nvar
if (typeof == 0){
nz <- fit$n.knots
zi <- fit$zi
}else{
nz <- 0
zi <- 0
}
if (typeof == 1){
nbintervR <- fit$nbintervR
nbintervDC <- fit$nbintervDC
ttt <- fit$time
tttdc <- fit$timedc
}else{
nbintervR <- 0
nbintervDC <- 0
ttt <- 0
tttdc <- 0
}
# recuperation des profils d'individus
m <- fit$call
m0 <- match.call()
if (!missing(newdata)){
if (length(colnames(eval(m$data)))!=length(colnames(eval(m0$newdata)))) stop("Your new dataset must have the same number of columns than the dataset used in the 'fit'")
if (any(colnames(eval(m$data))!=colnames(eval(m0$newdata)))) stop("Your new dataset must have the very same variables than the dataset used in the 'fit'")
}
if (!missing(newdata.Longi)){
if (length(colnames(eval(m$data.Longi)))!=length(colnames(eval(m0$newdata.Longi)))) stop("Your new dataset for longitudinal data must have the same number of columns than the dataset used in the 'fit'")
if (any(colnames(eval(m$data.Longi))!=colnames(eval(m0$newdata.Longi)))) stop("Your new dataset for longitudinal data must have the very same variables than the dataset used in the 'fit'")
}
if (is.null(m$recurrentAG)) recurrentAG <- FALSE
else recurrentAG <- TRUE
if(inherits(fit, "jointPenal") | inherits(fit, "trivPenal") | inherits(fit, "trivPenalNL")){
#m$formula.terminalEvent <- m$n.knots <- m$recurrentAG <- m$cross.validation <- m$kappa <- m$maxit <- m$hazard <- m$nb.int1 <-m$nb.int2 <- m$RandDist <- m$betaorder <- m$betaknots <- m$init.B <- m$LIMparam <- m$LIMlogl <- m$LIMderiv <- m$print.times <- m$init.Theta <- m$init.Alpha <- m$Alpha <- m$... <- NULL
m$formula.LongitudinalData <- m$formula.terminalEvent <- m$recurrentAG <- m$data.Longi <- m$n.knots <- m$random <- m$link <- m$id <- m$kappa <- m$maxit <- m$hazard <- m$nb.int <- m$betaorder <- m$betaknots <- m$init.B <- m$LIMparam <- m$LIMlogl <- m$LIMderiv <- m$left.censoring <- m$print.times <- m$init.Random <- m$init.Eta <- m$init.Alpha <- m$method.GH <- m$intercept <- m$n.nodes <- m$biomarker <- m$formula.KG <- m$formula.KD <- m$dose <- m$time.biomarker <- m$BoxCox <- m$init.Biomarker <- m$... <- m$RandDist <- NULL
m[[1]] <- as.name("model.frame")
if (!missing(newdata)) m[[3]] <- as.name(m0$newdata) # nouveau dataset
dataset <- eval(m, sys.parent())
typeofY <- attr(model.extract(dataset, "response"),"type")
Y <- model.extract(dataset, "response")
if (typeofY=="right"){
tt0 <- rep(0,dim(dataset)[1])
tt1 <- Y[,1]
c <- Y[,2]
}else{
tt0 <- Y[,1]
tt1 <- Y[,2]
c <- Y[,3]
}
tt0 <- as.numeric(tt0)
tt1 <- as.numeric(tt1)
c <- as.numeric(c)
class(m$formula) <- "formula"
special <- c("strata", "cluster", "subcluster", "terminal", "num.id", "timedep")
Terms <- terms(m$formula, special)#, data = m$data)
m$formula <- Terms
dropx <- NULL
tempc <- untangle.specials(Terms, "cluster", 1:10)
cluster <- strata(dataset[, tempc$vars], shortlabel = TRUE)
# numbers <- table(cluster)[order(unique(cluster))]
# newCluster <- rep(1:nsujet,numbers)
dropx <- c(dropx,tempc$terms)
tempterm <- untangle.specials(Terms, "terminal", 1:10)
terminal <- strata(dataset[, tempterm$vars], shortlabel = TRUE)
terminal <- as.numeric(as.character(terminal))
dropx <- c(dropx,tempterm$terms)
if (!is.null(dropx)) newTerms <- Terms[-dropx]
else newTerms <- Terms
X <- model.matrix(newTerms, dataset)
if (ncol(X) > 1) X <- X[, -1, drop = FALSE]
nva1 <- ncol(X)
if (!missing(newdata)){
nobs <- nrow(newdata)
nsujet <- length(unique(cluster))
}else{
nobs <- fit$n
nsujet <- fit$groups
}
if (!recurrentAG){
tt0dc <- aggregate(tt1,by=list(cluster),FUN=sum)[,2]#rep(0,nsujet)
tt1dc <- aggregate(tt1,by=list(cluster),FUN=sum)[,2]
}else{
tt0dc <- rep(0,nsujet)
tt1dc <- aggregate(tt1,by=list(cluster),FUN=function(x) x[length(x)])[,2]
}
cdc <- aggregate(terminal,by=list(cluster),FUN=function(x) x[length(x)])[,2]
m2 <- fit$call
m2$formula.LongitudinalData <- m2$n.knots <- m2$recurrentAG <- m2$cross.validation <- m2$kappa <- m2$maxit <- m2$hazard <- m2$nb.int1 <-m2$nb.int2 <- m2$RandDist <- m2$betaorder <- m2$betaknots <- m2$init.B <- m2$LIMparam <- m2$LIMlogl <- m2$LIMderiv <- m2$print.times <- m2$init.Theta <- m2$init.Alpha <- m2$Alpha <- m2$method.GH <- m2$intercept <- m2$init.Eta <- m2$data.Longi <- m2$init.Random <- m2$left.censoring <- m2$random <- m2$link <- m2$id <- m2$n.nodes <- m2$biomarker <- m2$formula.KG <- m2$formula.KD <- m2$dose <- m2$time.biomarker <- m2$BoxCox <- m2$init.Biomarker <- m2$... <- NULL
m2$formula[[3]] <- m2$formula.terminalEvent[[2]]
m2$formula.terminalEvent <- NULL
m2[[1]] <- as.name("model.frame")
if (!missing(newdata)) m2[[3]] <- as.name(m0$newdata) # nouveau dataset
datasetdc <- eval(m2, sys.parent())
class(m2$formula) <- "formula"
special2 <- c("strata", "timedep")
Terms2 <- terms(m2$formula, special2)#, data = m3$data)
X2 <- model.matrix(Terms2, datasetdc)
if (ncol(X2) > 1) X2 <- X2[, -1, drop = FALSE]
nva2 <- ncol(X2)
if (!is.null(ncol(X2))){
Xdc <- aggregate(X2[,1],by=list(cluster), FUN=function(x) x[length(x)])[,2]
if (ncol(X2)>1){
for (i in 2:ncol(X2)){
Xdc.i <- aggregate(X2[,i],by=list(cluster), FUN=function(x) x[length(x)])[,2]
Xdc <- cbind(Xdc,Xdc.i)
}
}
}else{
Xdc <- aggregate(X2,by=list(cluster), FUN=function(x) x[length(x)])[,2]
}
if (!missing(newdata) & length(fit$coef[1:(fit$nvarEnd+fit$nvarRec)])!=(ncol(X)+ncol(X2))) stop("Different covariates in model and newdata. Verify your dataset, be careful to the factor variables.")
}
if(inherits(fit, "trivPenal") | inherits(fit, "longiPenal")){
m2 <- fit$call
m2$formula <- m2$formula.terminalEvent <- m2$data <- m2$random <- m2$id <- m2$link <- m2$n.knots <- m2$kappa <- m2$maxit <- m2$hazard <- m2$nb.int <- m2$betaorder <- m2$betaknots <- m2$init.B <- m2$LIMparam <- m2$LIMlogl <- m2$LIMderiv <- m2$print.times <- m2$left.censoring <- m2$init.Random <- m2$init.Eta <- m2$method.GH <- m2$intercept <- m2$... <- NULL
if (!missing(newdata.Longi)){m2[[3]] <- as.name(m0$newdata.Longi) # nouveau dataset
data.Longi <- newdata.Longi
}else{data.Longi <- eval(m2$data.Longi) }
special <- c("strata", "cluster", "subcluster", "terminal","num.id","timedep")
#========= Longitudinal Data preparation =========================
class(m2$formula.LongitudinalData) <- "formula"
TermsY <- terms(m2$formula.LongitudinalData, special, data = data.Longi)
llY <- attr(TermsY, "term.labels")#liste des variables explicatives
ord <- attr(TermsY, "order")
#=========================================================>
name.Y <- as.character(attr(TermsY, "variables")[[2]])
yy <- data.Longi[,which(names(data.Longi)==name.Y)]
# on identifie les variables explicatives facteurs avec nombre de niveau plus que 2
ind.placeY <- which(llY%in%names(which(lapply(data.Longi[,which(names(data.Longi)%in%llY)],function(x) length(levels(x)))>2)))
vec.factorY <- NULL
vec.factorY <- c(vec.factorY,llY[ind.placeY])
mat.factorY <- matrix(vec.factorY,ncol=1,nrow=length(vec.factorY))
# Fonction servant a prendre les termes entre "as.factor"
vec.factorY <-apply(mat.factorY,MARGIN=1,FUN=function(x){
if (length(grep("as.factor",x))>0){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))
}else{
return(x)
}})
ind.placeY <- grep(paste(vec.factorY,collapse="|"),llY)
if(is.factor(data.Longi[,names(data.Longi)==llY[1]]))X_L<- as.numeric(data.Longi[,names(data.Longi)==llY[1]])-1
else X_L<- data.Longi[,names(data.Longi)==llY[1]]
if(length(llY)>1){
for(i in 2:length(llY)){
if(is.factor(data.Longi[,names(data.Longi)==llY[i]]))X_L<- cbind(X_L,as.numeric(data.Longi[,names(data.Longi)==llY[i]])-1)
else X_L<- cbind(X_L,data.Longi[,names(data.Longi)==llY[i]])
}}
#X_L<- data.Longi[,names(data.Longi)%in%(llY)]
if(sum(ord)>length(ord)){
for(i in 1:length(ord)){
if(ord[i]>1){
v1 <- strsplit(as.character(llY[i]),":")[[1]][1]
v2 <- strsplit(as.character(llY[i]),":")[[1]][2]
if(is.factor(data.Longi[,names(data.Longi)==v1]) && length(levels(data.Longi[,names(data.Longi)==v1]))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
if(is.factor(data.Longi[,names(data.Longi)==v2]) && length(levels(data.Longi[,names(data.Longi)==v2]))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
if(is.factor(data.Longi[,names(data.Longi)==v1]) || !is.factor(data.Longi[,names(data.Longi)==v2])){
X_L <- cbind(X_L,(as.numeric(data.Longi[,names(data.Longi)==v1])-1)*data.Longi[,names(data.Longi)==v2])
llY[i]<-paste(llY[i],levels(data.Longi[,names(data.Longi)==v1])[2],sep="")
}else if(!is.factor(data.Longi[,names(data.Longi)==v1]) || is.factor(data.Longi[,names(data.Longi)==v2])){
X_L <- cbind(X_L,data.Longi[,names(data.Longi)==v1]*(as.numeric(data.Longi[,names(data.Longi)==v2])-1))
llY[i]<-paste(llY[i],levels(data.Longi[,names(data.Longi)==v2])[2],sep="")
}else{
X_L <- cbind(X_L,data.Longi[,names(data.Longi)==v1]*data.Longi[,names(data.Longi)==v2])
}
}
}
}
if(dim(X_L)[2]!=length(llY))stop("The variables in the longitudinal part must be in the data.Longi")
X_L <- as.data.frame(X_L)
names(X_L) <- llY
Intercept <- rep(1,dim(X_L)[1])
if(fit$intercept)X_L <- cbind(Intercept,X_L)
X_Lall<- X_L
"%+%"<- function(x,y) paste(x,y,sep="")
if(length(vec.factorY) > 0){
for(i in 1:length(vec.factorY)){
X_L <- cbind(X_L[,-(which(names(X_L)==vec.factorY[i]))],model.matrix(as.formula("~"%+%0%+%"+"%+%paste(vec.factorY[i], collapse= "+")), data.Longi)[,-1])
}
vect.factY<-names(X_L)[which(!(names(X_L)%in%llY))]
occurY <- rep(0,length(vec.factorY))
for(i in 1:length(vec.factorY)){
#occur[i] <- sum(vec.factor[i] == vect.fact)
occurY[i] <- length(grep(vec.factorY[i],vect.factY))
}
}
if (ncol(X_L) == 0){
noVarY <- 1
}else{
noVarY <- 0
}
#=========================================================>
clusterY <- data.Longi$id
maxy_rep <- max(table(clusterY))
uni.cluster<-as.factor(unique(clusterY))
npred <- length(uni.cluster)
nva3<-ncol(X_L)
varY <- as.matrix(sapply(X_L, as.numeric))
#=======================================>
#======= Construction du vecteur des indicatrice
if(length(vec.factorY) > 0){
# ind.place <- ind.place -1
k <- 0
for(i in 1:length(vec.factorY)){
ind.placeY[i] <- ind.placeY[i]+k
k <- k + occurY[i]-1
}
}
if(fit$link=="Random-effects")link <- 1
if(fit$link=="Current-level") link <- 2
matzy <- NULL
names.matzy <- fit$names.re
matzy <- data.matrix(X_Lall[,which(names(X_Lall)%in%names.matzy)])
if(fit$leftCensoring==FALSE){s_cag_id = 0
s_cag = 0}else{
s_cag_id = 1
s_cag = fit$leftCensoring.threshold
}
}
if(inherits(fit, "longiPenal")){
#m$formula.terminalEvent <- m$n.knots <- m$recurrentAG <- m$cross.validation <- m$kappa <- m$maxit <- m$hazard <- m$nb.int1 <-m$nb.int2 <- m$RandDist <- m$betaorder <- m$betaknots <- m$init.B <- m$LIMparam <- m$LIMlogl <- m$LIMderiv <- m$print.times <- m$init.Theta <- m$init.Alpha <- m$Alpha <- m$... <- NULL
m$formula.LongitudinalData <- m$formula.terminalEvent <- m$recurrentAG <- m$data.Longi <- m$n.knots <- m$random <- m$link <- m$id <- m$kappa <- m$maxit <- m$hazard <- m$nb.int <- m$betaorder <- m$betaknots <- m$init.B <- m$LIMparam <- m$LIMlogl <- m$LIMderiv <- m$left.censoring <- m$print.times <- m$init.Random <- m$init.Eta <- m$init.Alpha <- m$method.GH <- m$intercept <- m$n.nodes <- m$... <- NULL
m[[1]] <- as.name("model.frame")
if (!missing(newdata)) m[[3]] <- as.name(m0$newdata) # nouveau dataset
dataset <- eval(m, sys.parent())
typeofY <- attr(model.extract(dataset, "response"),"type")
Y <- model.extract(dataset, "response")
if (typeofY=="right"){
tt0dc <- rep(0,dim(dataset)[1])
tt1dc <- Y[,1]
cdc <- Y[,2]
} else {
tt0dc <- Y[,1]
tt1dc <- Y[,2]
cdc <- Y[,3]
}
tt0dc <- as.numeric(tt0dc)
tt1dc <- as.numeric(tt1dc)
cdc <- as.numeric(cdc)
class(m$formula) <- "formula"
special <- c("strata", "cluster", "subcluster", "num.id", "timedep")
Terms <- terms(m$formula, special)#, data = m$data)
m$formula <- Terms
dropx <- NULL
newTerms <- Terms
Xdc <- model.matrix(newTerms, dataset)
if (ncol(Xdc) > 1) Xdc <- Xdc[, -1, drop = FALSE]
nva2 <- ncol(Xdc)
if (!missing(newdata)){
nsujet <- dim(newdata)[1]
}else{
# nobs <- fit$n
nsujet <- fit$groups
}
}
if(inherits(fit, "trivPenalNL")){
nva3 <- fit$nvarKG
nva4 <- fit$nvarKD
m3 <- fit$call # longitudinal (KG)
m3$formula <- m3$formula.terminalEvent <- m3$biomarker <- m3$formula.KD <- m3$dose <- m3$data <- m3$recurrentAG <- m3$random <- m3$id <- m3$link <- m3$n.knots <- m3$kappa <- m3$maxit <- m3$hazard <- m3$init.B <- m3$LIMparam <- m3$LIMlogl <- m3$LIMderiv <- m3$print.times <- m3$left.censoring <- m3$init.Random <- m3$init.Eta <- m3$init.Alpha <- m3$method.GH <- m3$n.nodes <- m3$init.GH <- m3$time.biomarker <- m3$BoxCox <- m3$... <- NULL
Names.data.Longi <- m3$data.Longi
if (!missing(newdata.Longi)){m3[[3]] <- as.name(m0$newdata.Longi) # nouveau dataset
data.Longi <- newdata.Longi
}else{data.Longi <- eval(m3$data.Longi) }
formula.KG <- fit$formula.KG
m4 <- fit$call # longitudinal (KD)
m4$formula <- m4$formula.terminalEvent <- m4$biomarker <- m4$formula.KG <- m4$dose <- m4$data <- m4$recurrentAG <- m4$random <- m4$id <- m4$link <- m4$n.knots <- m4$kappa <- m4$maxit <- m4$hazard <- m4$init.B <- m4$LIMparam <- m4$LIMlogl <- m4$LIMderiv <- m4$print.times <- m4$left.censoring <- m4$init.Random <- m4$init.Eta <- m4$init.Alpha <- m4$method.GH <- m4$n.nodes <- m4$init.GH <- m4$time.biomarker <- m4$BoxCox <- m4$... <- NULL
Y <- data.Longi[,which(names(data.Longi)==fit$biomarker)]
if(!is.null(formula.KG[3]) && formula.KG[3] != "1()"){
TermsKG <- if (missing(data.Longi)){
terms(formula.KG, special)
}else{
terms(formula.KG, special, data = data.Longi)
}
ord <- attr(TermsKG, "order") # longueur de ord=nbre de var.expli
#si pas vide tous si il ya au moins un qui vaut 1 on arrete
m2$formula.KG <- TermsKG
m2[[1]] <- as.name("model.frame") # m[[1]]=frailtypenal, il le remplace par model.frame en fait
if (NROW(m3) == 0)stop("No (non-missing) observations") #nombre ligne different de 0
llKG <- attr(TermsKG, "term.labels")#liste des variables explicatives
#=========================================================>
name.KG <- as.character(attr(TermsKG, "variables")[[2]])
KG <- data.Longi[,which(names(data.Longi)==name.KG)]
# on identifie les variables explicatives facteurs avec nombre de niveau plus que 2
ind.placeKG <- which(lapply(as.data.frame(data.Longi[,which(names(data.Longi)%in%llKG)]),function(x) length(levels(x)))>2)
defined.factor <- llKG[grep("factor",llKG)]
vec.factorKG.tmp <- NULL
if(length(defined.factor)>0){
mat.factorKG.tmp <- matrix(defined.factor,ncol=1,nrow=length(defined.factor))
# Fonction servant a prendre les termes entre "as.factor"
vec.factorKG.tmp <-apply(mat.factorKG.tmp,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos1,stop=pos2))])))>2)return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
else return(NaN)
}else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos3,stop=pos4))])))>2)return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
else return(NaN)
}else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos1,stop=pos2))])))>2 || length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos3,stop=pos4))])))>2)return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
else return(NaN)
}
}else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos1,stop=pos2))])))>2)return(substr(x,start=pos1,stop=pos2))
else return(NaN)
}
}else{
return(x)
}})
vec.factorKG.tmp <- vec.factorKG.tmp[which(vec.factorKG.tmp!="NaN")]
if(length(vec.factorKG.tmp)>0){
for(i in 1:length(vec.factorKG.tmp)){
if(length(grep(":",vec.factorKG.tmp[i]))==0){
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==vec.factorKG.tmp[i])])))>2)ind.placeKG <- c(ind.placeKG,which(llKG%in%paste("as.factor(",vec.factorKG.tmp[i],")",sep="")))
}
}}
}
ind.placeKG <- sort(ind.placeKG)
mat.factorKG2 <- matrix(llKG,ncol=1,nrow=length(llKG))
# Fonction servant a prendre les termes entre "as.factor"
llKG2 <-apply(mat.factorKG2,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0 && length(grep(":",x))==0 && unlist(strsplit(x,split=""))[length(unlist(strsplit(x,split="")))]==")"){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
x<-substr(x,start=pos1,stop=pos2)
return(paste(x,levels(as.factor(data.Longi[,which(names(data.Longi)==x)]))[2],sep=""))
}else{
return(x)
}})
# Fonction servant a prendre les termes entre "as.factor" - without the name of the level
llKG3 <-apply(mat.factorKG2,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0 && length(grep(":",x))==0 && unlist(strsplit(x,split=""))[length(unlist(strsplit(x,split="")))]==")"){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))
}else{
return(x)
}})
llKG.real.names <- llKG3
llKG3 <- llKG3[!llKG2%in%llKG]
if(is.factor(data.Longi[,names(data.Longi)==llKG.real.names[1]]))X_L<- as.numeric(data.Longi[,names(data.Longi)==llKG.real.names[1]])-1
else X_L<- data.Longi[,names(data.Longi)==llKG.real.names[1]]
if(length(llKG) == 1)X_L <- as.data.frame(X_L)
if(length(llKG)>1){
for(i in 2:length(llKG.real.names)){
if(is.factor(data.Longi[,names(data.Longi)==llKG.real.names[i]]))X_L<- cbind(X_L,as.numeric(data.Longi[,names(data.Longi)==llKG.real.names[i]])-1)
else X_L<- cbind(X_L,data.Longi[,names(data.Longi)==llKG.real.names[i]])
}
}
#X_L<- data.Longi[,names(data.Longi)%in%(llY)]
llKG.fin <- llKG.real.names
llKG <- llKG.real.names
if(sum(ord)>length(ord)){
for(i in 1:length(ord)){
if(ord[i]>1){
name_v1 <- strsplit(as.character(llKG[i]),":")[[1]][1]
name_v2 <- strsplit(as.character(llKG[i]),":")[[1]][2]
if(length(grep("factor",name_v1))>0){name_v1<-substring(name_v1,11,nchar(name_v1)-1)
v1 <- as.factor(data.Longi[,names(data.Longi)==name_v1])}
else{v1 <- data.Longi[,names(data.Longi)==name_v1]}
if(length(grep("factor",name_v2))>0){name_v2<-substring(name_v2,11,nchar(name_v2)-1)
v2 <- as.factor(data.Longi[,names(data.Longi)==name_v2])}
else{v2 <- data.Longi[,names(data.Longi)==name_v2]}
llKG[i] <- paste(name_v1,":",name_v2,sep="")
# if(is.factor(v1) && length(levels(v1))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
# if(is.factor(v2) && length(levels(v2))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
if(is.factor(v1) && !is.factor(v2)){
dummy <- model.matrix( ~ v1 - 1)
# if(length(levels(v1)>2))vec.factorY <- c(vec.factorY,paste(name_v1,":",name_v2,sep=""))
for(j in 2:length(levels(v1))){
X_L <- cbind(X_L,dummy[,j]*v2)
if(i>1 && i<length(llKG.fin))llKG.fin <- c(llKG.fin[1:(i-1+j-2)],paste(name_v1,".",levels(v1)[j],":",name_v2,sep=""),llKG.fin[(i+1+j-2):length(llKG.fin)])
else if(i==length(llKG.fin))llKG.fin <- c(llKG.fin[1:(i-1+j-2)],paste(name_v1,".",levels(v1)[j],":",name_v2,sep=""))
else llKG.fin <- c(paste(name_v1,".",levels(v1)[j],":",name_v2,sep=""),llKG.fin[(2+j-2):length(llKG.fin)])
}
}else if(!is.factor(v1) && is.factor(v2)){
dummy <- model.matrix( ~ v2 - 1)
# if(length(levels(v2)>2))vec.factorY <- c(vec.factorY,paste(name_v1,":",name_v2,sep=""))
for(j in 2:length(levels(v2))){
X_L <- cbind(X_L,dummy[,j]*v1)
if(i>1 && i<length(llKG.fin))llKG.fin <- c(llKG.fin[1:(i-1+j-2)],paste(name_v1,":",name_v2,levels(v2)[j],sep=""),llKG.fin[(i+1+j-2):length(llKG.fin)])
else if(i==length(llKG.fin))llKG.fin <- c(llKG.fin[1:(i-1+j-2)],paste(name_v1,":",name_v2,levels(v2)[j],sep=""))
else llKG.fin <- c(paste(name_v1,":",name_v2,levels(v2)[j],sep=""),llKG.fin[(2+j-2):length(llKG.fin)])
}
}else if(is.factor(v1) && is.factor(v2)){
dummy1 <- model.matrix( ~ v1 - 1)
dummy2 <- model.matrix( ~ v2 - 1)
# if(length(levels(v1)>2) || length(levels(v2)>2))vec.factorY <- c(vec.factorY,paste(name_v1,":",name_v2,sep=""))
for(j in 2:length(levels(v1))){
for(k in 2:length(levels(v2))){
X_L <- cbind(X_L,dummy1[,j]*dummy2[,k])
if(i>1 && i<length(llKG.fin))llKG.fin <- c(llKG.fin[1:(i-1+j-2+k-2)],paste(name_v1,levels(v1)[j],":",name_v2,levels(v2)[k],sep=""),llKG.fin[(i+1+j-2+k-2):length(llKG.fin)])
else if(i==length(llKG.fin))llKG.fin <- c(llKG.fin[1:(i-1+j-2+k-2)],paste(name_v1,levels(v1)[j],":",name_v2,levels(v2)[k],sep=""))
else llKG.fin <- c(paste(name_v1,levels(v1)[j],":",name_v2,levels(v2)[k],sep=""),llKG.fin[(2+j-2+k-2):length(llKG.fin)])
}
}
}else{
X_L <- cbind(X_L,v1*v2)
}
}
}
}
if(length(grep(":",llKG))>0){
for(i in 1:length(grep(":",llKG))){
if(length(levels(data.Longi[,which(names(data.Longi)%in%strsplit(llKG[grep(":",llKG)[i]],":")[[1]])[1]]))>2 || length(levels(data.Longi[,which(names(data.Longi)%in%strsplit(llKG[grep(":",llKG)[i]],":")[[1]])[2]]))>2){
ind.placeKG <- c(ind.placeKG,grep(":",llKG)[i])
# vec.factorY <- c(vec.factorY,llY[grep(":",llY)[i]])
}
}
}
vec.factorKG <- NULL
if(length(vec.factorKG.tmp)>0)vec.factorKG <- c(llKG[ind.placeKG],vec.factorKG.tmp)
else vec.factorKG <- c(vec.factorKG,llKG[ind.placeKG])
vec.factorKG <- unique(vec.factorKG)
mat.factorKG <- matrix(vec.factorKG,ncol=1,nrow=length(vec.factorKG))
# Fonction servant a prendre les termes entre "as.factor" et (AK 04/11/2015) interactions
vec.factorKG <-apply(mat.factorKG,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
}
}else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))}
}else{
return(x)
}})
for(i in 1:length(llKG.fin)){
if(sum(names(data.Longi)==llKG.fin[i])>0){
if(is.factor(data.Longi[,names(data.Longi)==llKG.fin[i]]) && length(levels(data.Longi[,names(data.Longi)==llKG.fin[i]]))==2){
llKG.fin[i] <- paste(llKG.fin[i],levels(data.Longi[,names(data.Longi)==llKG.fin[i]])[2],sep="")}
}
}
# llY <- llY.fin
X_L <- as.data.frame(X_L)
if(dim(X_L)[2]!=length(llKG.fin))stop("The variables in the longitudinal part must be in the data.Longi")
names(X_L) <- llKG.fin
X_Lall<- X_L
"%+%"<- function(x,y) paste(x,y,sep="")
if(length(vec.factorKG) > 0){
for(i in 1:length(vec.factorKG)){
if(length(grep(":",vec.factorKG[i]))==0){
factor.spot <- which(names(X_L)==vec.factorKG[i])
if(factor.spot<ncol(X_L)) X_L <- cbind(X_L[1:(factor.spot-1)],model.matrix(as.formula("~"%+%0%+%"+"%+%paste(vec.factorKG[i], collapse= "+")), model.frame(~.,data.Longi,na.action=na.pass))[,-1],X_L[(factor.spot+1):ncol(X_L)])
else X_L <- cbind(X_L[1:(factor.spot-1)],model.matrix(as.formula("~"%+%0%+%"+"%+%paste(vec.factorKG[i], collapse= "+")), model.frame(~.,data.Longi,na.action=na.pass))[,-1])
} }
}
varKG <- as.matrix(sapply(X_L, as.numeric))
nsujety<-nrow(X_L)
}
#=========================================================>
clusterY <- data.Longi$id
max_rep <- max(table(clusterY))
uni.clusterY<-as.factor(unique(clusterY))
if(is.null(formula.KG[3]) | formula.KG[3] == "1()"){
varKG <- c()#rep(0, dim(data.Longi)[1])
nsujety <- length(Y)
}
formula.KD <- fit$formula.KD
if(!is.null(formula.KD[3]) && formula.KD[3] != "1()"){
TermsKD <- if (missing(data.Longi)){
terms(formula.KD, special)
}else{
terms(formula.KD, special, data = data.Longi)
}
ord <- attr(TermsKD, "order") # longueur de ord=nbre de var.expli
#si pas vide tous si il ya au moins un qui vaut 1 on arrete
m4$formula.KD <- TermsKD
m4[[1]] <- as.name("model.frame") # m[[1]]=frailtypenal, il le remplace par model.frame en fait
if (NROW(m4) == 0)stop("No (non-missing) observations") #nombre ligne different de 0
llKD <- attr(TermsKD, "term.labels")#liste des variables explicatives
#=========================================================>
name.KD <- as.character(attr(TermsKD, "variables")[[2]])
KD <- data.Longi[,which(names(data.Longi)==name.KD)]
# on identifie les variables explicatives facteurs avec nombre de niveau plus que 2
ind.placeKD <- which(lapply(as.data.frame(data.Longi[,which(names(data.Longi)%in%llKD)]),function(x) length(levels(x)))>2)
defined.factor <- llKD[grep("factor",llKD)]
vec.factorKD.tmp <- NULL
if(length(defined.factor)>0){
mat.factorKD.tmp <- matrix(defined.factor,ncol=1,nrow=length(defined.factor))
# Fonction servant a prendre les termes entre "as.factor"
vec.factorKD.tmp <-apply(mat.factorKG.tmp,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos1,stop=pos2))])))>2)return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
else return(NaN)
}else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos3,stop=pos4))])))>2)return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
else return(NaN)
}else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos1,stop=pos2))])))>2 || length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos3,stop=pos4))])))>2)return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
else return(NaN)
}
}else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==substr(x,start=pos1,stop=pos2))])))>2)return(substr(x,start=pos1,stop=pos2))
else return(NaN)
}
}else{
return(x)
}})
vec.factorKD.tmp <- vec.factorKD.tmp[which(vec.factorKD.tmp!="NaN")]
if(length(vec.factorKD.tmp)>0){
for(i in 1:length(vec.factorKD.tmp)){
if(length(grep(":",vec.factorKD.tmp[i]))==0){
if(length(levels(as.factor(data.Longi[,which(names(data.Longi)==vec.factorKD.tmp[i])])))>2)ind.placeKD <- c(ind.placeKD,which(llKD%in%paste("as.factor(",vec.factorKD.tmp[i],")",sep="")))
}
}}
}
ind.placeKD <- sort(ind.placeKD)
mat.factorKD2 <- matrix(llKD,ncol=1,nrow=length(llKD))
# Fonction servant a prendre les termes entre "as.factor"
llKD2 <-apply(mat.factorKD2,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0 && length(grep(":",x))==0 && unlist(strsplit(x,split=""))[length(unlist(strsplit(x,split="")))]==")"){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
x<-substr(x,start=pos1,stop=pos2)
return(paste(x,levels(as.factor(data.Longi[,which(names(data.Longi)==x)]))[2],sep=""))
}else{
return(x)
}})
# Fonction servant a prendre les termes entre "as.factor" - without the name of the level
llKD3 <-apply(mat.factorKD2,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0 && length(grep(":",x))==0 && unlist(strsplit(x,split=""))[length(unlist(strsplit(x,split="")))]==")"){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))
}else{
return(x)
}})
llKD.real.names <- llKD3
llKD3 <- llKD3[!llKD2%in%llKD]
if(is.factor(data.Longi[,names(data.Longi)==llKD.real.names[1]]))X_L<- as.numeric(data.Longi[,names(data.Longi)==llKD.real.names[1]])-1
else X_L<- data.Longi[,names(data.Longi)==llKD.real.names[1]]
if(length(llKD)>1){
for(i in 2:length(llKD.real.names)){
if(is.factor(data.Longi[,names(data.Longi)==llKD.real.names[i]]))X_L<- cbind(X_L,as.numeric(data.Longi[,names(data.Longi)==llKD.real.names[i]])-1)
else X_L<- cbind(X_L,data.Longi[,names(data.Longi)==llKD.real.names[i]])
}}
#X_L<- data.Longi[,names(data.Longi)%in%(llY)]
llKD.fin <- llKD.real.names
llKD <- llKD.real.names
if(sum(ord)>length(ord)){
for(i in 1:length(ord)){
if(ord[i]>1){
name_v1 <- strsplit(as.character(llKD[i]),":")[[1]][1]
name_v2 <- strsplit(as.character(llKD[i]),":")[[1]][2]
if(length(grep("factor",name_v1))>0){name_v1<-substring(name_v1,11,nchar(name_v1)-1)
v1 <- as.factor(data.Longi[,names(data.Longi)==name_v1])}
else{v1 <- data.Longi[,names(data.Longi)==name_v1]}
if(length(grep("factor",name_v2))>0){name_v2<-substring(name_v2,11,nchar(name_v2)-1)
v2 <- as.factor(data.Longi[,names(data.Longi)==name_v2])}
else{v2 <- data.Longi[,names(data.Longi)==name_v2]}
llKD[i] <- paste(name_v1,":",name_v2,sep="")
# if(is.factor(v1) && length(levels(v1))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
# if(is.factor(v2) && length(levels(v2))>2)stop("Interactions not allowed for factors with 3 or more levels (yet)")
if(is.factor(v1) && !is.factor(v2)){
dummy <- model.matrix( ~ v1 - 1)
# if(length(levels(v1)>2))vec.factorY <- c(vec.factorY,paste(name_v1,":",name_v2,sep=""))
for(j in 2:length(levels(v1))){
X_L <- cbind(X_L,dummy[,j]*v2)
if(i>1 && i<length(llKD.fin))llKD.fin <- c(llKD.fin[1:(i-1+j-2)],paste(name_v1,".",levels(v1)[j],":",name_v2,sep=""),llKD.fin[(i+1+j-2):length(llKD.fin)])
else if(i==length(llKD.fin))llKD.fin <- c(llKD.fin[1:(i-1+j-2)],paste(name_v1,".",levels(v1)[j],":",name_v2,sep=""))
else llKD.fin <- c(paste(name_v1,".",levels(v1)[j],":",name_v2,sep=""),llKD.fin[(2+j-2):length(llKD.fin)])
}
}else if(!is.factor(v1) && is.factor(v2)){
dummy <- model.matrix( ~ v2 - 1)
# if(length(levels(v2)>2))vec.factorY <- c(vec.factorY,paste(name_v1,":",name_v2,sep=""))
for(j in 2:length(levels(v2))){
X_L <- cbind(X_L,dummy[,j]*v1)
if(i>1 && i<length(llKD.fin))llKD.fin <- c(llKD.fin[1:(i-1+j-2)],paste(name_v1,":",name_v2,levels(v2)[j],sep=""),llKD.fin[(i+1+j-2):length(llKD.fin)])
else if(i==length(llKD.fin))llKD.fin <- c(llKD.fin[1:(i-1+j-2)],paste(name_v1,":",name_v2,levels(v2)[j],sep=""))
else llKD.fin <- c(paste(name_v1,":",name_v2,levels(v2)[j],sep=""),llKD.fin[(2+j-2):length(llKD.fin)])
}
}else if(is.factor(v1) && is.factor(v2)){
dummy1 <- model.matrix( ~ v1 - 1)
dummy2 <- model.matrix( ~ v2 - 1)
# if(length(levels(v1)>2) || length(levels(v2)>2))vec.factorY <- c(vec.factorY,paste(name_v1,":",name_v2,sep=""))
for(j in 2:length(levels(v1))){
for(k in 2:length(levels(v2))){
X_L <- cbind(X_L,dummy1[,j]*dummy2[,k])
if(i>1 && i<length(llKD.fin))llKD.fin <- c(llKD.fin[1:(i-1+j-2+k-2)],paste(name_v1,levels(v1)[j],":",name_v2,levels(v2)[k],sep=""),llKG.fin[(i+1+j-2+k-2):length(llKG.fin)])
else if(i==length(llKD.fin))llKD.fin <- c(llKD.fin[1:(i-1+j-2+k-2)],paste(name_v1,levels(v1)[j],":",name_v2,levels(v2)[k],sep=""))
else llKD.fin <- c(paste(name_v1,levels(v1)[j],":",name_v2,levels(v2)[k],sep=""),llKD.fin[(2+j-2+k-2):length(llKD.fin)])
}
}
}else{
X_L <- cbind(X_L,v1*v2)
}
}
}
}
if(length(grep(":",llKD))>0){
for(i in 1:length(grep(":",llKD))){
if(length(levels(data.Longi[,which(names(data.Longi)%in%strsplit(llKD[grep(":",llKD)[i]],":")[[1]])[1]]))>2 || length(levels(data.Longi[,which(names(data.Longi)%in%strsplit(llKD[grep(":",llKD)[i]],":")[[1]])[2]]))>2){
ind.placeKD <- c(ind.placeKD,grep(":",llKD)[i])
# vec.factorY <- c(vec.factorY,llY[grep(":",llY)[i]])
}
}
}
vec.factorKD <- NULL
if(length(vec.factorKD.tmp)>0)vec.factorKD <- c(llKD[ind.placeKD],vec.factorKD.tmp)
else vec.factorKD <- c(vec.factorKD,llKD[ind.placeKD])
vec.factorKD <- unique(vec.factorKD)
mat.factorKD <- matrix(vec.factorKD,ncol=1,nrow=length(vec.factorKD))
# Fonction servant a prendre les termes entre "as.factor" et (AK 04/11/2015) interactions
vec.factorKD <-apply(mat.factorKD,MARGIN=1,FUN=function(x){
if (length(grep("factor",x))>0){
if(length(grep(":",x))>0){
if(grep('\\(',unlist(strsplit(x,split="")))[1]<grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep(":",unlist(strsplit(x,split="")))[1]
pos4 <- length(unlist(strsplit(x,split="")))
return(paste(substr(x,start=pos1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else if(grep("\\(",unlist(strsplit(x,split="")))[1]>grep(":",unlist(strsplit(x,split="")))[1] && length(grep('\\(',unlist(strsplit(x,split=""))))==1){
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[1]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=1,stop=pos2),substr(x,start=pos3,stop=pos4),sep=""))
}else{#both factors
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- grep(":",unlist(strsplit(x,split="")))[1]-2
pos3 <- grep("\\(",unlist(strsplit(x,split="")))[2]+1
pos4 <- length(unlist(strsplit(x,split="")))-1
return(paste(substr(x,start=pos1,stop=pos2),":",substr(x,start=pos3,stop=pos4),sep=""))
}
}else{
pos1 <- grep("r",unlist(strsplit(x,split="")))[1]+2
pos2 <- length(unlist(strsplit(x,split="")))-1
return(substr(x,start=pos1,stop=pos2))}
}else{
return(x)
}})
for(i in 1:length(llKD.fin)){
if(sum(names(data.Longi)==llKD.fin[i])>0){
if(is.factor(data.Longi[,names(data.Longi)==llKD.fin[i]]) && length(levels(data.Longi[,names(data.Longi)==llKD.fin[i]]))==2){
llKD.fin[i] <- paste(llKD.fin[i],levels(data.Longi[,names(data.Longi)==llKD.fin[i]])[2],sep="")}
}
}
# llY <- llY.fin
X_L <- as.data.frame(X_L)
if(dim(X_L)[2]!=length(llKD.fin))stop("The variables in the longitudinal part must be in the data.Longi")
names(X_L) <- llKD.fin
X_Lall<- X_L
"%+%"<- function(x,y) paste(x,y,sep="")
if(length(vec.factorKD) > 0){
for(i in 1:length(vec.factorKD)){
if(length(grep(":",vec.factorKD[i]))==0){
factor.spot <- which(names(X_L)==vec.factorKD[i])
if(factor.spot<ncol(X_L)) X_L <- cbind(X_L[1:(factor.spot-1)],model.matrix(as.formula("~"%+%0%+%"+"%+%paste(vec.factorKD[i], collapse= "+")), model.frame(~.,data.Longi,na.action=na.pass))[,-1],X_L[(factor.spot+1):ncol(X_L)])
else X_L <- cbind(X_L[1:(factor.spot-1)],model.matrix(as.formula("~"%+%0%+%"+"%+%paste(vec.factorKD[i], collapse= "+")), model.frame(~.,data.Longi,na.action=na.pass))[,-1])
} }
}
varKD <- as.matrix(sapply(X_L, as.numeric))
}else{
varKD <- c()#rep(0, dim(data.Longi)[1])
}
matzy <- NULL # here matzy is for time and dose
matzy <- cbind(data.Longi[,which(colnames(data.Longi)==fit$time.biomarker)],data.Longi[,which(colnames(data.Longi)==fit$dose)])
if(dim(matzy)[2] != 2)stop("Both information on the biomarker measurement times and dose must be included in the data.")
matzy <- as.matrix(matzy)
if(fit$link=="Random-effects")link <- 1
if(fit$link=="Current-level") link <- 2
if(fit$leftCensoring==FALSE){
s_cag_id = 0
s_cag = 0
}else{
s_cag_id = 1
s_cag = fit$leftCensoring.threshold
}
}
cat("\n")
cat("Calculating ... \n")
if(inherits(fit, 'jointPenal')){
if(fit$logNormal==0){
ans <- .Fortran(C_cvpl,
as.integer(nobs),
as.integer(nsujet),
as.integer(cluster),
as.integer(c),
as.integer(cdc),
as.integer(nva1),
as.integer(nva2),
as.double(X),
as.double(Xdc),
as.integer(typeof),
as.integer(nz),
as.double(zi),
as.double(ttt),
as.double(tttdc),
as.integer(nbintervR),
as.integer(nbintervDC),
as.integer(np),
as.double(b),
as.double(vopt),
as.double(tt0),
as.double(tt1),
as.double(tt0dc),
as.double(tt1dc),
as.integer(nt),
as.double(pred.times),
rl_cond=as.double(rep(0,nt)),
epoir=as.double(rep(0,nt)),
contribt=as.double(rep(0,nt*nsujet)),
atrisk=as.double(rep(0,nt))
)
}else{
# cat('logn...')
ans <- .Fortran(C_cvpl_logn,
as.integer(nobs),
as.integer(nsujet),
as.integer(cluster),
as.integer(c),
as.integer(cdc),
as.integer(nva1),
as.integer(nva2),
as.double(X),
as.double(Xdc),
as.integer(typeof),
as.integer(nz),
as.double(zi),
as.double(ttt),
as.double(tttdc),
as.integer(nbintervR),
as.integer(nbintervDC),
as.integer(np),
as.double(b),
as.double(vopt),
as.double(tt0),
as.double(tt1),
as.double(tt0dc),
as.double(tt1dc),
as.integer(nt),
as.double(pred.times),
rl_cond=as.double(rep(0,nt)),
epoir=as.double(rep(0,nt)),
contribt=as.double(rep(0,nt*nsujet)),
atrisk=as.double(rep(0,nt))
)
}}else if(inherits(fit, "longiPenal")){
ans <- .Fortran(C_cvpl_long,
as.integer(nsujet),
as.integer(1),
as.integer(length(clusterY)),
as.integer(0),
as.integer(clusterY),
as.integer(0),
as.integer(cdc),
as.double(yy),
as.integer(1),
as.integer(nva2),
as.integer(nva3),
as.integer(fit$ne_re),
as.integer(0),
as.integer(fit$netadc),
as.integer(link),
as.double(matrix(0,nrow=1,ncol=1)),
as.double(Xdc),
as.double(as.matrix(varY)),
as.double(matzy),
as.double(s_cag),
as.integer(s_cag_id),
as.integer(typeof),
as.integer(2),
as.integer(nz),
as.double(zi),
as.integer(np),
as.double(b),
as.double(vopt),
as.double(0),
as.double(0),
as.double(tt0dc),
as.double(tt1dc),
as.integer(nt),
as.double(pred.times),
rl_cond=as.double(rep(0,nt)),
epoir=as.double(rep(0,nt)),
contribt=as.double(rep(0,nt*nsujet)),
atrisk=as.double(rep(0,nt))
)
}else if(inherits(fit, "trivPenal")){
ans <- .Fortran(C_cvpl_long,
as.integer(nsujet),
as.integer(nobs),
as.integer(length(clusterY)),
as.integer(cluster),
as.integer(clusterY),
as.integer(c),
as.integer(cdc),
as.double(yy),
as.integer(nva1),
as.integer(nva2),
as.integer(nva3),
as.integer(fit$ne_re),
as.integer(fit$netar),
as.integer(fit$netadc),
as.integer(link),
as.double(X),
as.double(Xdc),
as.double(as.matrix(varY)),
as.double(matzy),
as.double(s_cag),
as.integer(s_cag_id),
as.integer(typeof),
as.integer(3),
as.integer(nz),
as.double(zi),
as.integer(np),
as.double(b),
as.double(vopt),
as.double(tt0),
as.double(tt1),
as.double(tt0dc),
as.double(tt1dc),
as.integer(nt),
as.double(pred.times),
rl_cond=as.double(rep(0,nt)),
epoir=as.double(rep(0,nt)),
contribt=as.double(rep(0,nt*nsujet)),
atrisk=as.double(rep(0,nt))
)
}else if(inherits(fit, "trivPenalNL")){
GH <- c(0,0)
GH[1] <- ifelse(fit$methodGH == "Standard", 0, 1)
# GH[1] <- 0
GH[2] <- fit$n.nodes
box_cox <- c(0,1)
box_cox[1] <- ifelse(fit$BoxCox == TRUE, 1, 0)
if(!is.null(fit$BoxCox_parameter))box_cox[2] <- fit$BoxCox_parameter
lappend <- function (lst, ...){
lst <- c(lst, list(...))
return(lst)
}
nodes <- gauss.quad(20,kind="hermite")$nodes
weights <- gauss.quad(20,kind="hermite")$weights*exp(nodes^2)
nodes2 <- gauss.quad(20,kind="hermite")$nodes
weights2 <- gauss.quad(20,kind="hermite")$weights*exp(nodes2^2)
tmp <- rep(list(nodes), fit$ne_re)
tmp <- lappend(tmp, nodes2)
nodes <- as.data.frame(do.call(expand.grid,tmp))
tmp <- rep(list(weights), fit$ne_re)
tmp <- lappend(tmp, weights2)
weights <- as.data.frame(do.call(expand.grid,tmp))
nodes <- sapply(nodes, as.double)
weights <- sapply(weights, as.double)
ans <- .Fortran(C_cvplnl,
as.integer(nsujet),
as.integer(nobs),
as.integer(length(clusterY)),
as.integer(cluster),
as.integer(clusterY),
as.integer(c),
as.integer(cdc),
as.double(Y),
as.integer(nva1),
as.integer(nva2),
as.integer(nva3),
as.integer(nva4),
as.integer(fit$ne_re),
as.integer(fit$random.which),
as.double(box_cox),
as.integer(fit$netar),
as.integer(fit$netadc),
as.integer(link),
as.double(X),
as.double(Xdc),
as.double(cbind(varKG,varKD)),
as.double(matzy),
as.double(s_cag),
as.integer(s_cag_id),
as.integer(typeof),
as.integer(nz),
as.double(zi),
as.integer(np),
as.double(b),
as.double(vopt),
as.double(tt0),
as.double(tt1),
as.double(tt0dc),
as.double(tt1dc),
as.integer(nt),
as.double(pred.times),
rl_cond=as.double(rep(0,nt)),
epoir=as.double(rep(0,nt)),
contribt=as.double(rep(0,nt*nsujet)),
atrisk=as.double(rep(0,nt)),
as.integer(GH),
as.double(fit$b_pred),
as.double(fit$weights),#weights),
as.double(fit$nodes),#nodes),#
as.integer(fit$n.nodes^fit$ne_re*20)#
)
}
out <- NULL
if (!missing(newdata)) out$data <- m0$newdata
else out$data <- fit$data
out$new.data <- !is.null(newdata)
out$pred.times <- pred.times
out$mpol <- ans$rl_cond
if(!missing(newdata) && any(out$mpol<0))stop("The program stopped abnormally. This may be related to the new datasets with insufficient information")
if (missing(newdata)) out$cvpol <- ans$epoir
out$IndivContrib <- matrix(ans$contribt,nrow=nsujet,ncol=nt)
out$AtRisk <- ans$atrisk
cat("Estimators of EPOCE computed for",length(pred.times),"times \n")
class(out) <- c("epoce")
out
}
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.