-
Notifications
You must be signed in to change notification settings - Fork 12
Description
Dear Mark,
I have a question stemming from one of your examples in "rstpm2: a simple guide", where the following code is reported:
library(rstpm2)
brcancer <- transform(brcancer, recyear=rectime / 365.24)
fit <- stpm2(Surv(recyear,censrec==1)~hormon, data=brcancer, df=4)
predHormon <- predict(fit, newdata=data.frame(hormon=0:1),
type="hazard", grid=TRUE, full=TRUE, se.fit=TRUE)
nrow(predHormon)I see that predHormon is a data.frame object with 598 rows. It contains distinct predictions for the values Hormon = 0 and Hormon = 1, meaning that we actually have 299 rows for each class of Hormon.
This can be confirmed by selecting only one of the two classes as below:
predHormon0 <- predict(fit, newdata=data.frame(hormon=0),
type="hazard", grid=TRUE, full=TRUE, se.fit=TRUE)
nrow(predHormon0)My question is, what is the reason for having exactly 299 observations as the result of the prediction?
I have noticed that this does not seem to be data-driven, as I reached exactly 299 observations for class after using different datasets and different model specifications for the flexible model, as in the example below:
colon3<- data.frame(
sex = ifelse(rstpm2::colon$sex == "Female", 0, 1),
status = ifelse(rstpm2::colon$status %in% c("Dead: cancer", "Dead: other"), 1, ifelse(colon$status == "Alive", 0, NA)),
persontime = as.numeric(rstpm2::colon$exit - colon$dx))
modvc <- stpm2(Surv(persontime, status) ~ sex, df = 3, tvc = list(sex = 3),
data = colon3)
predColon0 <- predict(modvc, newdata=data.frame(sex = 0),
type="hazard", grid=TRUE, full=TRUE, se.fit=TRUE)
nrow(predColon0)I am then interested in getting a plot coming from the prediction (example below); however, the output is not completely satisfying up to now (rather than a smooth curve, a sum of small segments are visible), and I am guessing that being able to increase the number of points coming from the prediction might help in this sense
exampleplot <- plot(modvc,newdata=data.frame(sex=1), type="hr",
exposed=function(data) transform(data, sex=0),
xlab="Time",
rug = F,
lwd = 3)
nrow(exampleplot)Thank you again for your help and for this extremely helpful package!
Luigi