Go back to Flu Mortality

library(bbmle)

library(deSolve)

data1 = c(106, 121, 136, 156, 183, 248, 279, 276, 239, 185, 160, 187, 150, 193, 153, 153, 138, 143, 119)

t1 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18)

base1 = c(105.1240, 106.3563, 107.5804, 108.7790, 109.9351, 111.0325, 112.0557, 112.9902, 113.8230, 114.5424, 115.1383, 115.6027, 115.9291, 116.1133, 116.1530, 116.0481, 115.8005, 115.4142, 114.8952)

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_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))

}

loglik_poisson = function(theta, data, mean) {

  sum(dpois(x=data, lambda=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)
      }
    
      l = model(times, theta)
    
      #"base" is the vector of baseline values
      mean = l[,1] + l[,2] + 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 theta0 contains no r,

  if (is.na(theta0['r']))
      theta0 = c(theta0, r = 10000)
    
  while (TRUE)
  {
      no_r = is.na(theta0['r']);
    
      if (!no_r)
      {
           p = SIR_fit_given_logL(times, data, base, theta0, model, loglik_negbin)
         
           print(p)
      
       }
         
      # if the fitted shape parameter r is too large (> 1000)
      # then it is approximately poisson
      # in this case we use poisson
    
      r1 = coef(p)["r"]
    
      if (no_r || r1[["r"]] > 1000) {
          # we need to remove r from theta
          if (!no_r) {
              old_theta0 = coef(p)
              theta0 = c()
              for (i in 1:length(old_theta0))
                  if (names(old_theta0[i]) != 'r')
                      theta0 = c(theta0, old_theta0[i])
          }
          p = SIR_fit_given_logL(times, data, base, theta0, model, loglik_poisson)
        
          print(p)
      }
    
  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))
  print(theta0)
  if (is.na(theta0['r']))
      theta0 = c(theta0, r = 10000)

  }

  return(list(fit=p, conf = conf, r=coef(p)["r"]))

}

c = SIR_fit(c(0:19), data1, base1, c(beta1 = 0.3, alpha1 = 0.2, N = 1377, I1 = 1, beta2 = 0.3, alpha2 = 0.2, I2 = 1), model_SIR)

print©

aicc = AICc(c"fit", nobs = length(data1))

print(aicc)