# test suite for Chapter 10 of Statistical Models in S [CXF April 1992]

# Puromycin tests

{Treated <- Puromycin[Puromycin$state == "treated",]; T}
{Untreated <- Puromycin[Puromycin$state == "untreated",]; T}

{Purfit1 <- nls( vel ~ Vm*conc/(K+conc), Treated, list( Vm=200, K=0.1)); T}
{Purfit2 <- nls( vel ~ Vm*conc/(K+conc), Untreated, list( Vm=200, K=0.1)); T}
{Purboth <- nls( vel ~ (Vm + delV*(state=="treated"))*conc/(K+conc), 
 Puromycin, list( Vm=160, delV=40, K=0.05)); T}

{combinedSS <- sum(Purfit1$res^2) + sum(Purfit2$res^2); T}
{Fval <- (sum(Purboth$res^2) - combinedSS) / (combinedSS/19); T}
{abs((1-pf(Fval, 1, 19)) - 0.205552) < .000001}
{rm( combinedSS, Fval, Untreated); T}

{Purfit1.summary <- summary(Purfit1); T}
{Purfit2.summary <- summary(Purfit2); T}
{Purboth.summary <- summary(Purboth); T}
{rm( Purfit1, Purfit2, Purfit1.summary, Purfit2.summary, Purboth.summary); T}

{Pur.pro <- Puromycin; T}
{parameters(Pur.pro) <- list(delV = 40); T}
{Pur.40 <- nls(vel ~ (Vm+delV*(state=="treated"))*conc/(K+conc), Pur.pro, 
 list( Vm=160, K=0.05)); T}
{abs(sum(Pur.40$res^2) - 2252.551) < .001}
{rm( Pur.40, Pur.pro); T}

{Pur.pro <- profile(Purboth); T}
{all(sort(names(Pur.pro)) == sort(c("Vm", "delV", "K")))}
{rm(Pur.pro); T}

{Purfit.pl <- nls(vel~conc/(K+conc), Treated, list(K=0.1), alg="plinear"); T}
{abs(sum(Purfit.pl$res^2) - 1195.449) < .001}
{rm(Purfit.pl); T}

{weighted.MM <- function ( resp, conc, Vm, K)
                        { pred <- (Vm*conc)/(K+conc)
                         (resp - pred)/sqrt(pred) }
 T
}

{Pur.wt <- nls( ~ weighted.MM( vel, conc, Vm, K), Treated, 
                  list( Vm = 200, K = 0.05)); T}
{abs(sum(Pur.wt$res^2) - 14.6) < .01}
{rm(Pur.wt, weighted.MM); T}

{Pur.TBS <- nls( ~ TBS( vel, Vm*conc/(K+conc), 0), Treated, 
                   list( Vm = 200, K = 0.05)); T}
{abs(sum(Pur.TBS$res^2) - .1677) < .00001}
{rm(Pur.TBS, Treated); T}

# pingpong tests

{param(pingpong, "p") <- 0; T}
{attach( pingpong, 1); T}
{D <- winner - loser; T}
{p <- sum(winner>loser)/length(winner); T}
{abs(p - .82234) < .000001}
{alpha <- log(p/(1-p))/mean(D); T}
{abs(alpha - .007661) < .000001}
{detach(1, save = "pingpong"); T}
{lprob <- function(lp)log(1+exp(lp))-lp; T}
{fit.alpha <- ms( ~ lprob(D*alpha), pingpong); T}
{abs(fit.alpha$value - 1127.63) < .01}
{rm( pingpong, fit.alpha, lprob); T}

# Lubricant tests

{Lub.mod <- function(temp, press, t2, t8, t9)
        {efac <- press * exp(( - temp)/(t8 + t9 * press^2))
         c(1/(t2 + temp), press, press^2, press^3, efac, efac * press^2)}
 T
}

{parameters(Lubricant) <- list( t2 = 100, t8 = 100, t9 = 0); T}
#{Lubfit <- nls( viscos ~ Lub.mod( tempC, pressure, t2, t8, t9), Lubricant, 
#  alg = "plinear"); T}
{Lubfit <- nls( viscos ~ Lub.mod( tempC, pressure, t2, t8, t9), Lubricant, 
  list(t2 = 100), alg = "plinear"); T}
{abs(Lubfit$parameters[["t2"]] - 219) < .1}
{Lubfit <- nls( viscos ~ Lub.mod( tempC, pressure, t2, t8, t9), Lubricant, 
 list(t2 = 219, t8 = 100), alg = "plinear"); T}
{(abs(Lubfit$parameters["t2"] - 209.5) < .1)}
{(abs(Lubfit$parameters["t8"] - 47.7) < .1)}
{Lubfit <- nls( viscos ~ Lub.mod( tempC, pressure, t2, t8, t9), Lubricant, 
 list(t2 = 210, t8 = 48, t9 = 0), alg = "plinear"); T}
{rm(Lubfit, Lub.mod, Lubricant); T}

# D() and deriv() test (just tests for bombing)

{temp <- D( substitute( log(1+exp(D*alpha)) - D*alpha), "alpha"); T}
{temp <- deriv( ~ log(1+exp(D*alpha)) - D*alpha, "alpha"); T}
{temp <- deriv( ~ log(1+exp(D*alpha)) - D*alpha, "alpha", c("D","alpha")); T}
{temp <- deriv( ~ log(1+exp(D*alpha)) - D*alpha, "alpha", 
 function(D=1,alpha=0)NULL); T}
{temp <- deriv( vel ~ Vm*(conc/(K+conc)), c("Vm","K")); T}
{rm(temp); T}