Go back to Flu Mortality
model_SIR = function(t, theta) {
SIR = function(Time, State, Pars) { with(as.list(c(State, Pars)), { DecreaseS = -1 * beta * S * I / N DecreaseI = -1 * alpha * I dS = DecreaseS dI = -1 * DecreaseS + DecreaseI return(list(c(dS, dI))) }) } 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 = alpha2, N = N) out1 = ode(y = state1, times = t, func = SIR, parms = pars1, method = 'ode45') out2 = ode(y = state2, times = t, func = SIR, parms = pars2, method = 'ode45') return(cbind(-diff(out1[,"S"]), -diff(out2[,"S"]))) })
}
loglik_poisson = function(data, mean) {
sum(dpois(x=data, lambda=mean, log=TRUE))
}
exp_fit_given_logL = function(times, data, 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) } l = model(times, theta) #"base" is the vector of baseline values mean = l[,1] + l[,2] + base # compute the negative log-likelihood L = -logL(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='BFGS', start=as.list(theta0), control=list(maxit=1e5*length(theta0), 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) }
p = exp_fit_given_logL(c(0:19), data2, c(beta1 = 0.3, alpha1 = 0.2, N = 1377, I1 = 1, beta2 = 0.3, alpha2 = 0.2, I2 = 1), model_SIR, loglik_poisson)