In fitting ODE model to data by maximum likelihood estimation, why fit transformed values of the actual parameters?

Question

I'm looking at some published R code in which the authors are fitting an ODE model to some data using maximum likelihood estimation. I don't understand why they are fitting transformed values of the parameters - and I hope to use this code for an exercise in a class, so I have to be able to explain why the code is written this way.

Here's the code:

# Necessary libraries
library(stats4)
library(deSolve)
library(chron)
library(bbmle)

# Read the data
ebola <- read.csv("data_WHO_SierraLeone_mle.csv",skip=3)

# Select only incidence counts between these two time points
begin <- chron("25 May 2014", format=c(dates = "day mon year"))
end <- chron("13 Sep 2015", format=c(dates = "day mon year"))
timepoints <- seq(begin,end,7)
ebola <- as.data.frame(cbind(timepoints,ebola[21:89,2:5]))

# Add last two entries of situation report to patient data base
ebola[68:69,4:5] <- ebola[68:69,2:3]
ebola <- ebola[,-2:-3]
ebola[,2] <- ebola[,2]+ebola[,3]
ebola[,3] <- cumsum(ebola[,2])
names(ebola) <- c("Date","Cases","Cumulative")

# Definition of the SEIR model
SEIR <- function(t, x, parms) {
    with(as.list(c(parms,x)),{
        # Exponential reduction in transmission rate
        ifelse(t < tau1, beta <- beta0, beta <- beta1 + (beta0-beta1)*exp(-k*(t-tau1)))
        # Density-dependent transmission
        N <- S + E + I + R + D
        dS <- - beta/N*S*I
        dE <- beta/N*S*I - sigma*E
        dI <- sigma*E - gamma*I
        dR <- (1-f)*gamma*I
        dD <- f*gamma*I
        dC <- sigma*E
        der <- c(dS,dE,dI,dR,dD,dC)
        list(der)
    })
}

# Negative log-likelihood
nll <- function(beta0,beta1,k,f,tau0,tau1,sigma,gamma,disp) {
    pars <- c(beta0=beta0,beta1=beta1,k=k,f=f,tau0=tau0,tau1=tau1,sigma=sigma,gamma=gamma,disp=disp)
    pars <- trans(pars)
    times <- c(0,min(data$times)+pars["tau0"]-7,data$times+pars["tau0"])
    simulation <- as.data.frame(ode(init,times,SEIR,parms=pars))
    ll <- sum(dnbinom(data$cases,size=pars["disp"]*diff(simulation$C)[-1],mu=diff(simulation$C)[-1],log=TRUE))
    return(-ll)
}

# Parameter transformation
trans <- function(pars) {
    pars["beta0"] <- exp(pars["beta0"])
    pars["beta1"] <- exp(pars["beta1"])
    pars["k"] <- exp(pars["k"])
    pars["f"] <- plogis(pars["f"])
    pars["tau0"] <- exp(pars["tau0"])
    pars["tau1"] <- exp(pars["tau1"])
    pars["disp"] <- exp(pars["disp"])
    return(pars)
}

# Prepare the data and set the initial values
data <- na.omit(ebola[c("Date","Cases")])
names(data) <- c("times","cases")
data$times <- data$times - data$times[1]
N <- 6.316e6      
init <- c(S = N - 1, E = 0, I = 1, R = 0, D = 0, C = 0)

# Fit the model to the data
fixed <- c(sigma = 1/11.4, gamma = 1/(15.3-11.4), f = qlogis(0.69), tau0 = log(32))
free <- c(beta0 = log(2.02/(15.3-11.4)), beta1 = log(0.5/(15.3-11.4)), tau1 = log(32), k = log(0.01), disp = log(0.18))
fit <- mle2(nll,start=as.list(free),fixed=as.list(fixed),method="Nelder-Mead",control=list(maxit=1e3))
summary(fit)
trans(coef(fit))

My expectation would have been that you would simply put in values for the free parameters (beta0, beta1, tau1, k, disp) and calculate the negative log likelihood, rather that putting in transformed parameters, un-transforming them, and then calculating the NLL. Not sure what difference this makes.

Attribution: https://doi.org/10.1371/journal.pntd.0004676 (S2 appendix contains the data that the model is being fit to)

Answer 1

By log transformation, one can ensure positivity of parameters. The optimizer can then work on the log-transformed scale within the range between (-Inf, Inf), while the parameters of the model are in the range (0, Inf).

Similarly to this, the plogis-transformation helps to keep f within the interval (0, 1).

This helps to stabilize optimization as it avoids unrealistic parameter values and potential numerical errors.

x <- seq(-5, 1, 0.01)
plot(x, exp(x), type="l", xlab="optimizer", ylab="ode model", main="exp")

x <- seq(-5, 5, 0.1)
plot(x, plogis(x), type="l", xlab="optimizer", ylab="ode model", main="sigmoidal")