# Time-course of drug amounts in the different compartments of a pharmacokinetic system

I am looking for a way to speed up applying my function to every individual in the data frame. Currently I am using the ddply but it is very slow for big data sets. To give a simple example, I am providing the code below:

The first part of the code generates an example data frame for only 2 subjects. The data frame contains dosing amounts AMT and parameters CL, V2,Q,V3,KA,F1 for each subject that are used in the function calculation.

#---------------------------------
#Run this part to generate the data frame to be processed by the function
dosetimes <- c(0,12)
sampletimes <- sort(unique(c(seq(from=0,to=100,by=1))))
#number of subjects
ID <- 1:10
#Make dataframe
df <- expand.grid("ID"=ID,"TIME"=sampletimes, "AMT"=0)
doserows <- subset(df, TIME%in%dosetimes)
doserows$AMT <- 100 #Add back dose information df <- rbind(df,doserows) df <- df[order(df$ID,df$TIME),] df <- subset(df, (TIME==0 & AMT==0)==F) df$CL[df$ID<=1] <- 2 df$CL[df$ID>=2] <- 2.5 df$V2[df$ID<=1] <- 10 df$V2[df$ID>=2] <- 15 df$Q [df$ID<=1] <- 20 df$Q [df$ID>=2] <- 20 df$V3[df$ID<=1] <- 30 df$V3[df$ID>=2] <- 30 df$KA[df$ID<=1] <- 0.2 df$KA[df$ID>=2] <- 0.25 df$F1[df$ID<=1] <- 1 df$F1[df$ID>=2] <- 0.5 #calculate micro-parameters df$k20 <- df$CL/df$V2
df$k23 <- df$Q/df$V2 df$k32 <- df$Q/df$V3


The following is the function that calculates Amounts (A1,A2,A3) by advancing the solution from one TIME point in the df to the next. This function is called by the ddply later. Please don't panic; these are just equations for A1,A2 and A3.

   #-----------------------------------
# This is the function to provess the data frame
TwoCompOral <- function(d){
#Accepts a NONMEM style data frame for 1 subject with columns for TIME, AMT,MDV,DV, CL, V2, Q, V3, KA & F1
#Returns a dataframe with populated columns for A1, A2, A3 and DV

#set initial values in the compartments
d$A1[d$TIME==0] <- d$AMT[d$TIME==0]*d$F1[1] # Amount in the absorption compartment at time zero. d$A2[d$TIME==0] <- 0 # Amount in the central compartment at time zero. d$A3[d$TIME==0] <- 0 # Amount in the peripheral compartment at time zero. #This loop calculates micro-rate constants based on individual's PK parameter values. #It uses these values to calculate macro-rate constants (Lambda1/lambda2). #Rate constants(micro- and macro), along other parameters, are used in the equations to calculate drug amounts in each compartment. #The loop advances the solution from one time interval to the next. k20 <- d$k20[1]
k23 <- d$k23[1] k32 <- d$k32[1]
KA <- d$KA[1] k30 <- 0 E2 <- k20+k23 E3 <- k32+k30 #calculate hybrid rate constants lambda1 = 0.5*((E2+E3)+sqrt((E2+E3)^2-4*(E2*E3-k23*k32))) lambda2 = 0.5*((E2+E3)-sqrt((E2+E3)^2-4*(E2*E3-k23*k32))) for(i in 2:nrow(d)) { t <- d$TIME[i]-d$TIME[i-1] A2last <- d$A2[i-1]
A3last <- d$A3[i-1] A1last <- d$A1[i-1]

A2term1 = (((A2last*E3+A3last*k32)-A2last*lambda1)*exp(-t*lambda1)-((A2last*E3+A3last*k32)-A2last*lambda2)*exp(-t*lambda2))/(lambda2-lambda1)
A2term2 = A1last*KA*(exp(-t*KA)*(E3-KA)/((lambda1-KA)*(lambda2-KA))+exp(-t*lambda1)*(E3-lambda1)/((lambda2-lambda1)*(KA-lambda1))+exp(-t*lambda2)*(E3-lambda2)/((lambda1-lambda2)*(KA-lambda2)))
d$A2[i] = A2term1+A2term2 #Amount in the central compartment A3term1 = (((A3last*E2+A2last*k23)-A3last*lambda1)*exp(-t*lambda1)-((A3last*E2+A2last*k23)-A3last*lambda2)*exp(-t*lambda2))/(lambda2-lambda1) A3term2 = A1last*KA*k23*(exp(-t*KA)/((lambda1-KA)*(lambda2-KA))+exp(-t*lambda1)/((lambda2-lambda1)*(KA-lambda1))+exp(-t*lambda2)/((lambda1-lambda2)*(KA-lambda2))) d$A3[i] = A3term1+A3term2  #Amount in the peripheral compartment

A1last = A1last*exp(-t*KA)
d$A1[i] = A1last + d$AMT[i]*d$F1[1] #Amount in the absorption compartment d$DV[i] <- d$A2[i]/d$V2[i]        #Concentration in the central compartment

}
d
}


The following is what I am using to process the data frame df I generated for the subjects:

#Apply the function to each ID
system.time(simdf <- ddply(df, .(ID), TwoCompOral))


All warning that I had before are fixed.

The problem is that the function TwoCompOral becomes very slow when applied using ddply to, for example, 1000 subjects and when the number of TIME points increases per subject (because of the for-loop). I want a way(s) to speed it up.

I want to know the following:

1. How can I apply the function TwoCompOral on the df using the dplyr package: probably it is faster than the plyr?
2. Would you suggest other ways?
• Welcome to Code Review! I've made a guess at a suitable title for this question. If I've guessed wrong, please edit the question to clarify the intentions of your code, so that we don't have to reverse-engineer it. MathJax is available, if you need to write mathematical expressions. – 200_success Mar 25 '15 at 2:15
• I have reopened this question, but have two concerns about it still. The request to use dplyr sounds a lot like "gimme the code", and if your intention is to use dplyr, then your code is currently not in a working state because it does not work as intended. The second concern I have is that while you have a lot of text describing low-level technical detail, you have nothing describing anything about the actual context of this code - what problem does it solve? I believe the question is now minimally on-topic, but it makes it unnecessarily hard to review. – rolfl Mar 27 '15 at 10:31
• Your overall code is slow for two reasons: (1) it needs to loop on all the IDs and (2) for each ID, it needs to run TwoCompOral which is slow because of its for loop. plyr and dplyr help with (1) but I am worried the real problem is with (2). You can test if (1) or (2) is the real problem by running TwoCompOral on an average individual, then multiply the runtime by the number of individuals. If you get a number close to the runtime you get with ddply for all individuals, then you know the problem is not with ddply but with (2). And dplyr won't help. – flodel Mar 27 '15 at 10:56
• So if TwoCompOral and its for loop are the real problem, how can we improve it? At a glance, it does not look easy as your system appears to be sequential (computation for row i requires the results of the same computations for row i-1). In R, there are only a few fast functions to help you with sequential processes: cumsum, cummax, filter, diff, head/tail, ... It's difficult but see if you can rewrite your code in a vectorized way using these functions. – flodel Mar 27 '15 at 11:01
• ...If not, then R might not be the best language for this task. There is always the Rcpp package for dealing with this within R but the learning curve is steep. – flodel Mar 27 '15 at 11:01