R/plot_efficiency_of_windows.R

In PopED: Population (and Individual) Optimal Experimental Design

#' 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)
}