inst/doc/RTMB-introduction.R

In RTMB: 'R' Bindings for 'TMB'

## ----setup, include=FALSE, message=FALSE--------------------------------------
##library(knitr)
library(RTMB)
set.seed(1)
formals(MakeADFun)$silent <- TRUE

## -----------------------------------------------------------------------------
library(RTMB)

## -----------------------------------------------------------------------------
data(ChickWeight)

## -----------------------------------------------------------------------------
parameters <- list(
    mua=0,          ## Mean slope
    sda=1,          ## Std of slopes
    mub=0,          ## Mean intercept
    sdb=1,          ## Std of intercepts
    sdeps=1,        ## Residual Std
    a=rep(0, 50),   ## Random slope by chick
    b=rep(0, 50)    ## Random intercept by chick
)

## -----------------------------------------------------------------------------
f <- function(parms) {
    getAll(ChickWeight, parms, warn=FALSE)
    ## Optional (enables extra RTMB features)
    weight <- OBS(weight)
    ## Initialize joint negative log likelihood
    nll <- 0
    ## Random slopes
    nll <- nll - sum(dnorm(a, mean=mua, sd=sda, log=TRUE))
    ## Random intercepts
    nll <- nll - sum(dnorm(b, mean=mub, sd=sdb, log=TRUE))
    ## Data
    predWeight <- a[Chick] * Time + b[Chick]
    nll <- nll - sum(dnorm(weight, predWeight, sd=sdeps, log=TRUE))
    ## Get predicted weight uncertainties
    ADREPORT(predWeight)
    ## Return
    nll
}

## -----------------------------------------------------------------------------
obj <- MakeADFun(f, parameters, random=c("a", "b"))

## -----------------------------------------------------------------------------
opt <- nlminb(obj$par, obj$fn, obj$gr)

## -----------------------------------------------------------------------------
sdr <- sdreport(obj)
sdr

## -----------------------------------------------------------------------------
set.seed(1)
chk <- checkConsistency(obj)
chk

## -----------------------------------------------------------------------------
osa <- oneStepPredict(obj, method="fullGaussian", discrete=FALSE)
qqnorm(osa$res); abline(0,1)

## -----------------------------------------------------------------------------
f2 <- function(parms) {
    getAll(ChickWeight, parms, warn=FALSE)
    ## Optional (enables extra RTMB features)
    weight <- OBS(weight)
    ## Random slopes
    a %~% dnorm(mean=mua, sd=sda)
    ## Random intercepts
    b %~% dnorm(mean=mub, sd=sdb)
    ## Data
    predWeight <- a[Chick] * Time + b[Chick]
    weight %~% dnorm(predWeight, sd=sdeps)
    ## Get predicted weight uncertainties
    ADREPORT(predWeight)
}

## -----------------------------------------------------------------------------
obj <- MakeADFun(f2, parameters, random=c("a", "b"))

## -----------------------------------------------------------------------------
f3 <- function(parms, data) {
    getAll(data, parms, warn=FALSE)
    ## Optional (enables extra RTMB features)
    weight <- OBS(weight)
    ## Random slopes
    a %~% dnorm(mean=mua, sd=sda)
    ## Random intercepts
    b %~% dnorm(mean=mub, sd=sdb)
    ## Data
    predWeight <- a[Chick] * Time + b[Chick]
    weight %~% dnorm(predWeight, sd=sdeps)
    ## Get predicted weight uncertainties
    ADREPORT(predWeight)
}

## -----------------------------------------------------------------------------
cmb <- function(f, d) function(p) f(p, d)

## -----------------------------------------------------------------------------
## Using the original ChickWeight
obj <- MakeADFun(cmb(f3, ChickWeight), parameters, random=c("a", "b"))
## Using a new dataset
ChickWeightNew <- transform(ChickWeight, weight=log(weight))
obj <- MakeADFun(cmb(f3, ChickWeightNew), parameters, random=c("a", "b"))

## -----------------------------------------------------------------------------
f_pred <- function(parms, data) {
    getAll(data, parms, warn=FALSE)
    ## Optional (enables extra RTMB features)
    weight <- OBS(weight)
    ## Random slopes
    a %~% dnorm(mean=mua, sd=sda)
    ## Random intercepts
    b %~% dnorm(mean=mub, sd=sdb)
    ## Data
    predWeight <- a[Chick] * Time + b[Chick]
    mis <- is.na(weight)                                 ## <-- Get missing responses
    weight[!mis] %~% dnorm(predWeight[!mis], sd=sdeps)   ## <-- Likelihood of non-missing
    ## Get predicted weight uncertainties
    ADREPORT(predWeight[mis])                            ## <-- Only report missing
}

## -----------------------------------------------------------------------------
newdata <- data.frame(Time=20:24, Diet="3", Chick="10", weight=NA)

## -----------------------------------------------------------------------------
augdata <- rbind(ChickWeight, newdata)
obj <- MakeADFun(cmb(f_pred, augdata), parameters, random=c("a", "b"))
opt <- nlminb(obj$par, obj$fn, obj$gr)
sdr <- sdreport(obj)
cbind(newdata, summary(sdr, select="report"))

## -----------------------------------------------------------------------------
newdata <- transform(newdata, Chick="51")
augdata <- rbind(ChickWeight, newdata)
parameters$a <- rep(0, nlevels(augdata$Chick))
parameters$b <- rep(0, nlevels(augdata$Chick))

## -----------------------------------------------------------------------------
par.old <- opt$par
hessian.old <- optimHess(par.old, obj$fn, obj$gr)

## -----------------------------------------------------------------------------
obj <- MakeADFun(cmb(f_pred, augdata), parameters, random=c("a", "b"))
sdr <- sdreport(obj, par.old, hessian.old)
cbind(newdata, summary(sdr, select="report"))