From 815635d9b29a2687ed0bec184736d9a65b81ce85 Mon Sep 17 00:00:00 2001 From: Andrew Hooker Date: Tue, 27 May 2014 14:35:15 +0200 Subject: [PATCH] changes to make package build work --- DESCRIPTION | 6 +- NAMESPACE | 1 + NEWS | 1 + R/LEDoptim.R | 4 +- R/RS_opt_gen.R | 7 + R/blockfinal.R | 1 + R/blockheader.R | 2 + R/calc_ofv_and_fim.R | 2 + R/ed_laplace_ofv.R | 7 +- R/log_prior_pdf.R | 3 + R/poped_optimize.R | 1 + man/LEDoptim.Rd | 116 ++++++------- man/PopED-internal.Rd | 1 - man/RS_opt_gen.Rd | 16 ++ man/blockfinal.Rd | 10 +- man/blockfinal_2.Rd | 151 ----------------- man/blockheader.Rd | 8 + man/blockheader_2.Rd | 155 ------------------ man/blockopt_2.Rd | 70 -------- man/calc_ofv_and_fim.Rd | 4 + man/ed_laplace_ofv.Rd | 6 + man/evaluate.fim.Rd | 1 + man/get_rse.Rd | 1 + man/log_prior_pdf.Rd | 42 +++++ man/poped_optimize.Rd | 46 +++++- .../examples_fcn_doc/examples_blockfinal.R | 2 +- .../examples_ed_laplace_ofv.R | 1 + .../examples_poped_optimize.R | 6 +- 28 files changed, 215 insertions(+), 456 deletions(-) delete mode 100644 man/blockfinal_2.Rd delete mode 100644 man/blockheader_2.Rd delete mode 100644 man/blockopt_2.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 942bfcce..7a93bb31 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,12 +4,8 @@ Title: Population (and individual) optimal Experimental Design Version: 0.1.1 Date: 2014-05-27 Depends: ggplot2 -Imports: MASS, mvtnorm +Imports: MASS, mvtnorm, nlme Suggests: testthat, Hmisc -Author: Andrew C. Hooker, Sebastian Ueckert (MATLAB version), - Marco Foracchia (O-Matrix version) and Joakim Nyberg (MATLAB version) - with contributions from Eric Stroemberg (MATLAB version). -Maintainer: Andrew C. Hooker Authors@R: c(person("Andrew C.","Hooker", email="andrew.hooker@farmbio.uu.se", role=c("aut","cre","trl","cph")), diff --git a/NAMESPACE b/NAMESPACE index 2a45c263..3fc72633 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,5 @@ import(ggplot2) importFrom(MASS, write.matrix) importFrom(mvtnorm, rmvnorm) +importFrom(nlme, fdHess) exportPattern("^[[:alpha:]]+") diff --git a/NEWS b/NEWS index e65a593f..34c23032 100644 --- a/NEWS +++ b/NEWS @@ -7,6 +7,7 @@ PopED 0.1.1 * New functionality to compute the ED OFV using the Laplace approximation. This can be orders of magnitude faster than the standard MC integration approach. + See '?ed_laplace_ofv' and '?evaluate.e.ofv.fim' * Added a general function to compute the FIM and OFV(FIM) for all avaialbale methods in PopED. See '?calc_ofv_and_fim'. diff --git a/R/LEDoptim.R b/R/LEDoptim.R index 15227789..44a70703 100644 --- a/R/LEDoptim.R +++ b/R/LEDoptim.R @@ -15,9 +15,11 @@ #' @inheritParams RS_opt_gen #' @inheritParams create.poped.database #' @inheritParams Doptim +#' @inheritParams calc_ofv_and_fim #' @param fim_init The initial value of the FIM. If set to zero then it is computed. #' @param ofv_init The inital OFV. If set to zero then it is computed. #' @param trflag Should the optimization be output to the screen and to a file? +#' @param use_RS should the fucntion use a random search algorithm? #' #' @family Optimize #' @@ -203,7 +205,7 @@ LEDoptim <- function(poped.db, fprintf('Starting BGFS minimization with OFV of %g \n', best_dmf) returnArgs <- bfgsb_min('ed_laplace_ofv', list(model_switch,aa,axt,poped.db$groupsize,ni, - xtopt,xopt,aopt,bpop,d,poped.db$sigma,docc_full,poped.db, + xtopt,xopt,aopt,bpopdescr,ddescr,poped.db$sigma,poped.db$param.pt.val$docc,poped.db, return_gradient=T, optxt=optxt, opta=optx, x=x), x_k,lb,ub,options) diff --git a/R/RS_opt_gen.R b/R/RS_opt_gen.R index 30f9dad7..faa89712 100644 --- a/R/RS_opt_gen.R +++ b/R/RS_opt_gen.R @@ -12,12 +12,19 @@ #' @inheritParams evaluate.fim #' @inheritParams Doptim #' @inheritParams create.poped.database +#' @inheritParams calc_ofv_and_fim #' @param cfaxt First step factor for sample times #' @param opt_xt Should the sample times be optimized? #' @param opt_a Should the continuous design variables be optimized? #' @param opt_x Should the discrete design variables be optimized? #' @param approx_type Approximation method for model, 0=FO, 1=FOCE, 2=FOCEI, 3=FOI. #' @param iter The number of iterations entered into the \code{blockheader_2} function. +#' @param header_flag Should the header text be printed out? +#' @param footer_flag Should the footer text be printed out? +#' @param out_file Which file should the output be directed to? A string, a file handle using +#' \code{\link{file}} or \code{""} will output to the screen. +#' @param compute_inv should the inverse of the FIM be used to compute expected RSE values? Often not needed +#' except for diagnostic purposes. #' @param ... arguments passed to \code{\link{evaluate.fim}} and \code{\link{ofv_fim}}. #' #' diff --git a/R/blockfinal.R b/R/blockfinal.R index 71d3a0f6..1d32687e 100644 --- a/R/blockfinal.R +++ b/R/blockfinal.R @@ -8,6 +8,7 @@ #' @inheritParams create.poped.database #' @inheritParams blockexp #' @inheritParams blockheader +#' @inheritParams RS_opt_gen #' @param fmf_init Initial FIM. #' @param dmf_init Initial OFV. #' @param param_cvs_init The inital design parameter RSE values. diff --git a/R/blockheader.R b/R/blockheader.R index 8e4dffb0..b4760983 100644 --- a/R/blockheader.R +++ b/R/blockheader.R @@ -7,11 +7,13 @@ #' @inheritParams blockexp #' @inheritParams Doptim #' @inheritParams create.poped.database +#' @inheritParams RS_opt_gen #' @param name The name used for the output file. Combined with \code{name_header} and \code{iter}. #' If \code{""} then output is to the screen. #' @param iter The last number in the name printed to the output file, combined with \code{name}. #' @param name_header The initial portion of the file name. #' @param file_path The path to where the file should be created. +#' @param header_flag Should the header text be printed out? #' @param ... Additional arguments passed to further functions. #' #' @family Helper diff --git a/R/calc_ofv_and_fim.R b/R/calc_ofv_and_fim.R index 99e9b1ce..fc78dfe7 100644 --- a/R/calc_ofv_and_fim.R +++ b/R/calc_ofv_and_fim.R @@ -8,6 +8,8 @@ #' @inheritParams evaluate.fim #' @inheritParams Doptim #' @inheritParams create.poped.database +#' @param ofv The current ofv. If other than zero then this values is simply returned unchanged. +#' @param fim The current FIM. If other than zero then this values is simply returned unchanged. #' @param use_laplace Should the Laplace method be used in calculating the expectation of the OFV? #' @param laplace.fim Should an E(FIM) be calculated when computing the Laplace approximated E(OFV). Typically #' the FIM does not need to be computed and, if desired, this calculation diff --git a/R/ed_laplace_ofv.R b/R/ed_laplace_ofv.R index 78ff21b6..14f00173 100644 --- a/R/ed_laplace_ofv.R +++ b/R/ed_laplace_ofv.R @@ -20,7 +20,10 @@ #' @param optxt If sampling times are optimized #' @param opta If continuous design variables are optimized #' @param aopto the continuous design variables -#' +#' @param method If 0 then use an optimization routine translated from poped code written in MATLAB to +#' optimize the parameters in the Laplace approximation. If 1 then use \code{\link{optim}} to compute both +#' k and the hessian of k (see Dodds et al, JPP, 2005 for more information). If 2 then use \code{\link{fdHess}} +#' to compute the hessian. #' @param return_gradient Should the gradient be returned. #' #' @return The FIM and the hessian of the FIM. @@ -224,7 +227,7 @@ ed_laplace_ofv <- function(model_switch,groupsize,ni,xtopto,xopto,aopto, ## Different hessian and gradient calculation if(method==2){ - k_vals <- nlme::fdHess(output$par, + k_vals <- fdHess(output$par, function(x) calc_k(x,model_switch,groupsize,ni,xtopto,xopto, aopto,bpopdescr,ddescr,covd,sigma,docc,poped.db,Engine, return_gradient=F)[["k"]]) diff --git a/R/log_prior_pdf.R b/R/log_prior_pdf.R index 13dee422..ba91954d 100644 --- a/R/log_prior_pdf.R +++ b/R/log_prior_pdf.R @@ -4,7 +4,10 @@ #' @inheritParams evaluate.fim #' @inheritParams create.poped.database #' @inheritParams Doptim +#' @param alpha A parameter vector. +#' @param return_hessian Should the hessian be returned? #' +#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R #' @example tests/testthat/examples_fcn_doc/examples_log_prior_pdf.R log_prior_pdf <- function(alpha, bpopdescr, ddescr,return_gradient=F,return_hessian=F){ diff --git a/R/poped_optimize.R b/R/poped_optimize.R index 61d88049..cedb0c94 100644 --- a/R/poped_optimize.R +++ b/R/poped_optimize.R @@ -11,6 +11,7 @@ #' @inheritParams Doptim #' @inheritParams create.poped.database #' @inheritParams Dtrace +#' @inheritParams calc_ofv_and_fim #' @param ... arguments passed to other functions. See \code{\link{Doptim}}. #' #' diff --git a/man/LEDoptim.Rd b/man/LEDoptim.Rd index 65fe5133..7ba767ba 100644 --- a/man/LEDoptim.Rd +++ b/man/LEDoptim.Rd @@ -1,7 +1,7 @@ % Generated by roxygen2 (4.0.1): do not edit by hand \name{LEDoptim} \alias{LEDoptim} -\title{Laplace approximated ED-optimization function} +\title{Optimization function for D-family, E-family and Laplace approximated ED designs} \usage{ LEDoptim(poped.db, model_switch = NULL, ni = NULL, xt = NULL, x = NULL, a = NULL, bpopdescr = NULL, ddescr = NULL, maxxt = NULL, @@ -19,9 +19,7 @@ LEDoptim(poped.db, model_switch = NULL, ni = NULL, xt = NULL, x = NULL, \item{trflag}{Should the optimization be output to the screen and to a file?} -\item{ls_step_size}{Number of grid points in the line search} - -\item{iter_tot}{Number of iterations of full Random search and full Stochastic Gradient if line search is not used.} +\item{use_RS}{should the fucntion use a random search algorithm?} \item{opt_xt}{Should the sample times be optimized?} @@ -31,6 +29,13 @@ LEDoptim(poped.db, model_switch = NULL, ni = NULL, xt = NULL, x = NULL, \item{...}{arguments passed to \code{\link{evaluate.fim}} and \code{\link{ofv_fim}}.} +\item{header_flag}{Should the header text be printed out?} + +\item{footer_flag}{Should the footer text be printed out?} + +\item{out_file}{Which file should the output be directed to? A string, a file handle using +\code{\link{file}} or \code{""} will output to the screen.} + \item{model_switch}{Matrix defining which response a certain sampling time belongs to.} \item{ni}{Vector defining the number of samples for each group.} @@ -77,9 +82,15 @@ D-family design (1) or ED-familty design (0) (with or without parameter uncertai }} \item{ddescr}{Matrix defining the diagnonals of the IIV (same logic as for the \code{bpopdescr}).} + +\item{use_laplace}{Should the Laplace method be used in calculating the expectation of the OFV?} + +\item{laplace.fim}{Should an E(FIM) be calculated when computing the Laplace approximated E(OFV). Typically +the FIM does not need to be computed and, if desired, this calculation +is done usng the standard MC integration technique, so can be slow.} } \description{ -Optimize the Laplace approximated ED-objective function. +Optimize the objective fucntion for D-family, E-family and Laplace approximated ED designs. Right now there is only one optimization algorithm used in this function \enumerate{ @@ -114,11 +125,19 @@ sfg <- function(x,a,bpop,b,bocc){ return(parameters) } +# Adding 10\% log-normal Uncertainty to fixed effects (not Favail) +bpop_vals <- c(CL=0.15, V=8, KA=1.0, Favail=1) +bpop_vals_ed_ln <- cbind(ones(length(bpop_vals),1)*4, # log-normal distribution + bpop_vals, + ones(length(bpop_vals),1)*(bpop_vals*0.1)^2) # 10\% of bpop value +bpop_vals_ed_ln["Favail",] <- c(0,1,0) +bpop_vals_ed_ln + ## -- Define initial design and design space poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", fg_file="sfg", fError_file="feps.add.prop", - bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), + bpop=bpop_vals_ed_ln, notfixed_bpop=c(1,1,1,0), d=c(CL=0.07, V=0.02, KA=0.6), sigma=c(0.01,0.25), @@ -129,78 +148,45 @@ poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", a=70, mina=0, maxa=100) -# warfarin optimize model +# warfain ed model \dontrun{ - ############## - # typically one will use poped_optimize - # This then calls Doptim for continuous optimization problems - ############## + LEDoptim(poped.db) + LEDoptim(poped.db,opt_xt=T,rsit=10) - # RS+SG+LS optimization of sample times - # optimization with just a few iterations - # only to check that things are working - output <- poped_optimize(poped.db,opt_xt=T, - rsit=5,sgit=5,ls_step_size=5) + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE) + + LEDoptim(poped.db,opt_xt=T,rsit=10,laplace.fim=TRUE) - # RS+SG+LS optimization of sample times - # (longer run time than above but more likely to reach a maximum) - output <- poped_optimize(poped.db,opt_xt=T) - get_rse(output$fmf,output$poped.db) - plot_model_prediction(output$poped.db) + LEDoptim(poped.db,opt_xt=T,rsit=10,use_laplace=FALSE) + ## testing header and footer + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + out_file="foobar.txt") - # Random search (just a few samples here) - rs.output <- poped_optimize(poped.db,opt_xt=1,opt_a=1,rsit=20, - bUseRandomSearch= 1, - bUseStochasticGradient = 0, - bUseBFGSMinimizer = 0, - bUseLineSearch = 0) + ff <- LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + trflag=FALSE) - # line search, DOSE and sample time optimization - ls.output <- poped_optimize(poped.db,opt_xt=1,opt_a=1, - bUseRandomSearch= 0, - bUseStochasticGradient = 0, - bUseBFGSMinimizer = 0, - bUseLineSearch = 1, - ls_step_size=10) + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + header_flag=FALSE) - # Stochastic gradient search, DOSE and sample time optimization - sg.output <- poped_optimize(poped.db,opt_xt=1,opt_a=1, - bUseRandomSearch= 0, - bUseStochasticGradient = 1, - bUseBFGSMinimizer = 0, - bUseLineSearch = 0, - sgit=20) + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + out_file="") - # BFGS search, DOSE and sample time optimization - bfgs.output <- poped_optimize(poped.db,opt_xt=1,opt_a=1, - bUseRandomSearch= 0, - bUseStochasticGradient = 0, - bUseBFGSMinimizer = 1, - bUseLineSearch = 0) - - ############## - # If you really want to you can use Doptim dirtectly - ############## - dsl <- downsizing_general_design(poped.db) - poped.db$optsw[2] <- 1 # sample time optimization - output <- Doptim(poped.db,dsl$ni, dsl$xt, dsl$model_switch, dsl$x, dsl$a, - dsl$bpop, dsl$d, dsl$maxxt, dsl$minxt,dsl$maxa,dsl$mina) + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + footer_flag=FALSE) -} - + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + footer_flag=FALSE, header_flag=FALSE) + + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + footer_flag=FALSE, header_flag=FALSE,out_file="foobar.txt") + + LEDoptim(poped.db,opt_xt=T,rsit=10,d_switch=TRUE, + footer_flag=FALSE, header_flag=FALSE,out_file="") -} -\references{ -\enumerate{ -\item M. Foracchia, A.C. Hooker, P. Vicini and A. Ruggeri, "PopED, a software for optimal -experimental design in population kinetics", Computer Methods and Programs in Biomedicine, 74, 2004. -\item J. Nyberg, S. Ueckert, E.A. Stroemberg, S. Hennig, M.O. Karlsson and A.C. Hooker, "PopED: An extended, -parallelized, nonlinear mixed effects models optimal design tool", -Computer Methods and Programs in Biomedicine, 108, 2012. } } \seealso{ diff --git a/man/PopED-internal.Rd b/man/PopED-internal.Rd index 182c2135..a1a69592 100644 --- a/man/PopED-internal.Rd +++ b/man/PopED-internal.Rd @@ -46,7 +46,6 @@ \alias{numel} \alias{new_ferror_file} \alias{min_function} -\alias{log_prior_pdf} \alias{linspace} \alias{line_search_uc} \alias{line_search} diff --git a/man/RS_opt_gen.Rd b/man/RS_opt_gen.Rd index 08de98f7..11539fc0 100644 --- a/man/RS_opt_gen.Rd +++ b/man/RS_opt_gen.Rd @@ -28,6 +28,16 @@ RS_opt_gen(poped.db, ni = NULL, xt = NULL, model_switch = NULL, \item{iter}{The number of iterations entered into the \code{blockheader_2} function.} +\item{header_flag}{Should the header text be printed out?} + +\item{footer_flag}{Should the footer text be printed out?} + +\item{out_file}{Which file should the output be directed to? A string, a file handle using +\code{\link{file}} or \code{""} will output to the screen.} + +\item{compute_inv}{should the inverse of the FIM be used to compute expected RSE values? Often not needed +except for diagnostic purposes.} + \item{...}{arguments passed to \code{\link{evaluate.fim}} and \code{\link{ofv_fim}}.} \item{poped.db}{A PopED database.} @@ -94,6 +104,12 @@ all a values are given the same max value} \item{d_switch}{\itemize{ \item \bold{******START OF CRITERION SPECIFICATION OPTIONS**********}} D-family design (1) or ED-familty design (0) (with or without parameter uncertainty)} + +\item{use_laplace}{Should the Laplace method be used in calculating the expectation of the OFV?} + +\item{laplace.fim}{Should an E(FIM) be calculated when computing the Laplace approximated E(OFV). Typically +the FIM does not need to be computed and, if desired, this calculation +is done usng the standard MC integration technique, so can be slow.} } \description{ Optimize the objective function using an adaptive random search algorithm. diff --git a/man/blockfinal.Rd b/man/blockfinal.Rd index b3f817c6..f760e6be 100644 --- a/man/blockfinal.Rd +++ b/man/blockfinal.Rd @@ -62,6 +62,14 @@ can also just supply the parameter values as a \code{c()}.} can also just supply the diagnonal parameter values (variances) as a \code{c()}.} \item{fn}{The file handle to write to.} + +\item{compute_inv}{should the inverse of the FIM be used to compute expected RSE values? Often not needed +except for diagnostic purposes.} + +\item{out_file}{Which file should the output be directed to? A string, a file handle using +\code{\link{file}} or \code{""} will output to the screen.} + +\item{footer_flag}{Should the footer text be printed out?} } \description{ Create some output to the screen and a text file that summarizes the problem you solved. @@ -109,7 +117,7 @@ poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", FIM <- evaluate.fim(poped.db) - <- det(FIM) +dmf <- det(FIM) blockfinal(fn="",fmf=FIM, diff --git a/man/blockfinal_2.Rd b/man/blockfinal_2.Rd deleted file mode 100644 index ac093269..00000000 --- a/man/blockfinal_2.Rd +++ /dev/null @@ -1,151 +0,0 @@ -\name{blockfinal_2} -\alias{blockfinal_2} -\title{Result function for optimization routines} -\usage{ -blockfinal_2(fn, fmf, dmf, groupsize, ni, xt, x, a, model_switch, bpop, d, docc, - sigma, m, poped.db, opt_xt = poped.db$optsw[2], opt_a = poped.db$optsw[4], - opt_x = poped.db$optsw[4], fmf_init = NULL, dmf_init = NULL, - param_cvs_init = NULL) -} -\arguments{ - \item{fmf_init}{Initial FIM.} - - \item{dmf_init}{Initial OFV.} - - \item{param_cvs_init}{The inital design parameter RSE - values.} - - \item{fmf}{The initial value of the FIM. If set to zero - then it is computed.} - - \item{dmf}{The inital OFV. If set to zero then it is - computed.} - - \item{ni}{A vector of the number of samples in each - group.} - - \item{xt}{A matrix of sample times. Each row is a vector - of sample times for a group.} - - \item{x}{A matrix for the discrete design variables. - Each row is a group.} - - \item{a}{A matrix of covariates. Each row is a group.} - - \item{model_switch}{A matrix that is the same size as xt, - specifying which model each sample belongs to.} - - \item{poped.db}{A PopED database.} - - \item{opt_xt}{Should the sample times be optimized?} - - \item{opt_a}{Should the continuous design variables be - optimized?} - - \item{opt_x}{Should the discrete design variables be - optimized?} - - \item{groupsize}{A vector of the numer of individuals in - each group.} - - \item{bpop}{Matrix defining the fixed effects, per row - (row number = parameter_number) we should have: \itemize{ - \item column 1 the type of the distribution for E-family - designs (0 = Fixed, 1 = Normal, 2 = Uniform, 3 = User - Defined Distribution, 4 = lognormal and 5 = truncated - normal) \item column 2 defines the mean. \item column 3 - defines the variance of the distribution (or length of - uniform distribution). } Can also just supply the - parameter values as a vector \code{c()}} - - \item{d}{Matrix defining the diagnonals of the IIV (same - logic as for the fixed efects). can also just supply the - parameter values as a \code{c()}.} - - \item{docc}{Matrix defining the IOV, the IOV variances - and the IOV distribution} - - \item{sigma}{Matrix defining the variances can - covariances of the residual variability terms of the - model. can also just supply the diagnonal parameter - values (variances) as a \code{c()}.} - - \item{m}{Number of groups/individuals} - - \item{fn}{The file handle to write to.} -} -\description{ -Create some output to the screen and a text file that -summarizes the problem you solved. -} -\examples{ -## Warfarin example from software comparison in: -## Nyberg et al., "Methods and software tools for design evaluation -## for population pharmacokinetics-pharmacodynamics studies", -## Br. J. Clin. Pharm., 2014. - -## Optimization using an additive + proportional reidual error to -## avoid sample times at very low concentrations (time 0 or very late samoples). -library(PopED) - -## find the parameters that are needed to define from the structural model -ff.PK.1.comp.oral.sd.CL - -## -- parameter definition function -## -- names match parameters in function ff -sfg <- function(x,a,bpop,b,bocc){ - parameters=c(CL=bpop[1]*exp(b[1]), - V=bpop[2]*exp(b[2]), - KA=bpop[3]*exp(b[3]), - Favail=bpop[4], - DOSE=a[1]) - return(parameters) -} - -## -- Define initial design and design space -poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", - fg_file="sfg", - fError_file="feps.add.prop", - bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), - notfixed_bpop=c(1,1,1,0), - d=c(CL=0.07, V=0.02, KA=0.6), - sigma=c(0.01,0.25), - groupsize=32, - xt=c( 0.5,1,2,6,24,36,72,120), - minxt=0, - maxxt=120, - a=70, - mina=0, - maxa=100) -# warfarin optimization model - - -FIM <- evaluate.fim(poped.db) -dmf <- det(FIM) - - -blockfinal_2(fn="",fmf=FIM, - dmf=dmf, - groupsize=poped.db$groupsize, - ni=poped.db$gni, - xt=poped.db$gxt, - x=poped.db$gx,a=poped.db$ga, - model_switch=poped.db$global_model_switch, - poped.db$param.pt.val$bpop, - poped.db$param.pt.val$d, - poped.db$docc, - poped.db$param.pt.val$sigma, - poped.db$m, - poped.db, - opt_xt=TRUE, - fmf_init=FIM, - dmf_init=dmf, - param_cvs_init=rbind(get_rse(FIM,poped.db,use_percent=FALSE))) - - -} -\seealso{ -Other Helper: \code{\link{blockexp}}; -\code{\link{blockheader_2}}; \code{\link{blockopt_2}} -} - diff --git a/man/blockheader.Rd b/man/blockheader.Rd index 1ed634e0..16739abe 100644 --- a/man/blockheader.Rd +++ b/man/blockheader.Rd @@ -22,6 +22,8 @@ If \code{""} then output is to the screen.} \item{file_path}{The path to where the file should be created.} +\item{header_flag}{Should the header text be printed out?} + \item{...}{Additional arguments passed to further functions.} \item{poped.db}{A PopED database.} @@ -60,6 +62,12 @@ can also just supply the parameter values as a \code{c()}.} \item{sigma}{Matrix defining the variances can covariances of the residual variability terms of the model. can also just supply the diagnonal parameter values (variances) as a \code{c()}.} + +\item{out_file}{Which file should the output be directed to? A string, a file handle using +\code{\link{file}} or \code{""} will output to the screen.} + +\item{compute_inv}{should the inverse of the FIM be used to compute expected RSE values? Often not needed +except for diagnostic purposes.} } \value{ fn A file handle (or \code{''} if \code{name=''}) diff --git a/man/blockheader_2.Rd b/man/blockheader_2.Rd deleted file mode 100644 index 26e1e27a..00000000 --- a/man/blockheader_2.Rd +++ /dev/null @@ -1,155 +0,0 @@ -\name{blockheader_2} -\alias{blockheader_2} -\title{Header function for optimization routines} -\usage{ -blockheader_2(name, iter, poped.db, e_flag = FALSE, - opt_xt = poped.db$optsw[2], opt_a = poped.db$optsw[4], - opt_x = poped.db$optsw[4], opt_samps = poped.db$optsw[1], - opt_inds = poped.db$optsw[5], fmf = 0, dmf = 0, bpop = NULL, - d = NULL, docc = NULL, sigma = NULL, - name_header = poped.db$strOutputFileName, - file_path = poped.db$strOutputFilePath, ...) -} -\arguments{ - \item{name}{The name used for the output file. Combined - with \code{name_header} and \code{iter}. If \code{""} - then output is to the screen.} - - \item{iter}{The last number in the name printed to the - output file, combined with \code{name}.} - - \item{name_header}{The initial portion of the file name.} - - \item{file_path}{The path to where the file should be - created.} - - \item{...}{Additional arguments passed to further - functions.} - - \item{poped.db}{A PopED database.} - - \item{opt_xt}{Should the sample times be optimized?} - - \item{opt_a}{Should the continuous design variables be - optimized?} - - \item{opt_x}{Should the discrete design variables be - optimized?} - - \item{fmf}{The initial value of the FIM. If set to zero - then it is computed.} - - \item{dmf}{The inital OFV. If set to zero then it is - computed.} - - \item{e_flag}{Shuould output be with uncertainty around - parameters?} - - \item{opt_samps}{Are the nuber of sample times per group - being optimized?} - - \item{opt_inds}{Are the nuber of individuals per group - being optimized?} - - \item{bpop}{Matrix defining the fixed effects, per row - (row number = parameter_number) we should have: \itemize{ - \item column 1 the type of the distribution for E-family - designs (0 = Fixed, 1 = Normal, 2 = Uniform, 3 = User - Defined Distribution, 4 = lognormal and 5 = truncated - normal) \item column 2 defines the mean. \item column 3 - defines the variance of the distribution (or length of - uniform distribution). } Can also just supply the - parameter values as a vector \code{c()}} - - \item{d}{Matrix defining the diagnonals of the IIV (same - logic as for the fixed efects). can also just supply the - parameter values as a \code{c()}.} - - \item{docc}{Matrix defining the IOV, the IOV variances - and the IOV distribution} - - \item{sigma}{Matrix defining the variances can - covariances of the residual variability terms of the - model. can also just supply the diagnonal parameter - values (variances) as a \code{c()}.} -} -\value{ -fn A file handle (or \code{''} if \code{name=''}) -} -\description{ -Create some output to the screen and a text file that -summarizes the problem you are tying to solve. -} -\examples{ -## Warfarin example from software comparison in: -## Nyberg et al., "Methods and software tools for design evaluation -## for population pharmacokinetics-pharmacodynamics studies", -## Br. J. Clin. Pharm., 2014. - -## Optimization using an additive + proportional reidual error to -## avoid sample times at very low concentrations (time 0 or very late samoples). -library(PopED) - -## find the parameters that are needed to define from the structural model -ff.PK.1.comp.oral.sd.CL - -## -- parameter definition function -## -- names match parameters in function ff -sfg <- function(x,a,bpop,b,bocc){ - parameters=c(CL=bpop[1]*exp(b[1]), - V=bpop[2]*exp(b[2]), - KA=bpop[3]*exp(b[3]), - Favail=bpop[4], - DOSE=a[1]) - return(parameters) -} - -## -- Define initial design and design space -poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", - fg_file="sfg", - fError_file="feps.add.prop", - bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), - notfixed_bpop=c(1,1,1,0), - d=c(CL=0.07, V=0.02, KA=0.6), - sigma=c(0.01,0.25), - groupsize=32, - xt=c( 0.5,1,2,6,24,36,72,120), - minxt=0, - maxxt=120, - a=70, - mina=0, - maxa=100) -# warfarin optimization model - - -FIM <- evaluate.fim(poped.db) -dmf <- det(FIM) - -blockheader_2(name="",iter=1,poped.db) - - -blockheader_2(name='', - iter=1, - poped.db, - e_flag=FALSE, - opt_xt=TRUE, - opt_a=TRUE,opt_x=poped.db$optsw[4], - opt_samps=poped.db$optsw[1],opt_inds=poped.db$optsw[5], - fmf=FIM,dmf=dmf, - bpop=poped.db$param.pt.val$bpop, - d=poped.db$param.pt.val$d, - docc=poped.db$docc,sigma=poped.db$param.pt.val$sigma) - - - - - - - - -} -\seealso{ -Other Helper: \code{\link{blockexp}}; -\code{\link{blockfinal_2}}; \code{\link{blockopt_2}} -} - diff --git a/man/blockopt_2.Rd b/man/blockopt_2.Rd deleted file mode 100644 index ebc97b1a..00000000 --- a/man/blockopt_2.Rd +++ /dev/null @@ -1,70 +0,0 @@ -\name{blockopt_2} -\alias{blockopt_2} -\title{Summarize your optimization settings for optimization routines} -\usage{ -blockopt_2(fn, poped.db, opt_method = "") -} -\arguments{ - \item{opt_method}{If "RS" (random search), "SG" - (stochastic gradient) or "DO" (discrete optimization) - then specifc output is produced.} - - \item{poped.db}{A PopED database.} - - \item{fn}{The file handle to write to.} -} -\description{ -Create some output to the screen and a text file that -summarizes the optimization settings you will use to -optimize. -} -\examples{ -## Warfarin example from software comparison in: -## Nyberg et al., "Methods and software tools for design evaluation -## for population pharmacokinetics-pharmacodynamics studies", -## Br. J. Clin. Pharm., 2014. - -## Optimization using an additive + proportional reidual error to -## avoid sample times at very low concentrations (time 0 or very late samoples). -library(PopED) - -## find the parameters that are needed to define from the structural model -ff.PK.1.comp.oral.sd.CL - -## -- parameter definition function -## -- names match parameters in function ff -sfg <- function(x,a,bpop,b,bocc){ - parameters=c(CL=bpop[1]*exp(b[1]), - V=bpop[2]*exp(b[2]), - KA=bpop[3]*exp(b[3]), - Favail=bpop[4], - DOSE=a[1]) - return(parameters) -} - -## -- Define initial design and design space -poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", - fg_file="sfg", - fError_file="feps.add.prop", - bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), - notfixed_bpop=c(1,1,1,0), - d=c(CL=0.07, V=0.02, KA=0.6), - sigma=c(0.01,0.25), - groupsize=32, - xt=c( 0.5,1,2,6,24,36,72,120), - minxt=0, - maxxt=120, - a=70, - mina=0, - maxa=100) -# warfarin optimization model - -blockopt_2(fn="",poped.db,opt_method="SG") -blockopt_2(fn="",poped.db,opt_method="RS") -blockopt_2(fn="",poped.db,opt_method="DO") -} -\seealso{ -Other Helper: \code{\link{blockexp}}; -\code{\link{blockfinal_2}}; \code{\link{blockheader_2}} -} - diff --git a/man/calc_ofv_and_fim.Rd b/man/calc_ofv_and_fim.Rd index f3e6be28..61d041ab 100644 --- a/man/calc_ofv_and_fim.Rd +++ b/man/calc_ofv_and_fim.Rd @@ -13,6 +13,10 @@ calc_ofv_and_fim(poped.db, ofv = 0, fim = 0, d_switch = poped.db$d_switch, use_laplace = poped.db$iEDCalculationType, laplace.fim = FALSE, ...) } \arguments{ +\item{ofv}{The current ofv. If other than zero then this values is simply returned unchanged.} + +\item{fim}{The current FIM. If other than zero then this values is simply returned unchanged.} + \item{use_laplace}{Should the Laplace method be used in calculating the expectation of the OFV?} \item{laplace.fim}{Should an E(FIM) be calculated when computing the Laplace approximated E(OFV). Typically diff --git a/man/ed_laplace_ofv.Rd b/man/ed_laplace_ofv.Rd index 3d934085..a89a29ef 100644 --- a/man/ed_laplace_ofv.Rd +++ b/man/ed_laplace_ofv.Rd @@ -21,6 +21,11 @@ ed_laplace_ofv(model_switch, groupsize, ni, xtopto, xopto, aopto, bpopdescr, \item{aopto}{the continuous design variables} +\item{method}{If 0 then use an optimization routine translated from poped code written in MATLAB to +optimize the parameters in the Laplace approximation. If 1 then use \code{\link{optim}} to compute both +k and the hessian of k (see Dodds et al, JPP, 2005 for more information). If 2 then use \code{\link{fdHess}} +to compute the hessian.} + \item{return_gradient}{Should the gradient be returned.} \item{model_switch}{A matrix that is the same size as xt, specifying which model each sample belongs to.} @@ -90,6 +95,7 @@ sfg <- function(x,a,bpop,b,bocc){ ###################### # Normal distribution ###################### +bpop_vals <- c(CL=0.15, V=8, KA=1.0, Favail=1) bpop_vals_ed_n <- cbind(ones(length(bpop_vals),1)*1, # normal distribution bpop_vals, ones(length(bpop_vals),1)*(bpop_vals*0.1)^2) # 10\% of bpop value diff --git a/man/evaluate.fim.Rd b/man/evaluate.fim.Rd index 15df35ca..661a18cb 100644 --- a/man/evaluate.fim.Rd +++ b/man/evaluate.fim.Rd @@ -136,6 +136,7 @@ FIM.7 <- evaluate.fim(poped.db,fim.calc.type=7) FIM.7 det(FIM.7) get_rse(FIM.7,poped.db,fim.calc.type=7) + } \seealso{ Other FIM: \code{\link{LinMatrixH}}; diff --git a/man/get_rse.Rd b/man/get_rse.Rd index 3c905b13..345e9051 100644 --- a/man/get_rse.Rd +++ b/man/get_rse.Rd @@ -125,6 +125,7 @@ FIM.7 <- evaluate.fim(poped.db,fim.calc.type=7) FIM.7 det(FIM.7) get_rse(FIM.7,poped.db,fim.calc.type=7) + } \seealso{ Other evaluate_design: \code{\link{evaluate.fim}}; diff --git a/man/log_prior_pdf.Rd b/man/log_prior_pdf.Rd index 70fa3d31..a814c5f6 100644 --- a/man/log_prior_pdf.Rd +++ b/man/log_prior_pdf.Rd @@ -7,6 +7,10 @@ log_prior_pdf(alpha, bpopdescr, ddescr, return_gradient = F, return_hessian = F) } \arguments{ +\item{alpha}{A parameter vector.} + +\item{return_hessian}{Should the hessian be returned?} + \item{bpopdescr}{Matrix defining the fixed effects, per row (row number = parameter_number) we should have: \itemize{ \item column 1 the type of the distribution for E-family designs (0 = Fixed, 1 = Normal, 2 = Uniform, @@ -23,6 +27,44 @@ log_prior_pdf(alpha, bpopdescr, ddescr, return_gradient = F, Compute the natural log of the PDF for the parameters in an E-family design } \examples{ +## Warfarin example from software comparison in: +## Nyberg et al., "Methods and software tools for design evaluation +## for population pharmacokinetics-pharmacodynamics studies", +## Br. J. Clin. Pharm., 2014. + +## Optimization using an additive + proportional reidual error to +## avoid sample times at very low concentrations (time 0 or very late samoples). +library(PopED) + +## find the parameters that are needed to define from the structural model +ff.PK.1.comp.oral.sd.CL + +## -- parameter definition function +## -- names match parameters in function ff +sfg <- function(x,a,bpop,b,bocc){ + parameters=c(CL=bpop[1]*exp(b[1]), + V=bpop[2]*exp(b[2]), + KA=bpop[3]*exp(b[3]), + Favail=bpop[4], + DOSE=a[1]) + return(parameters) +} + +## -- Define initial design and design space +poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", + fg_file="sfg", + fError_file="feps.add.prop", + bpop=c(CL=0.15, V=8, KA=1.0, Favail=1), + notfixed_bpop=c(1,1,1,0), + d=c(CL=0.07, V=0.02, KA=0.6), + sigma=c(0.01,0.25), + groupsize=32, + xt=c( 0.5,1,2,6,24,36,72,120), + minxt=0, + maxxt=120, + a=70, + mina=0, + maxa=100) # warfarin optimize model # Adding 40\% Uncertainty to fixed effects log-normal (not Favail) diff --git a/man/poped_optimize.Rd b/man/poped_optimize.Rd index dacf158f..b65a46de 100644 --- a/man/poped_optimize.Rd +++ b/man/poped_optimize.Rd @@ -16,7 +16,7 @@ poped_optimize(poped.db, ni = NULL, xt = NULL, model_switch = NULL, approx_type = poped.db$iApproximationMethod, bUseExchangeAlgorithm = poped.db$bUseExchangeAlgorithm, iter = 1, d_switch = poped.db$d_switch, ED_samp_size = poped.db$ED_samp_size, - bLHS = poped.db$bLHS, ...) + bLHS = poped.db$bLHS, use_laplace = poped.db$iEDCalculationType, ...) } \arguments{ \item{...}{arguments passed to other functions. See \code{\link{Doptim}}.} @@ -120,6 +120,8 @@ D-family design (1) or ED-familty design (0) (with or without parameter uncertai \item{opt_samps}{Are the nuber of sample times per group being optimized?} \item{opt_inds}{Are the nuber of individuals per group being optimized?} + +\item{use_laplace}{Should the Laplace method be used in calculating the expectation of the OFV?} } \description{ Optimize the objective function. @@ -173,7 +175,7 @@ poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", \dontrun{ ############## - # Optimization + # D-family Optimization ############## # below are a number of ways to optimize the problem @@ -233,6 +235,46 @@ poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", bUseBFGSMinimizer = 1, bUseLineSearch = 0) + + ############## + # E-family Optimization + ############## + + # Adding 10\% log-normal Uncertainty to fixed effects (not Favail) + bpop_vals <- c(CL=0.15, V=8, KA=1.0, Favail=1) + bpop_vals_ed_ln <- cbind(ones(length(bpop_vals),1)*4, # log-normal distribution + bpop_vals, + ones(length(bpop_vals),1)*(bpop_vals*0.1)^2) # 10\% of bpop value + bpop_vals_ed_ln["Favail",] <- c(0,1,0) + bpop_vals_ed_ln + + ## -- Define initial design and design space + poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.sd.CL", + fg_file="sfg", + fError_file="feps.add.prop", + bpop=bpop_vals_ed_ln, + notfixed_bpop=c(1,1,1,0), + d=c(CL=0.07, V=0.02, KA=0.6), + sigma=c(0.01,0.25), + groupsize=32, + xt=c( 0.5,1,2,6,24,36,72,120), + minxt=0, + maxxt=120, + a=70, + mina=0, + maxa=100) + + # ED optimization using Random search (just a few samples here) + output <- poped_optimize(poped.db,opt_xt=1,opt_a=1,rsit=10,d_switch=0) + get_rse(output$fmf,output$poped.db) + + # ED with laplace approximation, + # optimization using Random search (just a few samples here) + output <- poped_optimize(poped.db,opt_xt=1,opt_a=1,rsit=10, + d_switch=0,use_laplace=TRUE,laplace.fim=TRUE) + get_rse(output$fmf,output$poped.db) + + } } \references{ diff --git a/tests/testthat/examples_fcn_doc/examples_blockfinal.R b/tests/testthat/examples_fcn_doc/examples_blockfinal.R index 6c79736b..c3303b60 100644 --- a/tests/testthat/examples_fcn_doc/examples_blockfinal.R +++ b/tests/testthat/examples_fcn_doc/examples_blockfinal.R @@ -2,7 +2,7 @@ FIM <- evaluate.fim(poped.db) - <- det(FIM) +dmf <- det(FIM) blockfinal(fn="",fmf=FIM, diff --git a/tests/testthat/examples_fcn_doc/examples_ed_laplace_ofv.R b/tests/testthat/examples_fcn_doc/examples_ed_laplace_ofv.R index 0d5b544a..8530b175 100644 --- a/tests/testthat/examples_fcn_doc/examples_ed_laplace_ofv.R +++ b/tests/testthat/examples_fcn_doc/examples_ed_laplace_ofv.R @@ -24,6 +24,7 @@ sfg <- function(x,a,bpop,b,bocc){ ###################### # Normal distribution ###################### +bpop_vals <- c(CL=0.15, V=8, KA=1.0, Favail=1) bpop_vals_ed_n <- cbind(ones(length(bpop_vals),1)*1, # normal distribution bpop_vals, ones(length(bpop_vals),1)*(bpop_vals*0.1)^2) # 10% of bpop value diff --git a/tests/testthat/examples_fcn_doc/examples_poped_optimize.R b/tests/testthat/examples_fcn_doc/examples_poped_optimize.R index 1727cbad..6a24fff0 100644 --- a/tests/testthat/examples_fcn_doc/examples_poped_optimize.R +++ b/tests/testthat/examples_fcn_doc/examples_poped_optimize.R @@ -96,8 +96,10 @@ output <- poped_optimize(poped.db,opt_xt=1,opt_a=1,rsit=10,d_switch=0) get_rse(output$fmf,output$poped.db) - # ED with laplace approximation, optimization using Random search (just a few samples here) - output <- poped_optimize(poped.db,opt_xt=1,opt_a=1,rsit=10,d_switch=0,use_laplace=TRUE,laplace.fim=TRUE) + # ED with laplace approximation, + # optimization using Random search (just a few samples here) + output <- poped_optimize(poped.db,opt_xt=1,opt_a=1,rsit=10, + d_switch=0,use_laplace=TRUE,laplace.fim=TRUE) get_rse(output$fmf,output$poped.db)