-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathstan_emax.R
executable file
·227 lines (195 loc) · 7.26 KB
/
stan_emax.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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
#' Bayesian Emax model fit with Stan
#'
#' Run sigmoidal Emax model fit with formula notation
#'
#' @export
#' @param formula a symbolic description of variables for Emax model fit.
#' @param data an optional data frame containing the variables in the model.
#' @param gamma.fix a (positive) numeric or NULL to specify gamma (Hill coefficient) in the sigmoidal Emax model.
#' If NULL, gamma will be estimated from the data.
#' If numeric, gamma is fixed at the number provided.
#' Default = 1 (normal Emax model).
#' @param e0.fix a numeric or NULL to specify E0 in the Emax model.
#' If NULL, E0 will be estimated from the data.
#' If numeric, E0 is fixed at the number provided.
#' Default = NULL (estimate from the data).
#' @param emax.fix a numeric or NULL to specify Emax in the Emax model.
#' If NULL, Emax will be estimated from the data.
#' If numeric, Emax is fixed at the number provided.
#' Default = NULL (estimate from the data).
#' @param priors a named list specifying priors of parameters (ec50, emax, e0, gamma, sigma).
#' Each list item should be length 2 numeric vector, one corresponding to mean and
#' another corresponding to standard deviation.
#' Currently only supports normal distribution for priors.
#' @param param.cov a named list specifying categorical covariates on parameters (ec50, emax, e0).
#' Convert a column into factor if specific order of covariates are needed.
#' @param ... Arguments passed to [rstan::sampling] (e.g. iter, chains).
#' @return An object of class `stanemax`
#' @details The following structure is used for the Emax model:
#' \deqn{Response = e_0 + e_{max} \times exposure ^{\gamma} / (ec50 ^{\gamma} + exposure ^ {\gamma}) + \epsilon}{Response = e0 + emax * exposure ^ gamma / (ec50 ^ gamma + exposure ^ gamma) + epsilon}
#' \deqn{\epsilon \sim N(0, \sigma^2)}{epsilon ~ N(0, sigma ^ 2)}
#'
#' @examples
#' \dontrun{
#' data(exposure.response.sample)
#' fit1 <- stan_emax(response ~ exposure,
#' data = exposure.response.sample,
#' # the next line is only to make the example go fast enough
#' chains = 1, iter = 500, seed = 12345
#' )
#' print(fit1)
#'
#' # Set priors manually, also estimate gamma instead of the default of fix to 1
#' fit2 <- stan_emax(response ~ exposure,
#' data = exposure.response.sample, gamma.fix = NULL,
#' priors = list(
#' ec50 = c(100, 30), emax = c(100, 30), e0 = c(10, 5),
#' gamma = c(0, 3), sigma = c(0, 30)
#' ),
#' # the next line is only to make the example go fast enough
#' chains = 1, iter = 500, seed = 12345
#' )
#' print(fit2)
#'
#' data(exposure.response.sample.with.cov)
#' # Specify covariates
#' fit3 <- stan_emax(
#' formula = resp ~ conc, data = exposure.response.sample.with.cov,
#' param.cov = list(emax = "cov2", ec50 = "cov3", e0 = "cov1"),
#' # the next line is only to make the example go fast enough
#' chains = 1, iter = 500, seed = 12345
#' )
#' print(fit3)
#' }
#'
# Remove NA data, show warning
stan_emax <- function(formula, data,
gamma.fix = 1, e0.fix = NULL, emax.fix = NULL,
priors = NULL, param.cov = NULL, ...) {
out.prep <- stan_emax_prep(formula, data, gamma.fix, e0.fix, emax.fix, param.cov)
standata <- set_prior(out.prep$standata, priors)
out <- stan_emax_run(stanmodels$emax, standata = standata, ...)
out$prminput <- out.prep$prminput
return(out)
}
#' Bayesian Emax model fit with Stan for binary endpoint
#'
#' @export
#' @inheritParams stan_emax
#' @examples
#' \dontrun{
#' data(exposure.response.sample.binary)
#' fit1 <- stan_emax_binary(
#' y ~ conc,
#' data = exposure.response.sample.binary,
#' # the next line is only to make the example go fast enough
#' chains = 2, iter = 500, seed = 12345
#' )
#' print(fit1)
#'
#' # Specify covariates
#' fit2 <- stan_emax_binary(
#' formula = y ~ conc, data = exposure.response.sample.binary,
#' param.cov = list(emax = "sex"),
#' # the next line is only to make the example go fast enough
#' chains = 2, iter = 500, seed = 12345
#' )
#' print(fit2)
#' }
stan_emax_binary <- function(
formula, data,
gamma.fix = 1, e0.fix = NULL, emax.fix = NULL,
priors = NULL, param.cov = NULL, ...) {
out.prep <- stan_emax_prep(formula, data, gamma.fix, e0.fix, emax.fix, param.cov)
standata <- set_prior(out.prep$standata, priors, endpoint_type = "binary")
out <- stan_emax_run(stanmodels$emax_binary,
standata = standata,
out_class_name = "stanemaxbin",
...
)
out$prminput <- out.prep$prminput
return(out)
}
# Parse formula and put together stan data object
stan_emax_prep <- function(formula, data,
gamma.fix = 1, e0.fix = NULL, emax.fix = NULL, param.cov = NULL) {
check_param_cov(param.cov)
cov.levels <- covs_get_levels(data, param.cov)
df.model <- create_model_frame(formula, data, param.cov, cov.levels)
standata <-
create_standata(df.model, gamma.fix, e0.fix, emax.fix)
out.prep <- list()
out.prep$standata <- standata
out.prep$prminput <- list()
out.prep$prminput$formula <- formula
out.prep$prminput$df.model <- df.model
out.prep$prminput$cov.levels <- cov.levels
out.prep$prminput$param.cov <- param.cov
return(out.prep)
}
# Check param.cov input
check_param_cov <- function(param.cov = NULL) {
if (sum(!(names(param.cov) %in% c("emax", "ec50", "e0")))) {
stop("Covariates can be specified only to emax, ec50, or e0")
}
}
# Run Emax model
# @param stanmodel a Stan model object.
# @param standata a data file for model fit
# @param ... Arguments passed to `rstan::sampling` (e.g. iter, chains).
# @return An object of class `stanemax`
#
stan_emax_run <- function(stanmodel, standata,
out_class_name = c("stanemax", "stanemaxbin"),
...) {
out_class_name <- match.arg(out_class_name)
# Run stan model and prepare `stanemax` object
stanfit <- rstan::sampling(stanmodel, data = standata, ...)
out <- list(
stanfit = stanfit,
standata = standata
)
structure(out, class = out_class_name)
}
create_standata <- function(df.model, gamma.fix = 1, e0.fix = NULL, emax.fix = NULL) {
out <- list(
exposure = df.model$exposure,
response = df.model$response,
covemax = as.numeric(df.model$covemax),
covec50 = as.numeric(df.model$covec50),
cove0 = as.numeric(df.model$cove0),
n_covlev_emax = length(levels(df.model$covemax)),
n_covlev_ec50 = length(levels(df.model$covec50)),
n_covlev_e0 = length(levels(df.model$cove0)),
N = nrow(df.model),
gamma_fix_flg = 1,
gamma_fix_value = 1,
e0_fix_flg = 0,
e0_fix_value = 0,
emax_fix_flg = 0,
emax_fix_value = 0
)
if (!is.null(gamma.fix) && !is.na(gamma.fix)) {
if (!is.numeric(gamma.fix)) stop("gamma.fix must be NULL or numeric")
if (gamma.fix <= 0) stop("gamma.fix must be NULL or positive numeric")
out$gamma_fix_flg <- 1
out$gamma_fix_value <- gamma.fix
} else {
out$gamma_fix_flg <- 0
}
if (!is.null(e0.fix) && !is.na(e0.fix)) {
if (!is.numeric(e0.fix)) stop("e0.fix must be NULL or numeric")
out$e0_fix_flg <- 1
out$e0_fix_value <- e0.fix
} else {
out$e0_fix_flg <- 0
}
if (!is.null(emax.fix) && !is.na(emax.fix)) {
if (!is.numeric(emax.fix)) stop("emax.fix must be NULL or numeric")
out$emax_fix_flg <- 1
out$emax_fix_value <- emax.fix
} else {
out$emax_fix_flg <- 0
}
return(out)
}