Exercise J-2.3 Part A(i.): Our model and Initial Values

We use nonlinear regression to fit two-term clearance model: \[ y = \alpha_1e^{-\lambda_1t} + \alpha_2e^{-\lambda_2t} + \varepsilon \]

Before running our algorithm we must first come up with initial values. To do so we first suppose that our model is only a one-term clearance model. In other words, assume \[ y = \alpha_1e^{\lambda_1t} + \varepsilon \]

We take log of each side and arrive at \[ \log{y} = \log{\alpha_1} - \lambda_1t \]

We fit a linear model using the R command lm(log y ~ time) and come up with the following: \[ \lambda_1 \approx 0.1147\\ \log{\alpha_1} \approx 5.04 \implies \alpha_1 \approx 154 \]

Now with these estimates we go back to our original two-term clearance model: \[ y = 154e^{-0.1147t} + \alpha_2e^{-\lambda_2t} \]

We replace \(t\) and \(y\) with the data we have. The expression \(154e^{-0.1147t}\) becomes some constant \(\hat{y}\), and substracting it to the left side of the equation results in some constant \(c = y - \hat{y}\). Then we have the following:

\[ c = \alpha_2e^{-\lambda_2t} \]

We again take log of each side and run a linear model lm(log(c) ~ Time) in R. Again we come up with the following estimates: \[ \lambda_2 \approx 1.627\\ \log{\alpha_2} \approx 4.623 \implies \alpha_2 \approx 103.94 \] Therefore, we have the following initial value for our algorithm: \[ (\alpha_1, \alpha_2, \lambda_1, \lambda_2)_{initial} = (154, 103, 0.11, 1.627)^{T} \]

Exercise J-2.3 Part A(ii.): Run nls() and results

data = read.table("serum_conc.txt", header = TRUE)
attach(data)

model2 <- function(alpha1, alpha2, lam1, lam2, t) {
    E1 = exp(-lam1 * t)
    E2 = exp(-lam2 * t)
    f = alpha1 * E1 + alpha2 * E2
    grad = cbind(E1, E2, -alpha1 * t * E1, -alpha2 * t * E2)
    # grad = c(E1, E2, -alpha1*t*E1, -alpha2*t*E2)
    attr(f, "gradient") <- grad
    f
}

##### Estimate initial values of alpha1 and lambda1
l.y = log(Concentration + 1)
model.linear = lm(l.y ~ Time)
model.linear
## 
## Call:
## lm(formula = l.y ~ Time)
## 
## Coefficients:
## (Intercept)         Time  
##      5.0446      -0.1147
alpha1.init = exp(model.linear$coefficients[1])  #5.044588,  #155.1803
lambda1.init = -1 * model.linear$coefficients[2]  #0.1147291


y.hat = alpha1.init * exp(-1 * lambda1.init * Time)
res = Concentration - y.hat
l.r = log(res)
## Warning in log(res): NaNs produced
model.linear.2 = lm(l.r ~ Time)
alpha2.init = exp(model.linear.2$coefficients[1])  #103.9414
lambda2.init = -1 * model.linear.2$coefficients[2]  # 1.627351
##### 

Result2 = nls(Concentration ~ model2(alpha1, alpha2, lam1, lam2, Time), start = list(alpha1 = 155, 
    alpha2 = 103, lam1 = 0.11, lam2 = 1.627), trace = TRUE, nls.control(maxiter = 50, 
    tol = 1e-05, minFactor = 1/1024, printEval = TRUE))
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## 1133.086 :  155.000 103.000   0.110   1.627
## 68.25636 :  155.306180  84.501484   0.146178   1.145490
## 34.94291 :  162.8403329  80.9460229   0.1616855   1.2976770
## 34.38037 :  162.6290416  81.2265586   0.1618231   1.3070259
## 34.38026 :  162.5944895  81.2428148   0.1617869   1.3059539
## 34.38026 :  162.5984066  81.2411192   0.1617911   1.3060788
## 34.38026 :  162.5979478  81.2413203   0.1617906   1.3060643
summary(Result2)
## 
## Formula: Concentration ~ model2(alpha1, alpha2, lam1, lam2, Time)
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## alpha1 162.59795    6.94319  23.418 1.18e-08 ***
## alpha2  81.24132    6.10752  13.302 9.74e-07 ***
## lam1     0.16179    0.00878  18.427 7.75e-08 ***
## lam2     1.30606    0.19757   6.611 0.000168 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.073 on 8 degrees of freedom
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 3.017e-06

Our R console printout states our estimates for our parameters. In particular we have the following: \[ \alpha_1 \approx 162.59795\\ \lambda_1 \approx 0.16179\\ \alpha_2 \approx 81.24132\\ \lambda_2 \approx 1.30606\\ \implies y = 162.59795e^{-0.16179t} + 81.24132e^{-1.30606t} \]

Exercise J-2.3 Part B

Below we plot the data, the one-term clearance model fit from Lecture (in red) and our two-term clearance model with the above estimates (in blue)

plot(Time, Concentration)
t = seq(0, 50, 0.1)
lines(t, 162.59795 * exp(-0.16179 * t) + 81.24132 * exp(-1.30606 * t), col = "blue")
lines(t, 211.9203 * exp(-0.2357 * t), col = "red")  #one-term clearance model fit in class

It is easy to see that our blue model does much better.