data(PKPD)
#Fixed parameters
PKPD$Emax <- rep(log(1),dim(PKPD)[1])
IRmodel <- list( DiffEq = list(
dy1dt = ~ -ke*y1,
dy2dt = ~ Kin * (1-Emax*(y1/Vd)**gamma/(EC50**gamma+(y1/Vd)**gamma)) - Kin/BL * y2),
ObsEq = list(
PK = ~ y1/Vd,
PD = ~ y2),
States=c("y1","y2"),
Parms=c("ke","Vd","Kin","BL","Emax","EC50","gamma"),
LogParms=TRUE,
Init=list(0,"BL"),
JAC=FALSE,
SEQ=FALSE)
PKPDmodel <- nlmeODE(IRmodel,PKPD)
PKPD.nlme <- nlme(Conc ~ PKPDmodel(ke,Vd,Kin,BL,Emax,EC50,gamma,Time,Subject,Type),
data = PKPD, fixed=ke+Vd+Kin+BL+EC50+gamma~1, random = pdDiag(ke+Vd+EC50~1),
groups=~Subject,
weights=varIdent(form=~1|Type),
start=c(ke=log(0.1),Vd=log(10),Kin=log(5),BL=log(10),EC50=log(5),gamma=log(3)),
control=list(msVerbose=TRUE,tolerance=1e-3,pnlsTol=1e-1,msTol=1e-3,msMaxIter=20,pnlsMaxIter=20))
#Plot results
ni <- 100
TimeSim <- seq(from=0,to=24,length=ni)
TimeSim <- rep(rep(TimeSim,each=2),12)
SubjectSim <- rep(1:12,each=2*ni)
TypeSim <- rep(rep(c(1,2),ni),12)
IndCoef <- coef(PKPD.nlme)
IpredSim <- PKPDmodel( rep(IndCoef[,1],each=2*ni),
rep(IndCoef[,2],each=2*ni),
rep(IndCoef[,3],each=2*ni),
rep(IndCoef[,4],each=2*ni),
rep(rep(log(1),12),each=2*ni),
rep(IndCoef[,5],each=2*ni),
rep(IndCoef[,6],each=2*ni),
TimeSim,SubjectSim,TypeSim)
PopCoef <- fixef(PKPD.nlme)
PredSim <- PKPDmodel( rep(rep(PopCoef[1],12),each=2*ni),
rep(rep(PopCoef[2],12),each=2*ni),
rep(rep(PopCoef[3],12),each=2*ni),
rep(rep(PopCoef[4],12),each=2*ni),
rep(rep(log(1),12),each=2*ni),
rep(rep(PopCoef[5],12),each=2*ni),
rep(rep(PopCoef[6],12),each=2*ni),
TimeSim,SubjectSim,TypeSim)
plotIR <- as.data.frame(rbind(cbind(TimeSim,SubjectSim,PredSim,TypeSim,rep("Pred",2400)),
cbind(TimeSim,SubjectSim,IpredSim,TypeSim,rep("Ipred",2400)),
cbind(PKPD$Time,PKPD$Subject,PKPD$Conc,PKPD$Type,rep("Obs",288))
))
names(plotIR) <- c("Time","Subject","Conc","Type","Flag")
plotIR$Subject <- as.factor(as.numeric(as.character(plotIR$Subject)))
plotIR$Type <- as.factor(plotIR$Type)
plotIR$Flag <- as.factor(plotIR$Flag)
plotIR$Conc <- as.numeric(as.character(plotIR$Conc))
plotIR$Time <- as.numeric(as.character(plotIR$Time))
plotIRPK <- subset(plotIR,Type==1)
plotIRPD <- subset(plotIR,Type==2)
trellis.device(width=7, height=6.5,bg="transparent")
sb <- trellis.par.get("strip.background")
sb$col <- c(16,11,9,13,10,15,14)
trellis.par.set("strip.background",sb)
par(oma=c(0,1,1,0),mgp=c(3,1,.3))
xyplot (Conc~Time | Subject, data=plotIRPK,
layout=c(4,3),
aspect=1/1,
groups=Flag,
grid=TRUE,
xlab="Time since drug administration (hr)",
ylab="PK concentration (ng/mL)",
key=list(x=0.23,y=1.03,corner=c(0,1),transparent=TRUE,
text = list(c("Population", "Individual","Observed")),
lines = list(type=c("l","l","p"), pch=1, col=c(1,1,1), lty=c(1,2,1)),columns=3),
strip=function(...) strip.default(..., strip.names=c(FALSE,TRUE), style=1),
panel = function(x, y, groups,...) {
panel.grid(h=3,v=3,col="lightgray",lwd=0.7,...)
panel.superpose.2(x,y,groups,type=c("l","p","l"), col=c(1,1,1), lty=c(2,1,1),pch=1, lwd=1.4,...)},
par.strip.text=list(cex=1.0))
trellis.device(width=7, height=6.5,bg="transparent")
sb <- trellis.par.get("strip.background")
sb$col <- c(16,11,9,13,10,15,14)
trellis.par.set("strip.background",sb)
par(oma=c(0,1,1,0),mgp=c(3,1,.3))
xyplot (Conc~Time | Subject, data=plotIRPD,
layout=c(4,3),
aspect=1/1,
groups=Flag,
grid=TRUE,
xlab="Time since drug administration (hr)",
ylab="PD concentration (ng/mL)",
key=list(x=0.23,y=1.03,corner=c(0,1),transparent=TRUE,
text = list(c("Population", "Individual","Observed")),
lines = list(type=c("l","l","p"), pch=1, col=c(1,1,1), lty=c(1,2,1)),columns=3),
strip=function(...) strip.default(..., strip.names=c(FALSE,TRUE), style=1),
panel = function(x, y, groups,...) {
panel.grid(h=3,v=3,col="lightgray",lwd=0.7,...)
panel.superpose.2(x,y,groups,type=c("l","p","l"), col=c(1,1,1), lty=c(2,1,1),pch=1, lwd=1.4,...)},
par.strip.text=list(cex=1.0))
Run the code above in your browser using DataLab