

# Weighted Michaelis-Menten model
# with data = a list that can not be coerced to a data.frame 
Treated <- Puromycin[Puromycin$state == "treated", ]
TreatIrreg <- with(Treated,
        list(conc1=conc[1], conc.1=conc[-1], rate=rate))
sapply(TreatIrreg, length)

# Passing arguments via 'get'
weighted.MM0 <- function(Vm, K){  
  TI <- get("TreatIrreg")
  conc <- with(TI, c(conc1, conc.1))
  resp <- TI$rate
#
  pred <- (Vm * conc)/(K + conc)
  (resp - pred) / sqrt(pred)
}
gc()
start.time <- proc.time()
Pur.wt0 <- nls( ~ weighted.MM0(Vm, K), start = list(Vm = 200, K = 0.1),
                   trace = TRUE)
(et0 <- proc.time()-start.time)
# NOTE:  'get' works but is not recommended,
# because it's too subtle and inflexible.  

# Passing arguments using a list that can not be coerced to a data.frame
weighted.MM1 <- function(resp, conc1, conc.1, Vm, K){
  conc <- c(conc1, conc.1)
#
  pred <- (Vm * conc)/(K + conc)
  (resp - pred) / sqrt(pred)
}
gc()
start.time <- proc.time()
Pur.wt1 <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K),
     data=TreatIrreg, start = list(Vm = 200, K = 0.1), trace = TRUE)
(et1 <- proc.time()-start.time)

gc()
start.time <- proc.time()
Pur.wt.1port <- nls( ~ weighted.MM1(rate, conc1, conc.1, Vm, K),
                    data=TreatIrreg, start = list(Vm = 200, K = 0.1),
                    algorithm="port", trace = TRUE)
(et.1port <- proc.time()-start.time)

# Chambers and Hastie (1992) Statistical Models in S (Wadsworth, p. 537)
# say, "If the value of the right side [of formula] has an attribute
# called 'gradient' this should be a matrix with the number of rows
# equal to the length of the response and one column for each
# parameter.  
weighted.MM.gradient <- function(resp, conc1, conc.1, Vm, K){
  conc <- c(conc1, conc.1)
#
  K.conc <- (K+conc)
  dy.dV <- conc/K.conc
  dy.dK <- (-Vm*dy.dV/K.conc)
  pred <- Vm*dy.dV
# Avoid NAs   
  pred.5 <- sqrt(pmax(0, pred))
  dev <- (resp - pred) / pred.5
  Ddev <- (-0.5*(resp+pred)/(pred.5*pred))
  attr(dev, "gradient") <- (Ddev * cbind(Vm = dy.dV, K = dy.dK))
  dev 
}
gc()
start.time <- proc.time()
Pur.wt.grad <- nls(
   ~ weighted.MM.gradient(rate, conc1, conc.1, Vm, K), data=TreatIrreg, 
   start = list(Vm = 200, K = 0.1), trace = TRUE)
(et.grad <- proc.time()-start.time)

