In bbolker/fitsir: SIR fitting tools

SIR

library(numDeriv)
library(bbmle)
library(fitsir)
## these aren't exported, probably should be
findSSQ <- fitsir:::findSSQ
findSens <- fitsir:::findSens
fitsir.optim <- fitsir:::fitsir.optim
library(deSolve)
library(emdbook)
## source("../R/fitSIR_funs.R")
library("devtools")
load_all("..")
bombay2 <- setNames(bombay, c("tvec", "count"))
fpars <- coef(f1 <- fitsir(bombay2, start = startfun()))

jacobian(fitsir:::findSSQ, fpars, data = bombay2, SSQonly = TRUE)
findSens(bombay2, fpars, sensOnly = TRUE)

Jacobian seems to work fairly well. Can we find the Hessian?

## fully finite-difference (Richardson)
(hess.m <- hessian(fitsir:::findSSQ, fpars, data = bombay2, SSQonly = TRUE))
## finite difference of symbolic gradient
(hess.grad <- jacobian(fitsir:::findSens, fpars, data = bombay2, sensOnly = TRUE))

Here's the model:

$$ \begin{aligned} &\frac{dSSQ}{d\log \theta}\ &= \frac{dSSQ}{d\theta} \frac{d\theta}{d\log \theta}\ &= \frac{dSSQ}{dI} \frac{dI}{d\theta} \frac{d\theta}{d\log\theta}\ &= 2(I-\hat{I}) \frac{dI}{d\theta} \theta \end{aligned} $$

Second derivative:

$$ \begin{aligned} &\frac{d(dSSQ/d\log\theta_1)}{d\log\theta_2}\ &= \frac{d(dSSQ/d\log\theta_1)}{d\theta_2} \theta_2\ &= (2\frac{dI}{d\theta_2} \frac{dI}{d\theta_1} \theta_1 + 2(I-\hat{I}) \frac{d^2I}{d\theta_1\theta_2}\theta_1 + 2(I-\hat{I})\frac{dI}{d\theta_1} \frac{d\theta_1}{d\theta_2} ) \theta_2\ \end{aligned} $$

tvec <- seq(1,30,0.1)

pars <- c(
    beta = 0.5,
    gamma = 0.1,
    N = 2,
    i0 = 0.01
)

pars2 <- c(
    beta = 0.5,
    gamma = 0.100001,
    N = 2,
    i0 = 0.01
)

r1 <- SIR.detsim.hessian(tvec, pars)

r2 <- SIR.detsim.hessian(tvec, pars2)

matplot(cbind(r1$S, r2$S), type = "l")

plot((r2$S-r1$S)/1, type = "l")
lines(r1$nu_S_N, col = 2)
lines(r2$nu_S_N, col = 3)

nu_I_bg = (r2$nu_I_b-r1$nu_I_b)/0.000001

plot(nu_I_bg, type = "l")
lines(r1$nu_I_bg, col = 2)
lines(r2$nu_I_bg, col = 3)

True hessian function:

(hess.alg <- findHess(bombay2, fpars))

(hess.grad <- jacobian(findSens, fpars, data = bombay2, sensOnly = TRUE))

all.equal(hess.alg, hess.grad)

eigen(hess.alg)

I want to do this with negative log likelihood

grad(SIR.logLik(), fpars, count = bombay2$count)
findSens(bombay2, fpars, nll = TRUE, sensOnly = TRUE)

These should be identical... I'm guessing that we're not using the correct method for some reason...

$$ -l(x; \mu, \sigma^2) = \frac{1}{2} \log(\sigma^2) + \frac{(x-\mu)^2}{2\sigma^2} + C $$

Substitute $x = I$, $\mu = \hat{I}$, and $\sigma^2 = \sum_i \frac{(I_i - \hat{I}_i)^2}{n-1}$, and perform chain rule:

$$ \begin{aligned} \frac{d(nll)}{d\theta} &= \sum_i \frac{d}{d\theta} (\frac{1}{2} \log(\sigma^2) + \frac{(x_i-\mu_i)^2}{2\sigma^2})\ &= \sum_i (\frac{\mu_i-x_i}{\sigma^2} \frac{\partial \mu_i}{\partial \theta} + (\frac{1}{2\sigma^2} - \frac{(x_i-\mu_i)^2}{2 \sigma^4}) \frac{\partial \sigma^2}{\partial \theta})\ &= \sum_i (\frac{\mu_i-x_i}{\sigma^2} \frac{\partial \hat{I}_i}{\partial \theta} + (\frac{1}{2\sigma^2} - \frac{(x_i-\mu_i)^2}{2 \sigma^4}) \sum_j \frac{\partial \sigma^2}{\partial \hat{I}_i} \frac{\partial \hat{I}_i}{\partial \theta})\ \end{aligned} $$

