absoluteRisk is slow
Closed this issue · 7 comments
turgeonmaxime commented
This was brought up by @Jesse-Islam : absoluteRisk is pretty slow. The main reason, I think, is that we have a double loop: for each time point in time and each row in newdata, we are computing the survival function. I'm not sure how we can make it faster. Since we are repeatedly calling integrate, and it's written in Fortran, we may want to move the whole computation to compiled code...
Here's a reprex:
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
modelCB <- fitSmoothHazard(status ~ time + karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
time_vect <- seq(0,999,10)
microbenchmark::microbenchmark(
absoluteRisk(object = modelCB, time = time_vect,
newdata = veteran),
times = 3
)
#> Unit: seconds
#> expr min
#> absoluteRisk(object = modelCB, time = time_vect, newdata = veteran) 13.38638
#> lq mean median uq max neval
#> 13.44013 13.51187 13.49389 13.57461 13.65533 3Created on 2020-02-21 by the reprex package (v0.3.0)
sahirbhatnagar commented
There are a few things to consider here. Feel free to disagree:
- I looked at C++ libraries for integration. Armadillo has the trapezoidal integral. It takes as input vectors X and Y and computes the integral of Y with respect to spacing in X
- There's also the GSL library which has a series of numerical integration functions (including something similar to the Fortran code behind
stats::integrate) - As you already pointed out, the problem is not the integration itself. I believe it is the double loop, but also the repeated calls to the
predictfunction which is called by thelambdafunction. There is also a subset of adata.frameoccurring for each individual. - Trying to implement the exact same
absoluteRiskfunction in C++ is troublesome because it would involve somehow porting thelambdafunction which in turn calls thepredictfunction. It gets even more complicated withabsoluteRisk.glmnetwhich involves aformulafunction. I dont know C++ that well to figure all of this out. - With that being said, the speedups (if any), are likely to come from thinking harder about how we calculate the absolute risks in R.
- Do we need repeated calls to subset ? Is there a faster way to subset as shown here and here
- Is it possible to vectorize the predict function, i.e., pre-calculate all the values of lambda for each time interval. If the intervals are small enough, we can assume a linear function, and just take the area of the triangle?
- The most simple thing is obviously to reduce the number of points at which we calculate the absolute risk as shown in the following code:
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
library(splines)
plotRisks <- function(modelCB,
number_of_intervals = c(10, 50, 100, 200, 300),
min_time = 0, max_time = 999){
Risks <- lapply(number_of_intervals, function(i)
absoluteRisk(object = modelCB,
time = seq(min_time,max_time,i),
nsamp = 100,
newdata = veteran,
method = "numerical")
)
par(mar = c(4.1, 4.0, 3.1, 1.1), mfrow = c(2,3))
lapply(Risks, function(i) matplot(x = i[,1],y = i[,-1], ylim = c(0,1.1), type = "l",
ylab = "Risk", xlab = "Time", main = sprintf("Number of intervals:%0.0f",
i[,1][2]-i[,1][1])))
plot(x = seq(0,1000, length.out = 50), y = seq(0,1.1, length.out = 50), type = "n",
ylab = "Risk", xlab = "Time", main = "Mean Risk (averaged over individuals")
lapply(seq_along(Risks), function(i) {
lines(x = Risks[[i]][,1], y = rowMeans(Risks[[i]][,-1]), col = i, lwd = 1.5)
})
legend("bottomright", legend = paste0("# intervals = ",number_of_intervals),
col = seq_along(number_of_intervals), lwd = 1.5)
}
modelCB_linear <- fitSmoothHazard(status ~ time + karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
modelCB_splines <- fitSmoothHazard(status ~ bs(time, 8) + karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
modelCB_splines_interaction <- fitSmoothHazard(status ~ bs(time, 8)*karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
plotRisks(modelCB = modelCB_linear)
plotRisks(modelCB = modelCB_splines)
plotRisks(modelCB = modelCB_splines_interaction)
Created on 2020-02-22 by the reprex package (v0.3.0)
turgeonmaxime commented
I agree that if we can avoid going to compiled code, we should. I also agree that we should be able to speed up the R code. I don't think subset is really the main problem: we can probably speed absoluteRisk up a bit by extracting the row, saving it to a variable current_obs, and compute the integral at all time points. But the gains should be minimal, since we never really compute the absolute risk over more than a few hundred time points. But maybe I'm underestimating it.
More generally, I feel like if we shouldn't have to repeat the integrate call over all time points... Perhaps if we had a good coverage of the time axis between zero and the maximum time point in time, we could reuse all evaluations of the hazard function to compute all integrals at once (e.g. through cumsum). As a bonus, we would end up with only one call to predict; on the downside, we may lose accuracy if we don't pick the time points correctly.
I'll make the change with subset described in the first paragraph, since it's easy to do. We can see after that if we need to make any further improvements.
turgeonmaxime commented
As I expected, subsetting once for all time points doesn't change the computation time by much:
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
modelCB <- fitSmoothHazard(status ~ time + karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
time_vect <- seq(0,999,10)
microbenchmark::microbenchmark(
absoluteRisk(object = modelCB, time = time_vect,
newdata = veteran),
times = 3
)
#> Unit: seconds
#> expr min
#> absoluteRisk(object = modelCB, time = time_vect, newdata = veteran) 12.85449
#> lq mean median uq max neval
#> 13.00463 13.37788 13.15477 13.63958 14.12439 3Created on 2020-02-24 by the reprex package (v0.3.0)
turgeonmaxime commented
Here's a proof of concept for computing all integrals (at a given observation) as cumulative sum:
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
modelCB <- fitSmoothHazard(status ~ time + karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
time_vect <- seq(0,999,10)
# New proposal
estimate_risk_newtime2 <- function(object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100) {
time_ordered <- unique(c(0, sort(time)))
# Compute points at which we evaluate integral
knots <- seq(0, max(time_ordered),
length.out = (length(time_ordered) - 1) * nsamp)
knots <- unique(sort(c(knots, time_ordered)))
# Create matrix to store output
output <- matrix(NA, ncol = nrow(newdata) + 1, nrow = length(time_ordered))
output[,1] <- time_ordered
colnames(output) <- c("time", rep("", nrow(newdata)))
rownames(output) <- rep("", length(time_ordered))
output[1,-1] <- 0
for (j in 1:nrow(newdata)) {
# Extract current obs
current_obs <- newdata[j,,drop = FALSE]
newdata2 <- data.frame(current_obs, offset = rep_len(0, length(knots)),
row.names = as.character(1:length(knots)))
newdata2[modelCB$timeVar] <- knots
pred <- predict(modelCB, newdata2)
# Compute integral using trapezoidal rule
output[,j + 1] <- pracma::cumtrapz(knots, exp(pred))[knots %in% c(0, time_vect)]
}
# Transform to survival scale
output[,-1] <- exp(-output[,-1])
output[,-1] <- 1 - output[,-1]
return(output)
}
# Compare output of two approaches
lambda <- function(x, fit = modelCB, newdata = veteran) {
# Note: the offset should be set to zero when estimating the hazard.
newdata2 <- data.frame(newdata, offset = rep_len(0, length(x)),
row.names = as.character(1:length(x)))
newdata2[fit$timeVar] <- x
pred <- predict(fit, newdata2)
return(as.numeric(exp(pred)))
}
foo <- casebase:::estimate_risk_newtime(lambda, object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100)
bar <- estimate_risk_newtime2(object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100)
all.equal(foo, bar)
#> [1] TRUE
# Compare computation time
microbenchmark::microbenchmark(
"current" = casebase:::estimate_risk_newtime(lambda, object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100),
"proposal" = estimate_risk_newtime2(object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100),
times = 3
)
#> Unit: milliseconds
#> expr min lq mean median uq max
#> current 12361.9065 12381.5055 12413.9866 12401.1045 12440.0267 12478.9489
#> proposal 747.8593 755.7452 784.2248 763.6311 802.4075 841.1839
#> neval cld
#> 3 b
#> 3 aCreated on 2020-02-24 by the reprex package (v0.3.0)
turgeonmaxime commented
As we can see, it's much faster, and it seems to be just as accurate (although time is modelled linearly in this example, so we would need to check more complicated functions). There are most likely better ways of choosing the knots, but this seems to be a good first approximation. @sahirbhatnagar what do you think?
sahirbhatnagar commented
Wow! That's alot faster now.
I tried running a model with splines and it seemed to work as well. However, a warning is issued because I guess 0 will always be outside the boundary for the basis functions. You can get the boundary knots used from attr(modelCB$terms, "predvars")[[3]]$Boundary.knots, and then limit the knots to be inside these points. But you would have have to figure out a way to do this programatically. Or maybe just surpress these warnings, or let the user know that this shouldn't be a problem.
Is there actually any difference in terms of the new and old function? Since within each knot (for sufficiently small distances between knots), the hazard is constant, the trapzoidal rule should give an exact answer.
If you still think there is a difference, then maybe keep the old option (call it exact ?), but default to the new one.
Also, in the documentation for casebase::absoluteRisk, method argument says that it defaults to montecarlo. This should be changed to numerical.
library(casebase)
#> See example usage at http://sahirbhatnagar.com/casebase/
library(survival)
library(splines)
modelCB <- fitSmoothHazard(status ~ bs(time) + karno + diagtime + age + prior +
celltype + trt, data = veteran, ratio = 100)
#> 'time' will be used as the time variable
time_vect <- seq(0,999,10)
# New proposal
estimate_risk_newtime2 <- function(object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100) {
time_ordered <- unique(c(0, sort(time)))
# Compute points at which we evaluate integral
knots <- seq(0, max(time_ordered),
length.out = (length(time_ordered) - 1) * nsamp)
knots <- unique(sort(c(knots, time_ordered)))
# Create matrix to store output
output <- matrix(NA, ncol = nrow(newdata) + 1, nrow = length(time_ordered))
output[,1] <- time_ordered
colnames(output) <- c("time", rep("", nrow(newdata)))
rownames(output) <- rep("", length(time_ordered))
output[1,-1] <- 0
for (j in 1:nrow(newdata)) {
# Extract current obs
current_obs <- newdata[j,,drop = FALSE]
newdata2 <- data.frame(current_obs, offset = rep_len(0, length(knots)),
row.names = as.character(1:length(knots)))
newdata2[modelCB$timeVar] <- knots
# browser()
# modelCB$coefficients
# attr(modelCB$terms, "predvars")[[3]]$Boundary.knots
pred <- predict(modelCB, newdata2)
# Compute integral using trapezoidal rule
output[,j + 1] <- pracma::cumtrapz(knots, exp(pred))[knots %in% c(0, time_vect)]
}
# Transform to survival scale
output[,-1] <- exp(-output[,-1])
output[,-1] <- 1 - output[,-1]
return(output)
}
# Compare output of two approaches
lambda <- function(x, fit = modelCB, newdata = veteran) {
# Note: the offset should be set to zero when estimating the hazard.
newdata2 <- data.frame(newdata, offset = rep_len(0, length(x)),
row.names = as.character(1:length(x)))
newdata2[fit$timeVar] <- x
pred <- predict(fit, newdata2)
return(as.numeric(exp(pred)))
}
foo <- casebase:::estimate_risk_newtime(lambda, object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100)
bar <- estimate_risk_newtime2(object = modelCB, time = time_vect,
newdata = veteran, method = "numerical",
nsamp = 100)
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
#> Warning in bs(time, degree = 3L, knots = numeric(0), Boundary.knots =
#> c(0.00794587563723326, : some 'x' values beyond boundary knots may cause ill-
#> conditioned bases
all.equal(foo, bar)
#> [1] TRUE
par(mfrow = c(1,2))
matplot(foo[,1], foo[,-1], type = "l")
matplot(bar[,1], bar[,-1], type = "l")
Created on 2020-02-24 by the reprex package (v0.3.0)
turgeonmaxime commented
Oh yes, those warnings! I have a way to turn them off in the code, I just removed that part when prototyping the new code.
Thanks for checking with a spline model! I feel pretty confident that I can switch the code now to this faster version. I'll add this to my todo list, it shouldn't take too long.