/
evaluate_fim_map.R
179 lines (156 loc) · 7.17 KB
/
evaluate_fim_map.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
#' Compute the Bayesian Fisher information matrix
#'
#' Computation of the Bayesian Fisher information matrix for
#' individual parameters of a population model based on
#' Maximum A Posteriori (MAP) estimation of the empirical Bayes estimates (EBEs) in a population model
#'
#' @param poped.db A PopED database
#' @param num_sim_ids If \code{use_mc=TRUE}, how many individuals should be
#' simulated to make the computations.
#' @param use_mc Should the calculation be based on monte-carlo simulations. If
#' not then then a first order approximation is used
#' @param use_purrr If \code{use_mc=TRUE} then should the method use the package
#' purrr in calculations? This may speed up computations (potentially).
#' @param shrink_mat Should the shrinkage matrix be returned. Calculated as the
#' inverse of the Bayesian Fisher information matrix times the inverse of the
#' omega matrix (variance matrix of the between-subject variability).
#' @return The Bayesian Fisher information matrix for each design group
#' @export
#'
#' @example tests/testthat/examples_fcn_doc/warfarin_basic.R
#' @example tests/testthat/examples_fcn_doc/examples_shrinkage.R
#' @references \enumerate{
#' \item Combes, F. P., Retout, S.,
#' Frey, N., & Mentre, F. (2013). Prediction of shrinkage of individual
#' parameters using the Bayesian information matrix in non-linear mixed effect
#' models with evaluation in pharmacokinetics. Pharmaceutical Research, 30(9),
#' 2355-67. \doi{10.1007/s11095-013-1079-3}.
#' \item Hennig, S., Nyberg, J., Fanta, S., Backman, J.
#' T., Hoppu, K., Hooker, A. C., & Karlsson, M. O. (2012). Application of the
#' optimal design approach to improve a pretransplant drug dose finding design
#' for ciclosporin. Journal of Clinical Pharmacology, 52(3), 347-360.
#' \doi{10.1177/0091270010397731}.
#' }
evaluate_fim_map <- function(poped.db,
use_mc = FALSE,
num_sim_ids = 1000,
use_purrr = FALSE,
shrink_mat=F){
# if (poped.db$design$m > 1) {
# warning("Shrinkage should only be computed for a single arm, please adjust your script accordingly.")
# }
group <- type <- NULL
# tranform random effects to fixed effects
tmp_fg <- poped.db$model$fg_pointer
if(is.character(tmp_fg)) tmp_fg <- eval(parse(text=tmp_fg))
largest_bpop <- find.largest.index(tmp_fg,lab = "bpop")
largest_b <- find.largest.index(tmp_fg,lab = "b")
largest_eps <- find.largest.index(poped.db$model$ferror_pointer,lab = "epsi",mat=T,mat.row = F)
replace_fun <- function(txt){
#extract number
num <- stringr::str_extract(txt,"\\d+") %>% as.numeric() %>% "+"(.,largest_bpop)
stringr::str_replace(txt,"\\d+",paste0(num))
}
txt <- capture.output(body(tmp_fg))
txt <- stringr::str_replace_all(txt,"b\\[(\\d+)",replace_fun)
txt <- stringr::str_replace_all(txt,"b\\[","bpop\\[")
env <- environment()
body(tmp_fg,envir=env) <- parse(text=txt)
#environment(sfg_tmp) <- env
# fix population parameters
# add IIV and IOV as fixed effects
# change to one individual
# TODO: need to add iov
extra_bpop <- matrix(0,nrow = largest_b,ncol = 3)
bpop_tmp <-rbind(poped.db$parameters$bpop,extra_bpop)
notfixed_d <- poped.db$parameters$notfixed_d
non_zero_d <- poped.db$parameters$d[,2]!=0
notfixed_d_new <- notfixed_d | non_zero_d
notfixed_bpop_tmp <- c(rep(0,largest_bpop),notfixed_d_new)
poped.db_sh <- create.poped.database(poped.db,
fg_fun=tmp_fg,
bpop=bpop_tmp,
nbpop=nrow(bpop_tmp),
notfixed_bpop = notfixed_bpop_tmp,
d=matrix(nrow = 0,ncol = 3),
notfixed_d = matrix(nrow = 0,ncol = 3),
NumRanEff = 0,
notfixed_sigma = c(rep(0,largest_eps)),
groupsize=1,
mingroupsize = 1,
mintotgroupsize = 1)
# Compute FIM_map
fulld = getfulld(poped.db$parameters$d[notfixed_d_new,2,drop=F],poped.db$parameters$covd)
# get just groupwise values as well
db_list <- list()
#db_list[["all"]] <- poped.db_sh
num_groups <- poped.db_sh$design$m
for(i in 1:num_groups){
tmp_design <- poped.db_sh$design
tmp_design$m <- 1
for(nam in names(tmp_design)[names(tmp_design)!="m"]){
tmp_design[[nam]] <- tmp_design[[nam]][i,,drop=F]
}
tmp_db <- poped.db_sh
tmp_db$design <- tmp_design
db_list[[paste0("grp_",i)]] <- tmp_db
}
#out_list <- list()
out_df <- c()
#for(i in 1:1){
for(i in 1:length(db_list)){
poped.db_sh <- db_list[[i]]
if(use_mc){
# create list of eta values from pop model
if(any(size(fulld)==0)){
b_sim_matrix = zeros(num_sim_ids,0)
} else {
b_sim_matrix = mvtnorm::rmvnorm(num_sim_ids,sigma=fulld)
}
# create list of occasion eta values from pop model
NumOcc=poped.db$parameters$NumOcc
if(NumOcc!=0) warning("Shrinkage not computed for occasions\n")
fulldocc = getfulld(poped.db$parameters$docc[,2,drop=F],poped.db$parameters$covdocc)
bocc_sim_matrix = zeros(num_sim_ids*NumOcc,length(poped.db$parameters$docc[,2,drop=F]))
if(nrow(fulldocc)!=0) bocc_sim_matrix = mvtnorm::rmvnorm(num_sim_ids*NumOcc,sigma=fulldocc)
#now use these values to compute FIMmap
if(use_purrr){
fim_mean <- b_sim_matrix %>% t() %>% data.frame() %>% as.list() %>%
purrr::map(function(x){
x_tmp <- matrix(0,nrow = largest_b,ncol = 3)
x_tmp[,2] <- t(x)
bpop_tmp <-rbind(poped.db$parameters$bpop,x_tmp)
poped.db_sh_tmp <- create.poped.database(poped.db_sh, bpop=bpop_tmp)
evaluate.fim(poped.db_sh_tmp)
}) %>% simplify2array() %>% apply(1:2, mean)
} else {
fim_tmp <- 0
for(i in 1:num_sim_ids){
x_tmp <- matrix(0,nrow = largest_b,ncol = 3)
x_tmp[,2] <- t(b_sim_matrix[i,])
bpop_tmp <-rbind(poped.db$parameters$bpop,x_tmp)
poped.db_sh_tmp <- create.poped.database(poped.db_sh, bpop=bpop_tmp)
fim_tmp <- fim_tmp+evaluate.fim(poped.db_sh_tmp)
}
fim_mean <- fim_tmp/num_sim_ids
}
fim_map <- fim_mean + inv(fulld)
} else {
fim_map <- evaluate.fim(poped.db_sh) + inv(fulld)
}
poped_db_tmp <- poped.db
poped_db_tmp$parameters$notfixed_d <- notfixed_d_new
parnam <- get_parnam(poped_db_tmp)
rownames(fim_map) <- parnam[grepl("^(D\\[|d\\_)",parnam)]
colnames(fim_map) <- parnam[grepl("^(D\\[|d\\_)",parnam)]
out_df_tmp <- tibble::tibble(group=names(db_list[i]),fim_map=list(fim_map))
if(shrink_mat){
shrink_m <- inv(fim_map)%*%inv(fulld)
dimnames(shrink_m) <- dimnames(fim_map)
out_df_tmp$shrink_mat <- list(shrink_m)
}
out_df <- rbind(out_df,out_df_tmp)
#out_list[[names(db_list[i])]] <- list(shrinkage_var=shrink,shrinkage_sd=shrink_sd, rse=rse)
}
return(out_df)
}