Go back to Flu Mortality
library(bbmle)
library(deSolve)
dyn.load(“mymod.so”)
data8 = c(117, 121, 109, 132, 145, 130, 139, 115, 144, 131, 150, 145, 134, 171, 198, 234, 378, 373, 319, 234, 219, 203, 193, 188, 218, 176, 218, 181, 193, 190, 180, 167, 154, 154, 158, 122, 122, 136)
t8 = c(0:37)
base8 = c(106.3900, 106.8131, 107.3700, 108.0529, 108.8525, 109.7575, 110.7553, 111.8318, 112.9719, 114.1594, 115.3776, 116.6092, 117.8369, 119.0433, 120.2112, 121.3243, 122.3668, 123.3241, 124.1826, 124.9304, 125.5571, 126.0540, 126.4143, 126.6332, 126.7081, 126.6381, 126.4247, 126.0715, 125.5840, 124.9697, 124.2378, 123.3994, 122.4671, 121.4548, 120.3775, 119.2513, 118.0929, 116.9195)
model_SIR = function(t, theta) {
with(as.list(theta), { state1 = c(S = N - I1, I = I1) pars1 = c(beta = beta1, alpha = alpha1, N = N) state2 = c(S = N - I2, I = I2) pars2 = c(beta = beta2, alpha = alpha1, N = N) state3 = c(S = N - I3, I = I3) pars3 = c(beta = beta3, alpha = alpha1, N = N) out1 = ode(y = state1, times = t, func = "derivs", parms = pars1, jacfunc = "jac", dllname = "mymod", initfunc = "initmod", nout = 1, outnames = "S") out2 = ode(y = state2, times = t, func = "derivs", parms = pars2, jacfunc = "jac", dllname = "mymod", initfunc = "initmod", nout = 1, outnames = "S") out3 = ode(y = state3, times = t, func = "derivs", parms = pars3, jacfunc = "jac", dllname = "mymod", initfunc = "initmod", nout = 1, outnames = "S") return(cbind(-diff(out1[,"S"]), -diff(out2[,"S"]), -diff(out3[,"S"]))) })
}
loglik_negbin = function(theta, data, mean) {
# for negative bimonial, mean = r p /(1-p) r = theta[["r"]]; # if r < 0, return a tiny likelihood if (r < 0) { return(-1e10) } sum(dnbinom(x=data, size=r, mu=mean, log=TRUE))
}
SIR_fit_given_logL = function(times, data, base, theta0, model, logL) { # construct a general negative log-likelyhood function
L <- function() { # reconstruct the vector of named pairs of # parameters from the list of arguments l=length(theta0) theta=c() vars=names(theta0) for (i in 1:l) { item=c(x=get(vars[i])) names(item) <- vars[i] theta <-c(theta, item) } N = theta[["N"]] if(N < sum(data)) { return(1e10) } for(i in 1:length(theta)) { if(theta[[i]] < 0) { return(1e10) } } l = model(times, theta) #"base" is the vector of baseline values mean = l[,1] + l[,2] + l[,3] + base # compute the negative log-likelihood L = -logL(theta, data, mean) # if not a number, return a very small likelihood if (is.na(L)) { return(1e10) } return(L) }
# replace the input arguments of L by the list of parameters formals(L)←as.list(theta0) # mle2 # sometimes confint may find a better solution. # In this case we need to refit the model with the better # solution as a starting paoint, because the parameter # names in the returned better fit can be wrong fit = mle2(L, method = “BFGS”, start = as.list(theta0), control = list(maxit = 1e6*length(theta0))) # fit = mle2(L, method = “Nelder-Mead”, start = as.list(theta0), control = list(maxit = 1e6*length(theta0), ndeps = 1e-4, reltol = 1e-10)) # if there is only one parameter to be fitted # the returned confidence interval is a vector # change it to a matrix return(fit) }
SIR_fit = function(times, data, base, theta0, model) {
# we start with negative binomial if (is.na(theta0['r'])) { theta0 = c(theta0, r = 128) } while (TRUE) { p = SIR_fit_given_logL(times, data, base, theta0, model, loglik_negbin) conf = confint(p) if(is.vector(conf)) { conf = matrix(conf, nrow = 1) } if(is.matrix(conf)) { break } theta0 = coef(conf) names(theta0) = names(coef(p)) } return(list(fit=p, conf = conf, r=coef(p)["r"]))
}
c = SIR_fit(c(0:38), data8, base8, c(beta1 = 7.4, alpha1 = 7, N = (sum(data8)*1.166), I1 = 1, beta2 = 7.6, I2 = 1, beta3 = 7.01, I3 = 1), model_SIR)
print©
aicc = AICc(c"fit", nobs = length(data8))
print(aicc)