user2739128
user2739128

Reputation: 1

mistake in multivePenal but not in frailtyPenal

The libraries used are: library(survival) library(splines) library(boot) library(frailtypack) and the function used is in the library frailty pack.
In my data I have two recurrent events(delta.stable and delta.unstable) and one terminal event (delta.censor). There are some time-varying explanatory variables, like unemployment rate(u.rate) (is quarterly) that's why my dataset has been splitted by quarters. Here there is a link to the subsample used in the code just below, just in case it may be helpful to see the mistake. https://www.dropbox.com/s/spfywobydr94bml/cr_05_males_services.rda The problem is that it takes a lot of time running until the warning message appear.

Main variables of the Survival function are: I have two recurrent events: delta.unstable (unst.): takes value one when the individual find an unstable job. delta.stable (stable): takes value one when the individual find a stable job. And one terminal event delta.censor (d.censor): takes value one when the individual has death, retired or emigrated.

row     id contadorbis unst. stable d.censor    .t0 .t
1   78  1   0   1   0   0   88
2   101 2   0   1   0   0   46
3   155 3   0   1   0   0   27
4   170 4   0   0   0   0   61
5   170 4   1   0   0   61  86
6   213 5   0   0   0   0   92
7   213 5   0   0   0   92  182
8   213 5   0   0   0   182 273
9   213 5   0   0   0   273 365
10  213 5   1   0   0   365 394
11  334 6   0   1   0   0   6
12  334 7   1   0   0   0   38
13  369 8   0   0   0   0   27
14  369 8   0   0   0   27  119
15  369 8   0   0   0   119 209
16  369 8   0   0   0   209 300
17  369 8   0   0   0   300 392

When I apply multivePenal I obtain the following message:

Error en aggregate.data.frame(as.data.frame(x), ...) : 
arguments must have same length
Además: Mensajes de aviso perdidos

In Surv(.t0, .t, delta.stable) : Stop time must be > start time, NA created
#### multivePenal function
fit.joint.05_malesP<multivePenal(Surv(.t0,.t,delta.stable)~cluster(contadorbis)+terminal(as.factor(delta.censor))+event2(delta.unstable),formula.terminalEvent=~1, formula2=~as.factor(h.skill),data=cr_05_males_serv,Frailty=TRUE,recurrentAG=TRUE,cross.validation=F,n.knots=c(7,7,7), kappa=c(1,1,1), maxit=1000, hazard="Splines")

I have checked if Surv(.t0,.t,delta.stable) contains NA, and there are no NA's.

In addition, when I apply for the same data the function frailtyPenal for both possible combinations, the function run well and I get results. I take one week looking at this and I do not find the key. I would appreciate some of light to this problem.

#delta unstable+death

enter code here

fit.joint.05_males<-frailtyPenal(Surv(.t0,.t,delta.unstable)~cluster(id)+u.rate+as.factor(h.skill)+as.factor(m.skill)+as.factor(non.manual)+as.factor(municipio)+as.factor(spanish.speakers)+ as.factor(no.spanish.speaker)+as.factor(Aged.16.19)+as.factor(Aged.20.24)+as.factor(Aged.25.29)+as.factor(Aged.30.34)+as.factor(Aged.35.39)+ as.factor(Aged.40.44)+as.factor(Aged.45.51)+as.factor(older61)+ as.factor(responsabilities)+
terminal(delta.censor),formula.terminalEvent=~u.rate+as.factor(h.skill)+as.factor(m.skill)+as.factor(municipio)+as.factor(spanish.speakers)+as.factor(no.spanish.speaker)+as.factor(Aged.16.19)+as.factor(Aged.20.24)+as.factor(Aged.25.29)+as.factor(Aged.30.34)+as.factor(Aged.35.39)+as.factor(Aged.40.44)+as.factor(Aged.45.51)+as.factor(older61)+ as.factor(responsabilities),data=cr_05_males_services,n.knots=12,kappa1=1000,kappa2=1000,maxit=1000, Frailty=TRUE,joint=TRUE, recurrentAG=TRUE)

###Be patient. The program is computing ... 
###The program took 2259.42 seconds 

