## nonlinear mixed models --- 3-part formulas ---
## 1. basic nonlinear fit. Use stats::SSlogis for its
## implementation of the 3-parameter logistic curve.
## "SS" stands for "self-starting logistic", but the
## "self-starting" part is not currently used by nlmer ... 'start' is
## necessary
startvec <- c(Asym = 200, xmid = 725, scal = 350)
(nm1 <- nlmer(circumference ~ SSlogis(age, Asym, xmid, scal) ~ Asym|Tree,
Orange, start = startvec))
## 2. re-run with "quick and dirty" PIRLS step
(nm1a <- update(nm1, nAGQ = 0L))
## 3. Fit the same model with a user-built function:
## a. Define formula
nform <- ~Asym/(1+exp((xmid-input)/scal))
## b. Use deriv() to construct function:
nfun <- deriv(nform,namevec=c("Asym","xmid","scal"),
function.arg=c("input","Asym","xmid","scal"))
nm1b <- update(nm1,circumference ~ nfun(age, Asym, xmid, scal) ~ Asym | Tree)
## 4. User-built function without using deriv():
## derivatives could be computed more efficiently
## by pre-computing components, but these are essentially
## the gradients as one would derive them by hand
nfun2 <- function(input, Asym, xmid, scal) {
value <- Asym/(1+exp((xmid-input)/scal))
grad <- cbind(Asym=1/(1+exp((xmid-input)/scal)),
xmid=-Asym/(1+exp((xmid-input)/scal))^2*1/scal*
exp((xmid-input)/scal),
scal=-Asym/(1+exp((xmid-input)/scal))^2*
-(xmid-input)/scal^2*exp((xmid-input)/scal))
attr(value,"gradient") <- grad
value
}
stopifnot(all.equal(attr(nfun(2,1,3,4),"gradient"),
attr(nfun(2,1,3,4),"gradient")))
nm1c <- update(nm1,circumference ~ nfun2(age, Asym, xmid, scal) ~ Asym | Tree)
Run the code above in your browser using DataLab