modified code [🇷 for BE/BA]

posted by martin  – Austria, 2008-10-16 14:29 (6445 d 14:13 ago) – Posting: # 2543
Views: 26,151

dear ace and yung-jin lee !

I suggest not to use the absolute value. please find below a new modification of the ace's code including some lines to prevent errors which I encountered in a simulation study (for strange PK profiles)

best regards

# function for selection of time points
tp <- function(dat){

   dat <- dat[order(dat$time),] # modification

   m <- which.max(dat$conc)
   f2 <-  function(m) return(cbind((nrow(dat)-m+1),abs(extractAIC(lm(log(conc)~time,dat[m:nrow(dat),])))[2],summary(lm(log(conc)~time,dat[m:nrow(dat),]))$adj.r.squared))
   overview <- as.data.frame(do.call(rbind,lapply((m+1):(nrow(dat)-2),f2)))
   names(overview) <- c("n","AIC","adjR2")
   n_ARS=0
   r.adj=0
   for (i in (nrow(dat)-2):(which.max(dat$conc)+1)) {
      check <- summary(lm(log(conc)~time,dat[i:nrow(dat),]))$adj.r.squared # modification
      if (!is.na(check) && (r.adj - check <(0.0001))) {
         n_ARS = nrow(dat)-i+1
         r.adj = summary(lm(log(conc)~time,dat[i:nrow(dat),]))$adj.r.squared
      }
   }   

   n_TTT_ARS=0
   r.adj2=0
   for (i in (nrow(dat)-2):(min(seq_along(dat$time)[dat$time>=dat$time[which.max(dat$conc)]*2]))) {
      check <- summary(lm(log(conc)~time,dat[i:nrow(dat),]))$adj.r.squared # modification
      if (!is.na(check) && (r.adj2 - check < (0.0001))) {
         n_TTT_ARS = nrow(dat)-i+1
         r.adj2 = summary(lm(log(conc)~time,dat[i:nrow(dat),]))$adj.r.squared
      }
   }   
   
   # start modification
   require(PK)
   leedat <- dat[c((m+1):nrow(dat)),]
   n_lee <- nrow(leedat)
   if(nrow(leedat) >= 5){
       l <- lee(conc=leedat$conc, time=leedat$time, method='ols', lt=FALSE, points=3)
       n_lee <- sum(dat$time> l$chgpt)
       if(is.na(n_lee)){n_lee <- nrow(leedat)}
   }
   # end modification

   n_TTT <- sum(dat$time> (dat$time[which.max(dat$conc)]*2))
   n_AIC <- overview$n[which.min(overview$AIC)] # modification

   #plot(l,log="y")
   print(overview)
   cat("n")
   return(c(TTT=n_TTT, last3=3, AIC=n_AIC, ARS=n_ARS,TTT_ARS=n_TTT_ARS,lee=n_lee))
}

Complete thread:

UA Flag
Activity
 Admin contact
23,653 posts in 4,991 threads, 1,570 registered users;
119 visitors (0 registered, 119 guests [including 25 identified bots]).
Forum time: 04:43 CEST (Europe/Vienna)

In theory, there is no difference between theory and practice.
But, in practice, there is.    Jan L.A. van de Snepscheut

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5