library("datasets")
attach(Indometh)

# --------------------------------------------------------------------
# Calculate the rate of decay (k) for each person and the global
# average.
# --------------------------------------------------------------------

# this is our function to calculate k
# we keep it separate so that the logic is clear
k.calc <- function(times, concs) {
	rv <- lm(log(concs) ~ times);
	-1 * as.numeric(rv$coefficients["times"])
}

# here's where our results will go
ks <- list()

# calculate k for each subject no.
all.times <- tapply(time, Subject, c)
all.concs <- tapply(conc, Subject, c)
for (i in sort(names(all.times)))
	ks[[i]] <- k.calc(all.times[[i]], all.concs[[i]])

# calculate the global k value
ks[["global"]] <- k.calc(time, conc)

# calculate the half-life
ks.hl <- log(2) / as.numeric(ks)


# --------------------------------------------------------------------
# this will plot the Indometh data along with a nice 
# regression curve.
# --------------------------------------------------------------------

# plot the original curves
plot(time[which(Subject==1)], conc[which(Subject==1)], "l", col="red", main="Pharmacokinetics of Indomethicin", xlab="time (hr)", ylab="concentration (mcg/ml)")
points(time[which(Subject==2)], conc[which(Subject==2)], "b", col="red")
points(time[which(Subject==3)], conc[which(Subject==3)], "b", col="red")
points(time[which(Subject==4)], conc[which(Subject==4)], "b", col="red")
points(time[which(Subject==5)], conc[which(Subject==5)], "b", col="red")
points(time[which(Subject==6)], conc[which(Subject==6)], "b", col="red")

# do regression
fit <- lm(log(conc) ~ time)

# get snappy new data points
time.pts <- time[which(Subject==1)]
time.coeff <- as.numeric(fit$coefficients["time"])
fit.pts <- exp(time.coeff * time.pts)

# plot the new curve
lines(time.pts, fit.pts)

# release the data
detach(Indometh)