#delta stable+death
fit.joint.05_males<frailtyPenal(Surv(.t0,.t,delta.stable)~cluster(id)+u.rate+as.factor(h.skill)+as.factor(m.skill)+as.factor(non.manual)+as.factor(municipio)+as.factor(spanish.speakers)+as.factor(no.spanish.speaker)+as.factor(Aged.16.19)+as.factor(Aged.20.24)+as.factor(Aged.25.29)+as.factor(Aged.30.34)+as.factor(Aged.35.39)+as.factor(Aged.40.44)+as.factor(Aged.45.51)+as.factor(older61)+as.factor(responsabilities)+terminal(delta.censor),formula.terminalEvent=~u.rate+as.factor(h.skill)+as.factor(m.skill)+as.factor(municipio)+as.factor(spanish.speakers)+as.factor(no.spanish.speaker)+as.factor(Aged.16.19)+as.factor(Aged.20.24)+as.factor(Aged.25.29)+as.factor(Aged.30.34)+as.factor(Aged.35.39)+as.factor(Aged.40.44)+as.factor(Aged.45.51)+as.factor(older61)+as.factor(responsabilities),data=cr_05_males_services,n.knots=12,kappa1=1000,kappa2=1000,maxit=1000, Frailty=TRUE,joint=TRUE, recurrentAG=TRUE)
###The program took 3167.15 seconds

Upvotes: 0

Views: 574

Answers (3)

Todd D
Todd D

Reputation: 278

I encountered the same error and arrived at this solution.

frailtyPenal() will not accept data.frames of different length. The data.frame used in Surv and data.frame named in data= in frailtyPenal must be the same length. I used a Cox regression to identify the incomplete cases, reset the survival object to exclude the missing cases and, finally, run frailtyPenal:

library(survival)
library(frailtypack)

data(readmission)

#Reproduce the error

#change the first start time to NA
  readmission[1,3] <- NA 

#create a survival object with one missing time
  surv.obj1 <- with(readmission, Surv(t.start, t.stop, event))

#observe the error
  frailtyPenal(surv.obj1 ~ cluster(id) + dukes, 
               data=readmission, 
               cross.validation=FALSE,
               n.knots=10,
               kappa=1,
               hazard="Splines")

#repair by resetting the surv object to omit the missing value(s)

#identify NAs using a Cox model
  cox.na <- coxph(surv.obj1 ~ dukes, data = readmission)

#remove the NA cases from the original set to create complete cases
  readmission2 <- readmission[-cox.na$na.action,]

#reset the survival object using the complete cases
  surv.obj2 <- with(readmission2, Surv(t.start, t.stop, event))

#run frailtyPenal using the complete cases dataset and the complete cases Surv object
  frailtyPenal(surv.obj2 ~ cluster(id) + dukes,
               data = readmission2, 
               cross.validation = FALSE,
               n.knots = 10,
               kappa = 1,
               hazard = "Splines")

Upvotes: 0

user2739128
user2739128

Reputation: 1

The authors of the package have told me that the function is not finished yet, so perhaps that is the reason that it is not working well.

Upvotes: 0

Henrik
Henrik

Reputation: 67818

Because you neither provide information about the packages used, nor the data necessary to run multivepenal or frailtyPenal, I can only help you with the Surv part (because I happened to have that package loaded).

The Surv warning message you provided (In Surv(.t0, .t, delta.stable) : Stop time must be > start time, NA created) suggests that something is strange with your variables .t0 (the time argument in Surv, refered to as 'start time' in the warning), and/or .t (time2 argument, 'Stop time' in the warning). I check this possibility with a simple example

# read the data you feed `Surv` with
df <- read.table(text = "row     id contadorbis unst. stable d.censor    .t0 .t
1   78  1   0   1   0   0   88
2   101 2   0   1   0   0   46
3   155 3   0   1   0   0   27
4   170 4   0   0   0   0   61
5   170 4   1   0   0   61  86
6   213 5   0   0   0   0   92
7   213 5   0   0   0   92  182
8   213 5   0   0   0   182 273
9   213 5   0   0   0   273 365
10  213 5   1   0   0   365 394
11  334 6   0   1   0   0   6
12  334 7   1   0   0   0   38
13  369 8   0   0   0   0   27
14  369 8   0   0   0   27  119
15  369 8   0   0   0   119 209
16  369 8   0   0   0   209 300
17  369 8   0   0   0   300 392", header = TRUE)

# create survival object
mysurv <- with(df, Surv(time = .t0, time2 = .t, event = stable))
mysurv

# create a new data set where one .t for some reason is less than .to
# on row five .t0 is 61, so I set .t to 60
df2 <- df
df2$.t[df2$.t == 86] <- 60

# create survival object using new data which contains at least one Stop time that is less than Start time
mysurv2 <- with(df2, Surv(time = .t0, time2 = .t, event = stable))

# Warning message:
#  In Surv(time = .t0, time2 = .t, event = stable) :
#  Stop time must be > start time, NA created
# i.e. the same warning message as you got

# check the survival object
mysurv2

# as you can see, the fifth interval contains NA
# I would recommend you check .t0 and .t in your data set carefully 
# one way to examine rows where Stop time (.t) is less than start time (.t0) is:
df2[which(df2$.t0 > df2$.t), ]

I am not familiar with multivepenal but it seems that it does not accept a survival object which contains intervals with NA, whereas might frailtyPenal might do so.

Upvotes: 1

Related Questions