SC_HW4_01

(a)

data_hw4 <- read.table("C:/RClass/serum_conc.txt", header=TRUE)
attach(data_hw4)

model_hw <- 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)
  attr(f,'gradient') <- grad
  f
}

### Initial Value
l_y = log(Concentration+1)
lm(l_y ~ Time)

## 
## Call:
## lm(formula = l_y ~ Time)
## 
## Coefficients:
## (Intercept)         Time  
##      5.0446      -0.1147

### alpha1(ini)=exp(5.0446)=155.1822 , lam1(ini)=-(-0.1147)
l_y1 = log(Concentration-(155.1822)*exp((-0.1147)*Time))

## Warning in log(Concentration - (155.1822) * exp((-0.1147) * Time)): NaNs
## produced

lm(l_y1 ~ Time) 

## 
## Call:
## lm(formula = l_y1 ~ Time)
## 
## Coefficients:
## (Intercept)         Time  
##       4.644       -1.628

### alpha2(ini)=exp(4.644)=103.9594 , lam2(ini)=-(-1.628)

### nls() Funtion
### model_hw(alpha1=155.1822, alpha2=103.9594, lam1=0.1147, lam2=1.628, t=seq(0,50,.1))

Result_hw <- nls(Concentration ~ model_hw(alpha1, alpha2, lam1, lam2, Time),
            start=list(alpha1=155.1822, alpha2=103.9594, lam1=0.1147, 
            lam2=1.628), trace = TRUE, nls.control(maxiter = 50, tol = 1e-5, 
            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

## 888.4771    (4.71e+00): par = (155.1822 103.9594 0.1147 1.628)
## 56.46580    (8.09e-01): par = (156.7234 83.59623 0.1491313 1.18452)
## 34.54546    (6.93e-02): par = (162.7865 81.07917 0.16178 1.303858)
## 34.38027    (4.97e-04): par = (162.6057 81.23797 0.1617988 1.306314)
## 34.38026    (5.76e-05): par = (162.5971 81.2417 0.1617897 1.306037)
## 34.38026    (6.70e-06): par = (162.5981 81.24125 0.1617908 1.306069)

summary(Result_hw)

## 
## Formula: Concentration ~ model_hw(alpha1, alpha2, lam1, lam2, Time)
## 
## Parameters:
##         Estimate Std. Error t value Pr(>|t|)    
## alpha1 162.59810    6.94316  23.418 1.18e-08 ***
## alpha2  81.24125    6.10751  13.302 9.74e-07 ***
## lam1     0.16179    0.00878  18.427 7.75e-08 ***
## lam2     1.30607    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: 5 
## Achieved convergence tolerance: 6.695e-06

\[y=(162.59810)e^{-0.16179t}+(81.24125)e^{-1.30607t}\]

(b)

plot(Time, Concentration)
### One component model(in class/red)
t = seq(0,50,.1)
lines(t,211.9203*exp(-0.2357*t),col='red')
### Two component model(blue)
alpha1 = 162.59810 
alpha2 = 81.24125 
lam1 = 0.16179 
lam2 = 1.30607
E1 = exp(-lam1*t) 
E2 = exp(-lam2*t)
t_hw = seq(0,50,.1)
lines(t_hw,(alpha1*E1)+(alpha2*E2),col='blue')

=> The two component model(blue) fits better than the one component model(red).