R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("nlmeODE-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('nlmeODE') Loading required package: odesolve Loading required package: nlme Loading required package: lattice > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "PKPDmodels" > > ### * PKPDmodels > > flush(stderr()); flush(stdout()) > > ### Name: PKPDmodels > ### Title: Library of PK/PD models > ### Aliases: PKPDmodels > ### Keywords: models > > ### ** Examples > > > ######################################## > ### Pharmacokinetics of Theophylline ### > ######################################## > > data(Theoph) > > TheophODE <- Theoph > TheophODE$Dose[TheophODE$Time!=0] <- 0 > TheophODE$Cmt <- rep(1,dim(TheophODE)[1]) > > OneComp <- list(DiffEq=list( + dy1dt = ~ -ka*y1 , + dy2dt = ~ ka*y1-ke*y2), + ObsEq=list( + c1 = ~ 0, + c2 = ~ y2/CL*ke), + Parms=c("ka","ke","CL"), + States=c("y1","y2"), + Init=list(0,0)) > > TheophModel <- nlmeODE(OneComp,TheophODE) > > #Remove '#' below to run the estimation > > #Theoph.nlme <- nlme(conc ~ TheophModel(ka,ke,CL,Time,Subject), > # data = TheophODE, fixed=ka+ke+CL~1, random = pdDiag(ka+CL~1), > # start=c(ka=0.5,ke=-2.5,CL=-3.2), > # control=list(returnObject=TRUE,msVerbose=TRUE,tolerance=1e-1,pnlsTol=1e-1,msTol=1e-1), > # verbose=TRUE) > > #plot(augPred(Theoph.nlme,level=0:1)) > > ######################################### > ### Pharmacokinetics of Indomethacine ### > ######################################### > > data(Indometh) > > TwoComp <- list(DiffEq=list( + dy1dt = ~ -(k12+k10)*y1+k21*y2 , + dy2dt = ~ -k21*y2 + k12*y1), + ObsEq=list( + c1 = ~ y1, + c2 = ~ 0), + States=c("y1","y2"), + Parms=c("k12","k21","k10","start"), + Init=list("start",0)) > > IndomethModel <- nlmeODE(TwoComp,Indometh) > > #Remove '#' below to run the estimation > > #Indometh.nlme <- nlme(conc ~ IndomethModel(k12,k21,k10,start,time,Subject), > # data = Indometh, fixed=k12+k21+k10+start~1, random = pdDiag(start+k12+k10~1), > # start=c(k12=-0.05,k21=-0.15,k10=-0.10,start=0.70), > # control=list(msVerbose=TRUE,tolerance=1e-1,pnlsTol=1e-1,msTol=1e-1), > # verbose=TRUE) > > #plot(augPred(Indometh.nlme,level=0:1)) > > ################################################################# > ### Absorption model with estimation of time/rate of infusion ### > ################################################################# > > OneCompAbs <- list(DiffEq=list( + dA1dt = ~ -ka*A1, + dA2dt = ~ ka*A1 - CL/V1*A2), + ObsEq=list( + SC= ~0, + C = ~ A2/V1), + States=c("A1","A2"), + Parms=c("ka","CL","V1","F1"), + Init=list(0,0)) > > ID <- rep(seq(1:18),each=11) > Time <- rep(seq(0,100,by=10),18) > Dose <- c(rep(c(100,0,0,100,0,0,0,0,0,0,0),6),rep(c(100,0,0,0,0,0,0,100,0,0,0),6),rep(c(100,0,0,0,0,0,0,0,0,0,0),6)) > Rate <- c(rep(rep(0,11),6),rep(c(5,rep(0,10)),6),rep(rep(0,11),6)) > Cmt <- c(rep(1,6*11),rep(c(2,0,0,0,0,0,0,1,0,0,0),6),rep(2,6*11)) > Conc <- rep(0,18*11) > > Data <- as.data.frame(list(ID=ID,Time=Time,Dose=Dose,Rate=Rate,Cmt=Cmt,Conc=Conc)) > > SimData <- groupedData( Conc ~ Time | ID, + data = Data, + labels = list( x = "Time", y = "Concentration")) > > OneCompAbsModel <- nlmeODE(OneCompAbs,SimData) > > kaSim <- rep(log(rep(0.05,18))+0.3*rnorm(18),each=11) > CLSim <- rep(log(rep(0.5,18))+0.2*rnorm(18),each=11) > V1Sim <- rep(log(rep(10,18))+0.1*rnorm(18),each=11) > F1Sim <- rep(log(0.8),18*11) > > SimData$Sim <- OneCompAbsModel(kaSim,CLSim,V1Sim,F1Sim,SimData$Time,SimData$ID) > > SimData$Conc <- SimData$Sim + 0.3*rnorm(dim(SimData)[1]) > > Data <- groupedData( Conc ~ Time | ID, + data = SimData, + labels = list( x = "Time", y = "Concentration")) > > plot(Data,aspect=1/1) > > #Estimation of model parameters > OneCompAbsModel <- nlmeODE(OneCompAbs,Data) > > #Remove '#' below to run the estimation > > #fit1 <- nlme(Conc ~ OneCompAbsModel(ka,CL,V1,F1,Time,ID), > # data = Data, fixed=ka+CL+V1+F1~1, random = pdDiag(ka+CL+V1~1), > # start=c(ka=log(0.05),CL=log(0.5),V1=log(10.0),F1=log(0.8)), > # control=list(msVerbose=TRUE,tolerance=1e-3,pnlsTol=1e-1,msTol=1e-3), > # verbose=TRUE) > #plot(augPred(fit1,level=0:1,length.out=300),aspect=1/1) > > #Estimation of rate of infusion > Data$Rate[Data$Rate==5] <- -1 > > OneCompAbs <- list(DiffEq=list( + dA1dt = ~ -ka*A1, + dA2dt = ~ ka*A1 - CL/V1*A2), + ObsEq=list( + SC= ~0, + C = ~ A2/V1), + States=c("A1","A2"), + Parms=c("ka","CL","V1","F1","Rate"), + Init=list(0,0)) > > OneCompAbsModel <- nlmeODE(OneCompAbs,Data) > > #Remove '#' below to run the estimation > > #fit2 <- nlme(Conc ~ OneCompAbsModel(ka,CL,V1,F1,Rate,Time,ID), > # data = Data, fixed=ka+CL+V1+F1+Rate~1, random = pdDiag(ka+CL+V1~1), > # start=c(ka=log(0.05),CL=log(0.5),V1=log(10.0),F1=log(0.8),Rate=log(5)), > # control=list(msVerbose=TRUE,tolerance=1e-3,pnlsTol=1e-1,msTol=1e-3), > # verbose=TRUE) > > #plot(augPred(fit2,level=0:1,length.out=300),aspect=1/1) > > #Estimation of length of infusion > Data$Rate[Data$Rate==-1] <- -2 > > OneCompAbs <- list(DiffEq=list( + dA1dt = ~ -ka*A1, + dA2dt = ~ ka*A1 - CL/V1*A2), + ObsEq=list( + SC= ~0, + C = ~ A2/V1), + States=c("A1","A2"), + Parms=c("ka","CL","V1","F1","Tcrit"), + Init=list(0,0)) > > OneCompAbsModel <- nlmeODE(OneCompAbs,Data) > > #Remove '#' below to run the estimation > > #fit3 <- nlme(Conc ~ OneCompAbsModel(ka,CL,V1,F1,Tcrit,Time,ID), > # data = Data, fixed=ka+CL+V1+F1+Tcrit~1, random = pdDiag(ka+CL+V1~1), > # start=c(ka=log(0.05),CL=log(0.5),V1=log(10.0),F1=log(0.8),Tcrit=log(20)), > # control=list(msVerbose=TRUE,tolerance=1e-3,pnlsTol=1e-1,msTol=1e-3), > # verbose=TRUE) > > ############################################################ > ### Simulation and simultaneous estimation of PK/PD data ### > ############################################################ > > PoolModel <- list( + DiffEq=list( + dy1dt = ~ -ke*y1, + dy2dt = ~ krel * (1-Emax*(y1/Vd)**gamma/(EC50**gamma+(y1/Vd)**gamma)) * y3 - kout * y2, + dy3dt = ~ Kin - krel * (1-Emax*(y1/Vd)**gamma/(EC50**gamma+(y1/Vd)**gamma))*y3), + ObsEq=list( + PK = ~ y1/Vd, + PD = ~ y2, + Pool = ~ 0), + States=c("y1","y2","y3"), + Parms=c("ke","Vd","Kin","kout","krel","Emax","EC50","gamma"), + Init=list(0,"Kin/kout","Kin/krel")) > > ID <- rep(seq(1:12),each=2*12) > Time <- rep(rep(c(0,0.25,0.5,0.75,1,2,4,6,8,10,12,24),each=2),12) > Dose <- rep(c(100,rep(0,23)),12) > Cmt <- rep(rep(c(1,2),12),12) > Type <- rep(rep(c(1,2),12),12) > Conc <- rep(0,2*12*12) > > Data <- as.data.frame(list(ID=ID,Time=Time,Dose=Dose,Cmt=Cmt,Type=Type,Conc=Conc)) > > SimData <- groupedData( Conc ~ Time | ID/Type, + data = Data, + labels = list( x = "Time", y = "Concentration")) > > PKPDpoolModel <- nlmeODE(PoolModel,SimData,JAC=FALSE) > > keSim <- rep(log(rep(0.05,12))+0.1*rnorm(12),each=2*12) > VdSim <- rep(log(rep(10,12))+0.01*rnorm(12),each=2*12) > EC50Sim <- rep(log(rep(5,12))+0.1*rnorm(12),each=2*12) > KinSim <- rep(log(5),2*12*12) > koutSim <- rep(log(0.5),2*12*12) > krelSim <- rep(log(2),2*12*12) > EmaxSim <- rep(log(1),2*12*12) > gammaSim <- rep(log(3),2*12*12) > > SimData$Sim <- PKPDpoolModel(keSim,VdSim,KinSim,koutSim,krelSim,EmaxSim,EC50Sim,gammaSim,SimData$Time,SimData$ID,SimData$Type) > SimData$Conc[SimData$Type==1] <- SimData$Sim[SimData$Type==1] + 0.1*rnorm(length(SimData[SimData$Type==1,1])) > SimData$Conc[SimData$Type==2] <- SimData$Sim[SimData$Type==2] + 0.01*rnorm(length(SimData[SimData$Type==2,1])) > > Data <- groupedData( Conc ~ Time | ID/Type, + data = SimData, + labels = list( x = "Time", y = "Concentration")) > > plot(Data,display=1,aspect=1/1) > > #Fixed parameters > Data$Emax <- rep(log(1),dim(Data)[1]) > > #Estimation of model parameters > PKPDpoolModel <- nlmeODE(PoolModel,Data,JAC=FALSE) > > #Remove '#' below to run the estimation > > #PKPDpool.nlme <- nlme(Conc ~ PKPDpoolModel(ke,Vd,Kin,kout,krel,Emax,EC50,gamma,Time,ID,Type), > # data = Data, fixed=ke+Vd+Kin+kout+krel+EC50+gamma~1, random = pdDiag(ke+Vd+EC50~1), > # groups=~ID, > # weights=varIdent(form=~1|Type), > # start=c(ke=log(0.05),Vd=log(10),Kin=log(5),kout=log(0.5),krel=log(2),EC50=log(5),gamma=log(3)), > # control=list(msVerbose=TRUE,tolerance=1e-1,pnlsTol=1e-1,msTol=1e-1,msMaxIter=20,pnlsMaxIter=20), > # verbose=TRUE) > > > > > ### *