Demonstration of the ‘reslifefsr’ function
This function exclusively uses flexsurvreg() objects as the source of the data. We will show how to use the function below. First, we’ll start off with a simple example and show how the mean produced is equal to the integration.
We run a flexsurvreg() with a single covariate for the ‘gengamma.orig’ distribution. We then run the reslifefsr() function, using the flexsurvreg() object. We have a life of 1, which is a scalar value at which we calculate the residual life, and specify that we only want the mean.
library(flexsurv)
fs1 <- flexsurvreg(Surv(recyrs, censrec) ~ 1, data = bc,dist = "gengamma.orig")
r <- reslifefsr(obj = fs1, life = 1, type = 'mean')
r
#> [1] 7.194142We can verify the accuracy of our derivation using the integrate function.
b <- exp(as.numeric(fs1$coefficients[1]))
a <- exp(as.numeric(fs1$coefficients[2]))
k <- exp(as.numeric(fs1$coefficients[3]))
# numerical integration
f <- function(x) {
return(pgengamma.orig(x, shape = b, scale = a, k = k, lower.tail = FALSE))
}
mx_interal <- integrate(f, 1, Inf)$value/pgengamma.orig(1, shape = b, scale = a, k = k, lower.tail = FALSE)
mx_interal
#> [1] 7.194142As you can see, the mean values are the same.
Demonstration with multiple covariates
What if the flexsurvreg() object has multiple covariates? We will see this below. We will use the ‘bc’ dataset contained within ‘flexsurv’ with an added covariate ‘age’.
data(bc)
newbc <- bc
newbc$age <- rnorm(dim(bc)[1], mean = 65-scale(newbc$recyrs, scale=FALSE),sd = 5)Generating a flexsurvreg object:
We can now use the reslifefsr() function to generate residual life values. There are 686 values, since the length of the input is 686. We are only going to show 10, in order to save space.
head(reslifefsr(obj = fsr2, life= 4), 10)
#> [1] 10.991160 7.445017 11.888460 8.855456 7.320557 5.101370 7.224315
#> [8] 8.983472 4.498345 10.005028If we want to look at the median, we can specify that.
head(reslifefsr(obj = fsr2,life= 4, type = 'median'),10)
#> [1] 8.893677 5.931759 9.646690 7.106529 5.828371 3.995473 5.748458 7.213403
#> [9] 3.502050 8.067468We can look at all 3 types now (mean, median, and percentile).
head(reslifefsr(obj = fsr2, life= 4, p = .8, type = 'all'),10)
#> mean median percentile
#> 1 10.991160 8.893677 17.531507
#> 2 7.445017 5.931759 11.945311
#> 3 11.888460 9.646690 18.939333
#> 4 8.855456 7.106529 14.172096
#> 5 7.320557 5.828371 11.748439
#> 6 5.101370 3.995473 8.225509
#> 7 7.224315 5.748458 11.596157
#> 8 8.983472 7.213403 14.373850
#> 9 4.498345 3.502050 7.263311
#> 10 10.005028 8.067468 15.981979Notice how when ‘all’ is used, the output is a dataframe with the name of each value as the column names.
But what if we want to supply our own data and make predictions? Since the ‘flexsurv’ regression object generates parameters, we can add our own data points and get the residual life values for each.
Let’s start by generating some data. This data has the same two levels as ‘fsr2’ above (group, age). Note that the data is placed in a data frame. The new data must be inputted as a data frame for the function to work.
We run that through the reslifefsr() function and get a 3 row data frame
reslifefsr(obj = fsr2,life= 4,p = .6, type = 'all', newdata= newdata)
#> mean median percentile
#> 1 12.300164 9.992529 12.53526
#> 2 31.905635 26.564496 32.84741
#> 3 7.674191 6.122261 7.74986Now we can show off some features of the function, specifically how it handles different input data.
Lets change the order of the variables in the input data and rerun. It will show the exact same results as above.
newdata = data.frame(group, age)
reslifefsr(obj = fsr2,life= 4,p = .6, type = 'all', newdata= newdata)
#> mean median percentile
#> 1 12.300164 9.992529 12.53526
#> 2 31.905635 26.564496 32.84741
#> 3 7.674191 6.122261 7.74986This is a handy feature of the function that can enhance the user experience. Another feature of this is the ability to handle data with extra columns. Let’s show that feature below.
extra = c(100, 100, 100)
newdata2 = data.frame(age, group, extra)
reslifefsr(obj = fsr2,life= 4,p = .6, type = 'all', newdata= newdata2)
#> mean median percentile
#> 1 12.300164 9.992529 12.53526
#> 2 31.905635 26.564496 32.84741
#> 3 7.674191 6.122261 7.74986The extra column ‘extra’ is ignored by the function.
Are there scenarios where the package would return an error? Yes, in fact, there are a few.
The first is if not enough data is provided, which is shown below.
newdata3 = data.frame(group)
reslifefsr(obj = fsr2, life = 4,p = .6, type = 'all', newdata= newdata3)
#> Warning in if (fsoutput$covdata$covnames != colnames(newdata)) {: the condition
#> has length > 1 and only the first element will be used
#> mean median percentile
#> 1 87.22487 73.51266 90.10934
#> 2 150.23904 127.02587 155.29444
#> 3 47.56366 39.84435 49.06227Since only ‘group’ is provided, predictions cannot be generated. The function needs values for ‘age’ in order to make predictions.
Another scenario could happen with factor data. The levels for ‘group’ are ‘Good’, ‘Medium’, and ‘Poor’. But what if a level was inputted that isn’t one of those three? Let’s see.
group = c("Medium", 'Good', "Terrible")
newdata = data.frame(age, group)
reslifefsr(obj = fsr2,life= 4,p = .6, type = 'all', newdata= newdata)
#> Error in weibull_mlr(obj, life, p, type, newdata): Incorrect Level EnteredSince ‘Terrible’ is not a valid level, the function errors and returns a message notifying the user of this error.
It is also possible to plot the mean residual life for a set of values from a ‘flexsurv’ object.
We take the first line from newbc and create a blank vector of length 10. We then loop the function to create the MRLs for the first 10 integers, which gives us the expected survival time given the patient has survived y years.
We can then plot these values.
newdata4 = newbc[1,]
mrl = rep(0,10)
for (i in 1:10){
mrl[i] = reslifefsr(obj = fsr2,life= i,p = .6, type = 'mean', newdata= newdata4)
}
plot(x=c(1:10), y=mrl,type="l", xlab = 'Years Survived', ylab = 'Mean Residual Life', main = 'MRL over Time')
That is an overview of how to use the flexsurvreg() dependent function. Next, we will show how to get residual life values by inputting your own parameters.
