-
Notifications
You must be signed in to change notification settings - Fork 0
/
gompertz_class.R
144 lines (133 loc) · 5.51 KB
/
gompertz_class.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
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
## Gompertz survival objects for simulations
## by JJAV 20231128
#' Factory of SURVIVAL objects with Gompertz distributions
#'
#' Creates a SURVIVAL object with an Gompertz distribution.
#'
#' @section Parameters:
#'
#' To create an exponential survival object the following
#' options are available:
#'
#' _`scale`_ and _`shape`_ to specify the canonical parameter of the distribution, or
#'
#' _`surv`_, _`t`_ and _`shape`_ for the proportion surviving (no events) at time t and shape, or
#'
#' _`fail`_ and _`t`_ and _`shape`_ for the proportion failing (events) at time t and shape.
#'
#' scale = -log(surv)·shape/(exp(shape·t))
#'
#' scale = -log(1-fail)·shape/(exp(shape·t))
#'
#' The parameters should be spell correctly as partial matching is not available
#'
#' @param ... Parameters to define the distribution. See the Parameters for details
#' @return a SURVIVAL object of the Gompertz distribution family. See the
#' documentation of `s_factory` for the methods available for SURVIVAL objects
#' @export
#' @importFrom stats runif
#' @examples
#' s_gompertz(scale = 1, shape = 1.5)
#' s_gompertz(surv = 0.4, t = 2, shape = 1.5)
#' s_gompertz(fail = 0.6, t = 2, shape = 1.5)
s_gompertz <- function(...) {
params <- list(...)
nparam <- names(params)
# This function is the factory of the class
.factory_gompertz <- function(scale, shape) {
iCum_Hfx <- function(H){
stopifnot("Must be positive number" = all(H >= 0))
1/shape*log(shape/scale*H+1)
#log(shape/scale*H+1)/shape
}
structure(
list(
distribution = "GOMPERTZ",
params = list(scale = scale, shape = shape),
sfx = function(t) {
stopifnot("Must be positive number" = all(t >= 0))
exp(scale/shape*(1-exp(shape*t)))
#exp(-scale/shape*(exp(shape*t)-1))
},
hfx = function(t) {
stopifnot("Must be positive number" = all(t >= 0))
scale*exp(shape*t)
#scale*exp(shape*t)
},
Cum_Hfx = function(t) {
stopifnot("Must be positive number" = all(t >= 0))
scale/shape*(exp(shape*t)-1)
# scale/shape*(exp(shape*t)-1)
},
invCum_Hfx=iCum_Hfx,
rsurv = function(n){
stopifnot("Must be positive number" = all(n > 0))
iCum_Hfx(-log(runif(n)))
},
rsurvhr = function(hr){
stopifnot("hr must be numeric" = is.numeric(hr))
stopifnot("Must be positive numbers > 0" = all(hr > 0))
# Following Bender, Augustin and Blettner 2005
iCum_Hfx(-log(runif(length(hr)))/hr)
},
rsurvaft = function(aft){
stopifnot("aft must be numeric" = is.numeric(aft))
stopifnot("aft must be positive numbers > 0" = all(aft > 0))
iCum_Hfx(-log(runif(length(aft))))/aft
},
rsurvah = function(aft,hr){
stopifnot("aft must be numeric" = is.numeric(aft))
stopifnot("hr must be numeric" = is.numeric(hr))
stopifnot("aft and hr must be of the same length" = length(aft)==length(hr) )
stopifnot("aft must be positive numbers > 0" = all(aft > 0))
stopifnot("hr must be positive numbers > 0" = all(hr > 0))
iCum_Hfx(-log(runif(length(aft)))/hr)/aft
}
),
class = c("SURVIVAL")
)
}
# Definition based on scale and shape
if (length(nparam == 2) &
all(c("scale","shape") %in% nparam)) {
stopifnot("scale should be a single number" = is_single_number(params$scale))
stopifnot("scale must be greater than 0" = params$scale > 0)
stopifnot("shape should be a single number" = is_single_number(params$shape))
stopifnot("shape can't be 0 " = params$shape != 0)
return(.factory_gompertz(params$scale, params$shape))
}
# Definition based in proportion surviving, time and shape
if(
length(nparam == 3) &
all(c("surv","t","shape") %in% nparam)) {
stopifnot("surv must be a single number" = is_single_number(params$surv))
stopifnot("surv must be greater than 0" = params$surv > 0)
stopifnot("surv must be smaller than 1" = params$surv < 1)
stopifnot("t must be a single number" = is_single_number(params$t))
stopifnot("t must be greater than 0" = params$t > 0)
stopifnot("shape should be a single number" = is_single_number(params$shape))
stopifnot("shape can't be 0 " = params$shape != 0)
scale = -log(params$surv)*params$shape/(exp(params$shape*params$t)-1)
return(.factory_gompertz(scale, params$shape))
}
# Definition based on proportion failing and time
if(
length(nparam == 3) &
all(c("fail","t","shape") %in% nparam)) {
stopifnot("fail must be a single number" = is_single_number(params$fail))
stopifnot("fail must be greater than 0" = params$fail > 0)
stopifnot("fail must be lower than 1" = params$fail < 1)
stopifnot("t must be a single number" = is_single_number(params$t))
stopifnot("t must be greater than 0" = params$t > 0)
stopifnot("shape should be a single number" = is_single_number(params$shape))
stopifnot("shape can't be 0 " = params$shape != 0)
scale = -log((1-params$fail))*params$shape/(exp(params$shape*params$t)-1)
return(.factory_gompertz(scale, params$shape))
}
message(
"Valid parameters to define a Gompertz distribution are: \n",
"scale and shape: for the canonical parameters of the distribution, or\n",
"surv, t and shape: for the surviving proportion (no events) at time t and shape, or\n",
"fail, t and shape: for the failure proportion (events) at time t and shape, or \n")
stop("Not valid parameters")
}