#' Plot the efficiency of windows
#'
#' Function plots the efficiency of windows around the sample time points.
#' The function samples from a uniform distribution around the sample time
#' points for each group (or each individual with \code{deviate_by_id=TRUE},
#' with slower calculation times) and compares the results with the
#' design defined in \code{poped.db}. The maximal and minimal allowed values for all design variables as
#' defined in poped.db are respected (e.g. poped.db$design_space$minxt and
#' poped.db$design_space$maxxt).
#'
#' @param poped.db A poped database
#' @param iNumSimulations The number of design simulations to make within the
#' specified windows.
#' @inheritParams Doptim
#' @inheritParams create.poped.database
#' @param ... Extra arguments passed to \code{evaluate.fim}
#' @param xt_windows The distance on one direction from the optimal sample
#' times. Can be a number or a matrix of the same size as the xt matrix found
#' in \code{poped.db$design$xt}.
#' @param xt_plus The upper distance from the optimal sample times (xt +
#' xt_plus). Can be a number or a matrix of the same size as the xt matrix
#' found in \code{poped.db$design$xt}.
#' @param xt_minus The lower distance from the optimal sample times (xt -
#' xt_minus). Can be a number or a matrix of the same size as the xt matrix
#' found in \code{poped.db$design$xt}.
#' @param y_eff Should one of the plots created have efficiency on the y-axis?
#' @param y_rse Should created plots include the relative standard error of each
#' parameter as a value on the y-axis?
#' @param mean_line Should a mean value line be created?
#' @param mean_color The color of the mean value line.
#' @param deviate_by_id Should the computations look at deviations per
#' individual instead of per group?
#' @param parallel Should we use parallel computations (T/F)?
#' Other options can be defined in this function and passed
#' to \code{\link{start_parallel}}. See especially
#' the options \code{dlls} and \code{mrgsolve_model} from that function
#' if you have a model defined with
#' compiled code or are using mrgsolve.
#' @param seed The random seed to use.
#' @return A \link[ggplot2]{ggplot} object.
#'
#' @family evaluate_design
#' @family Simulation
#' @family Graphics
#'
#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
#' @example
#' tests/testthat/examples_fcn_doc/examples_plot_efficiency_of_windows.R
#' @export
#' @import ggplot2
plot_efficiency_of_windows <- function(poped.db,
xt_windows=NULL,
xt_plus=xt_windows,
xt_minus=xt_windows,
iNumSimulations=100,
y_eff = TRUE,
y_rse = TRUE,
ofv_calc_type = poped.db$settings$ofv_calc_type,
mean_line=TRUE,
mean_color="red",
#a_windows=NULL,x_windows=NULL,
#d_switch=poped.db$settings$d_switch,
deviate_by_id=FALSE,
parallel=F,
#parallel_type=NULL,
#num_cores = NULL,
#mrgsolve_model=NULL,
seed=round(runif(1,0,10000000)),
...){
if(!is.null(seed)) set.seed(seed)
design = poped.db$design
design_space = poped.db$design_space
#design$xt <- poped.db$gxt
#design$x <- poped.db$gx
#design$a <- poped.db$design$a
p = get_fim_size(poped.db) #%Size of FIM
xt_val = design$xt
a_val = design$a
x_val = design$x
ref_fmf <- evaluate.fim(poped.db,xt=xt_val,a=a_val,x=x_val,...)
#NumRows = size(xt_windows,1)*size(xt_windows,2)/2+size(a_windows,1)*size(a_windows,2)/2+size(x_windows,1)*size(x_windows,2)/2+p+1
if(y_eff){
eff = zeros(1,iNumSimulations)
# if(ofv_calc_type==4) {
# d_eff = zeros(1,iNumSimulations)
# }
}
if(y_rse) rse <- zeros(iNumSimulations,p)
fulld = getfulld(design$d[,2],design$covd)
fulldocc = getfulld(design$docc[,2,drop=F],design$covdocc)
if(!is.null(xt_windows) || !is.null(xt_minus) || !is.null(xt_plus)){
xt_val_min <- xt_val
xt_val_max <- xt_val
if(!is.null(xt_minus)) xt_val_min <- xt_val - xt_minus
if(!is.null(xt_plus)) xt_val_max <- xt_val + xt_plus
xt_val_min[xt_val_min<design_space$minxt] = design_space$minxt[xt_val_min<design_space$minxt]
xt_val_max[xt_val_max>design_space$maxxt] = design_space$maxxt[xt_val_max>design_space$maxxt]
}
# ------------ generate and evaluate new parameter sets
gen_eff <- function (it=1,...) {
if(!is.null(xt_windows) || !is.null(xt_minus) || !is.null(xt_plus)){
if(deviate_by_id){
for(grp in 1:size(xt_val)[1]){
groupsize <- design$groupsize[grp]
MS = design$model_switch
xt_tmp <- matrix(xt_val[grp,],nrow = groupsize,ncol = length(xt_val[grp,]),byrow = T)
a_tmp <- matrix(a_val[grp,],nrow = groupsize,ncol = length(a_val[grp,]),byrow = T)
x_tmp <- matrix(x_val[grp,],nrow = groupsize,ncol = length(x_val[grp,]),byrow = T)
MS_tmp <- matrix(MS[grp,],nrow = groupsize,ncol = length(MS[grp,]),byrow = T)
xt_val_max_tmp <- matrix(xt_val_max[grp,],nrow = groupsize,ncol = length(xt_val_max[grp,]),byrow = T)
xt_val_min_tmp <- matrix(xt_val_min[grp,],nrow = groupsize,ncol = length(xt_val_min[grp,]),byrow = T)
xt_tmp <- rand(size(xt_tmp,1),size(xt_tmp,2))*abs(xt_val_max_tmp - xt_val_min_tmp)+xt_val_min_tmp
grouped_xt <- design_space$G_xt[grp,]
if(anyDuplicated(grouped_xt)){
for(k in unique(grouped_xt)){
tmp.lst <- xt_tmp[,grouped_xt==k]
if(size(tmp.lst)[2]!=1){
tmp <- sample(1:size(tmp.lst)[2],1)
xt_tmp[,grouped_xt==k]=tmp.lst[,tmp]
}
}
}
if(grp==1){
xt_new <- xt_tmp
x_new <- x_tmp
a_new <- a_tmp
MS_new <- MS_tmp
} else {
xt_new <- rbind(xt_new,xt_tmp)
x_new <- rbind(x_new,x_tmp)
a_new <- rbind(a_new,a_tmp)
MS_new <- rbind(MS_new,MS_tmp)
}
}
design_new <- create_design(xt=xt_new,
groupsize=1,
x=x_new,
a=a_new,
model_switch=MS_new)
poped.db_new <- poped.db
poped.db_new$design <- design_new
fmf <- evaluate.fim(poped.db_new,...)
} else {
xt_new <- rand(size(xt_val,1),size(xt_val,2))*abs(xt_val_max - xt_val_min)+xt_val_min
grouped_xt <- design_space$G_xt
if(anyDuplicated(grouped_xt)){
#if((poped.db$design_space$bUseGrouped_xt)){
for(j in 1:size(xt_val,1) ){
for(k in unique(grouped_xt[j,])){
#for(k in min(poped.db$design_space$G_xt[j,]):max(poped.db$design_space$G_xt[j,])){
tmp.lst <- xt_new[j,design_space$G_xt[j,]==k]
if(length(tmp.lst)!=1){
tmp <- sample(tmp.lst,1)
xt_new[j,design_space$G_xt[j,]==k]=tmp
}
}
}
}
xt_val=xt_new
fmf <- evaluate.fim(poped.db,xt=xt_val,...)
}
}
# for(j in 1:size(xt_val,1) ){#Get simulated xt
# for(k in 1:size(xt_val,2)){
# xt_val[j,k] = xt_val[j,k] + rand*(xt_val_min[j,k]-xt_val_max(j,(k-1)*2+1))
# }
# }
# for(j in 1:size(a_windows,1) ){#Get simulated a
# for(k in 1:size(a_windows,2)/2){
# a_val(j,k) = a_windows(j,(k-1)*2+1) + rand*(a_windows(j,(k)*2)-a_windows(j,(k-1)*2+1))
# }
# }
# for(j in 1:size(x_windows,1) ){#Get simulated x
# for(k in 1:size(x_windows,2)/2){
# x_val(j,k) = x_windows(j,(k-1)*2+1) + rand*(x_windows(j,(k)*2)-x_windows(j,(k-1)*2+1))
# }
# }
# if((poped.db$design_space$bUseGrouped_xt)){
# xt_val=group_matrix(xt_val,poped.db$design_space$G_xt)
# }
# if((poped.db$design_space$bUseGrouped_a)){
# a_val=group_matrix(a_val,poped.db$design_space$G_a)
# }
# if((poped.db$design_space$bUseGrouped_x)){
# x_val=group_matrix(x_val,poped.db$design_space$G_x)
# }
#
# if((!bParallel)){
#if((poped.db$settings$d_switch)){
# returnArgs <- mftot(model_switch,poped.db$design$groupsize,ni,xt_val,x_val,a_val,bpop[,2],fulld,design$sigma,fulldocc,poped.db)
# fmf <- returnArgs[[1]]
# poped.db <- returnArgs[[2]]
eff <- NULL
if(y_eff){
tmp1 <- ofv_criterion(ofv_fim(fmf,poped.db,ofv_calc_type=ofv_calc_type),
p,poped.db,ofv_calc_type=ofv_calc_type)
tmp2 <- ofv_criterion(ofv_fim(ref_fmf,poped.db,ofv_calc_type=ofv_calc_type),
p,poped.db,ofv_calc_type=ofv_calc_type)
eff = tmp1/tmp2*100
# if(ofv_calc_type==4) {
# tmp1 <- ofv_criterion(ofv_fim(fmf,poped.db,ofv_calc_type=1),
# p,poped.db,ofv_calc_type=1)
# tmp2 <- ofv_criterion(ofv_fim(ref_fmf,poped.db,ofv_calc_type=1),
# p,poped.db,ofv_calc_type=1)
# d_eff[1,i] = tmp1/tmp2
# }
}
rse <- NULL
if(y_rse) rse <- get_rse(fmf,poped.db)
# } else {
# eff(i) = (ofv_fim(fmf,poped.db)/optdetfim)
# }
# } else {
# returnArgs <- ed_mftot(model_switch,poped.db$design$groupsize,ni,xt_val,x_val,a_val,bpop,d,poped.db$parameters$covd,design$sigma,poped.db$parameters$docc,poped.db)
# fmf <- returnArgs[[1]]
# ED_det <- returnArgs[[2]]
# poped.db <- returnArgs[[3]]
# if((bNormalizedEfficiency)){
# eff(i) = ofv_efficiency(ofv_criterion(ED_det,p,poped.db),ofv_criterion(optdetfim,p,poped.db))
# } else {
# eff(i) = (ED_det/optdetfim)
# }
# }
#
# if((bStoreInFile)){
# tmp_xt= reshape_matlab(xt_val',size(xt_windows,1)*size(xt_windows,2)/2,1)
# tmp_a = reshape_matlab(a_val',size(a_windows,1)*size(a_windows,2)/2,1)
# tmp_x = reshape_matlab(x_val',size(x_windows,1)*size(x_windows,2)/2,1)
# fileData(i,1:length(tmp_xt)) = tmp_xt
# fileData(i,1+length(tmp_xt):length(tmp_xt)+length(tmp_a)) = tmp_a
# fileData(i,1+length(tmp_a)+length(tmp_xt):length(tmp_xt)+length(tmp_a)+length(tmp_x)) = tmp_x
# try
# imf = inv(fmf)
# returnArgs <- get_cv(diag_matlab(imf),bpop,d,poped.db$parameters$docc,design$sigma,poped.db)
# params <- returnArgs[[1]]
# params_cv <- returnArgs[[2]]
# catch
# params_cv = zeros(0,p)
# }
# fileData(i,1+length(tmp_a)+length(tmp_xt)+length(tmp_x):length(tmp_a)+length(tmp_xt)+length(tmp_x)+p)=params_cv
# fileData(i,NumRows) = eff(i)
# }
# } else {
# designsin = update_designinlist(designsin,poped.db$design$groupsize,ni,xt_val,x_val,a_val,i,0)
# }
# }
#
# if((bParallel) ){#Execute everything in parallel instead
# designsout = execute_parallel(designsin,poped.db)
# for(i in 1:iNumSimulations){
# ED_det=designsout[[i]].ofv
# fmf=designsout[[i]]$FIM
#
# if((bNormalizedEfficiency)){
# eff(i) = ofv_efficiency(ofv_criterion(ED_det,p,poped.db),ofv_criterion(optdetfim,p,poped.db))
# } else {
# eff(i) = (ED_det/optdetfim)
# }
#
# if((bStoreInFile)){
# tmp_xt= reshape_matlab(designsin[[i]].xt',size(xt_windows,1)*size(xt_windows,2)/2,1)
# tmp_a = reshape_matlab(designsin[[i]].a',size(a_windows,1)*size(a_windows,2)/2,1)
# tmp_x = reshape_matlab(designsin[[i]].x',size(x_windows,1)*size(x_windows,2)/2,1)
# fileData(i,1:length(tmp_xt)) = tmp_xt
# fileData(i,1+length(tmp_xt):length(tmp_xt)+length(tmp_a)) = tmp_a
# fileData(i,1+length(tmp_a)+length(tmp_xt):length(tmp_xt)+length(tmp_a)+length(tmp_x)) = tmp_x
# try
# imf = inv(fmf)
# returnArgs <- get_cv(diag_matlab(imf),bpop,d,poped.db$parameters$docc,design$sigma,poped.db)
# params <- returnArgs[[1]]
# params_cv <- returnArgs[[2]]
# catch
# params_cv = zeros(0,p)
# }
# fileData(i,1+length(tmp_a)+length(tmp_xt)+length(tmp_x):length(tmp_a)+length(tmp_xt)+length(tmp_x)+p)=params_cv
# fileData(i,NumRows) = eff(i)
# }
# }
# }
#
#
# min_eff = min(eff)
# max_eff = max(eff)
# mean_eff = mean(eff)
#
# if((bStoreInFile)){
# csvwrite(strFileName,fileData)
# }
#
# if((bPlot)){
# figure
# ##hold on
# if((bNormalizedEfficiency)){
# title('Normalized Efficiency of samples from sampling/covariate - windows')
# } else {
# title('Efficiency of samples from sampling/covariate - windows')
# }
# ylim(matrix(c(100*min(eff)-10,max(100,100*max(eff))),nrow=1,byrow=T))
# plot(1:iNumSimulations,eff*100,'-b')
# xlabel('Sample number')
# if((bNormalizedEfficiency)){
# ylabel('Normalized Efficiency (%)')
# } else {
# ylabel('Efficiency (%)')
# }
# ##hold off
# }
#
# return( eff min_eff mean_eff max_eff )
return(c(Efficiency=eff,rse))
}
if(parallel){
parallel <- start_parallel(parallel,
seed=seed,
#parallel_type=parallel_type,
#num_cores=num_cores,
#mrgsolve_model=mrgsolve_model,
...)
on.exit(if(parallel && (attr(parallel,"type")=="snow")) parallel::stopCluster(attr(parallel,"cluster")))
}
it_seq <- 1:iNumSimulations
if(parallel && (attr(parallel,"type")=="multicore")){
res <- parallel::mclapply(it_seq,gen_eff,mc.cores=attr(parallel, "cores"))
} else if(parallel && (attr(parallel,"type")=="snow")){
res <- parallel::parLapply(attr(parallel, "cluster"),it_seq,gen_eff)
} else {
res <- lapply(it_seq,gen_eff)
}
df <- tibble::as_tibble(t(mapply(c,res,drop=0)))
df$drop <- NULL
df$sample <- it_seq
values <- NULL
ind <- NULL
#parameter_names_ff <- codetools::findGlobals(eval(parse(text=poped.db$model$ff_pointer)),merge=F)$variables
df_stack <- cbind(df$sample,stack(df,select=-sample))
names(df_stack) <- c("sample","values","ind")
if(y_eff){
levs <- levels(df_stack$ind)
# if(ofv_calc_type==4) {
# df_stack$ind <- factor(df_stack$ind,levels=c("D-Efficiency","Efficiency",levs[-c(grep("Efficiency",levs))]))
# } else {
df_stack$ind <- factor(df_stack$ind,levels=c("Efficiency",levs[-c(grep("Efficiency",levs))]))
# }
#levels(df_stack$ind) <- c("Efficiency",levs[-c(grep("Efficiency",levs))])
}
p <- ggplot(data=df_stack,aes(x=sample,y=values, group=ind))
p <- p+geom_line()+geom_point() + facet_wrap(~ind,scales="free")
# Compute sample mean and standard deviation in each group
if(mean_line){
ds <- df_stack %>% dplyr::group_by(ind) %>% dplyr::summarise(mean=mean(values),sd=stats::sd(values))
#ds <- dplyr::summarise_(dplyr::group_by(df_stack, ind), mean=~mean(values), sd=~sd(values))
p <- p + geom_hline(data = ds, aes(yintercept = mean),colour = mean_color)
}
#p
#p <- ggplot(data=df,aes(x=sample,y=efficiency)) #+ labs(colour = NULL)
#p <- ggplot(data=df,aes(x=Efficiency)) #+ labs(colour = NULL)
#p+geom_histogram()
#p <- p+geom_line()+geom_point()+ylab("Normalized Efficiency (%)")
p <- p+ylab("Value in %")+xlab("Simulation number of design sampled from defined windows")
#p <- p + stat_summary(fun.y = "mean", geom="hline")
return(p)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.