nllSens <- function(data, params){
    t <- data$tvec
    tpars <- trans.pars(params)
    r <- SIR.detsim.hessian(t, tpars)

    with(as.list(c(r, tpars)),{

        n <- length(t)
        count <- data$count
        sigma2 <- sum((I-count)^2)/(n-1)

        sensVec <- c(beta, gamma, N, i0^2*exp(-qlogis(i0)))

        derVec <- c("nu_I_b", "nu_I_g", "nu_I_N", "nu_I_i")

        findDeriv <- function(i){
            d1 <- get(derVec[i])

            deriv <- sum((I-count)/sigma2 * d1 * sensVec[i] + (1/(2*sigma2) - ((I - count)^2)/(2*sigma2^2)) * sum(2 * (I - count)/(n-1) * d1 * sensVec[i]))
            return(deriv)
        }

        sens <- c(1:4)
        return(unlist(lapply(sens, findDeriv)))
    })
}

Does this work well?

##They certainly look similar...
grad(SIR.logLik(), fpars, count = bombay2$count)
nllSens(bombay2, fpars)

##Can I take a look at the hessian matrix?
hessian(SIR.logLik(), fpars, count = bombay2$count)
j <- jacobian(nllSens, fpars, data = bombay2)

fpars2 <- coef(f2 <- fitsir(bombay2, start = startfun(data = bombay2, auto = TRUE)))

##These look different now...
grad(SIR.logLik(), fpars2, count = bombay2$count)
nllSens(bombay2, fpars2)

fpars3 <- fpars2

fpars3[["logit.i"]] <- fpars2[["logit.i"]] + 1e-8

(SIR.logLik()(fpars3, bombay2$count)-SIR.logLik()(fpars2, bombay2$count))/1e-8

Can I find the hessian?

nllHess <- function(data, params){
    t <- data$tvec
    tpars <- trans.pars(params)
    r <- SIR.detsim.hessian(t, tpars)

    with(as.list(c(r, tpars)),{

        n <- length(t)
        count <- data$count
        sigma2 <- sum((I-count)^2)/(n-1)

        sensVec <- c(beta, gamma, N, i0^2*exp(-qlogis(i0)))

        sensVec2 <- c(1, 1, 1, 2*i0*exp(-qlogis(i0)) - i0^2 * exp(-qlogis(i0)) * 1/(i0-i0^2))

        derVec <- c("nu_I_b", "nu_I_g", "nu_I_N", "nu_I_i")

        derMat <- matrix(c("nu_I_bb", "nu_I_bg", "nu_I_bN", "nu_I_bi",
                           "nu_I_bg", "nu_I_gg", "nu_I_gN", "nu_I_gi",
                           "nu_I_bN", "nu_I_gN", "nu_I_NN", "nu_I_Ni",
                           "nu_I_bi", "nu_I_gi", "nu_I_Ni", "nu_I_ii"), 4, 4)

        findDeriv <- function(i1, i2){
            d1 <- get(derVec[i1])
            d2 <- get(derVec[i2])

            db <- get(derMat[i1,i2])

            if(i1 == i2){
                d12 <- sensVec2[i1]
            }else{
                d12 <- 0
            }

            diff_I <- I - count

            dsigma_I <- 2 * diff_I/(n-1)

            dsigma_theta1 <- sum(dsigma_I * d1)
            dsigma_theta2 <- sum(dsigma_I * d2)

            deriv <- sum((1/sigma2 * d2 - diff_I/sigma2^2 * dsigma_theta2) * d1 +
                diff_I/sigma2 * db +
                (-1/(2*sigma2^2) * dsigma_theta2 + (diff_I)^2/sigma2^3 * dsigma_theta2 - diff_I/(sigma2^2) * d2) * dsigma_theta1 +
                (1/(2*sigma2^2) - diff_I^2/(2*sigma2^2)) * sum(2 * (d2 * d1 + diff_I * db))) * sensVec[i1] * sensVec[i2] +
                sum(diff_I/sigma2 * d1 +
                        (1/(2*sigma2) - (diff_I^2)/(2*sigma2^2)) * dsigma_theta1) * d12 * sensVec[i2]

            return(deriv)
        }
        m <- matrix(NA, 4, 4)

        for(i in 1:4){
            for(j in 1:4){
                m[i,j] = findDeriv(i,j)
            }
        }

        return(m)
    })
}

Test:

jacobian(nllSens, fpars, data = bombay2)
nllHess(bombay2, fpars)

hessian(SIR.logLik(), fpars2, count = bombay2$count)
jacobian(nllSens, fpars2, data = bombay2)