-
Notifications
You must be signed in to change notification settings - Fork 0
/
importance_survreg.R
80 lines (74 loc) · 2.75 KB
/
importance_survreg.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
# Copyright 2019 Robert Carnell
#' Create a variable importance plot for a survreg model
#'
#' @inheritParams importance
#' @param model_data the data used to fit the model
#' @param dict a plotting dictionary for models terms
#' @param nperm the number of permutations used to calculate the importance
#'
#' @inherit importance return
#' @export
#'
#' @seealso \code{\link{importance}}
#'
#' @import survival
#' @import ggplot2
#'
#' @examples
#' model_final <- survival::survreg(survival::Surv(futime, fustat) ~ ecog.ps*rx + age,
#' data = survival::ovarian,
#' dist = "weibull")
#' imp <- importance(model_final, survival::ovarian, nperm = 50)
#' plot(imp)
importance.survreg <- function(model_final, model_data, dict=NA, nperm = 500, ...)
{
# model_final <- survival::survreg(survival::Surv(futime, fustat) ~ ecog.ps*rx + age, survival::ovarian, dist="weibull")
# model_data <- survival::ovarian
# nperm <- 100
# dict = NA
#
# w <- sample(1:7, size = nrow(survival::ovarian), replace = TRUE)
# model_final <- survival::survreg(survival::Surv(futime, fustat) ~ ecog.ps*rx + age, survival::ovarian, dist = "weibull", weights = w)
otherVariables <- list(...)
vars <- rownames(attr(formula(model_final), "factors"))[-1]
n <- nrow(model_data)
baseLikelihood <- model_final$loglik[2]
if (!is.null(model_final$weights)) {
w <- model_final$weights
}
# randomly permute each variable
importances <- numeric(length(vars))
for (j in seq_along(vars))
{
v <- vars[j]
model_data_new <- model_data
temp <- numeric(nperm)
for (i in 1:nperm)
{
model_data_new[, v] <- model_data_new[sample(1:n, n, replace = FALSE), v]
if (is.null(model_final$weights))
{
model_new <- survival::survreg(formula(model_final), model_data_new,
dist = model_final$dist)
} else
{
# Note that if we included "weights = model_final$weights" then the context
# does not travel into the function. We need to re-assign that vector
# to w before passing in
model_new <- survival::survreg(formula(model_final), model_data_new,
dist = model_final$dist,
weights = w)
}
temp[i] <- model_new$loglik[2]
}
importances[j] <- mean(temp)
}
importances_final <- 1 - exp(pmin(0, importances - baseLikelihood))
dat2 <- data.frame(variable = vars,
value = importances_final)
dat2 <- dat2[order(dat2$value, decreasing = TRUE),]
dat2$variable <- factor(dat2$variable, levels = rev(dat2$variable))
return(structure(list(data = dat2,
type = "survreg"),
class = "importance_plot"))
}