Summer Project: Flu Mortality

Ten Years Graph

baseline_plot.pdf

First Five Years Graph

baseline_plot_1st_5yr.pdf

Second Five Years Graph

baseline_plot_2nd_5yr.pdf

Best Model Results

Best Results Table

1. Create a model that adds 2 SIR epidemic curves with the fitted baseline, and fit to the P&I mortality curve, instead of the excess mortality curve.
2. Test the new code using mle2, and compare its results to the old code
3. Test a model with 3 epidemics, and compare these models with AICc
4. On the cluster, run a sequence of jobs for each season, with 2 different optimization methods ('BFGS' and 'Nelder-Mead'), for both the model with 2 SIR epidemic with different alphas and the one with 2 SIR epidemics and a single alpha value, and models with 3 SIR epidemics
   • For Nelder-Mead, try maxit=1e5 first, if it is worse than BFGS or did not converge, try maxit=1e6 or larger, and decrease the reltol
   • For the 2-SIR+baseline model with 2 alpha parameters, compare the alphas across seasons.

Method

Code using bbmle package

Test R code for a single SIR epidemic

Test results for a single SIR epidemic

Test results for a single SIR epidemic (with MLE2 method)

Single SIR model for 1989-1990 epidemic

Single SIR model for 1989-1990 epidemic (with MLE2 method)

Weston's new R code for MLE 2 method

Double SIR Epidemic Model

Weston's R code for the two SIR model

Weston's R code for the two SIR model (with MLE2 method)

Weston's R code for the two SIR model (with MLE2 method and Negative Binomial/Poisson log-likelihood functions)

Triple SIR Epidemic Model

Weston's R code for the three SIR model (with MLE2 method)

Weston's R code for the three SIR model (with MLE2 method and Negative Binomial/Poisson log-likelihood functions)

Double SIR model for 1989-1990 epidemic with "BFGS" method (different alpha values)

Double SIR model for 1989-1990 epidemic with "SANN" method (different alpha values)

Double SIR model for 1989-1990 epidemic with "Nelder-Mead" method (different alpha values)

Double SIR model for 1989-1990 epidemic with "BFGS" method (same alpha values)

AICc to compare the "BFGS" method model with the same and different alpha values

Double SIR model for 1989-1990 epidemic with "BFGS" method (with alpha values equal to 1/7) (Junling: Note that alpha=1/7 means the infectious period is 7 weeks, not 7 days, which means you are not likely to get reasonable fit, as illustrated by your results)

From the MLE2 Method

Double SIR model for 1989-1990 epidemic with "BFGS" method (different alpha values, mle2)

Double SIR model for 1989-1990 epidemic with "BFGS" method (different alpha values, mle2, Neg Bin/Poi log-lik)

Three SIR model for 1989-1990 epidemic with "BFGS" method (different alpha values, mle2)

Three SIR model for 1989-1990 epidemic with "BFGS" method (different alpha values, mle2, Neg Bin/Poi log-lik)

Double SIR model for 1992-1993 epidemic with "BFGS" method (different alpha values)

Double SIR model for 1992-1993 epidemic with "BFGS" method (same alpha values)

First eight of the models (1990-1991 same except different data)

1989-1990 Epidemic 2 SIR model baseline Different Alpha BFGS

1989-1990 Epidemic 2 SIR model baseline Different Alpha NM

1989-1990 Epidemic 2 SIR model baseline Same Alpha BFGS

1989-1990 Epidemic 2 SIR model baseline Same Alpha NM

1989-1990 Epidemic 3 SIR model baseline Different Alpha BFGS

1989-1990 Epidemic 3 SIR model baseline Different Alpha NM

1989-1990 Epidemic 3 SIR model baseline Same Alpha BFGS

1989-1990 Epidemic 3 SIR model baseline Same Alpha NM

1989-1990 Epidemic 2 SIR model baseline Different Alpha BFGS (r = 1)

Important Data from 1989-1999

This code restricts all of the parameters to be positive, as well as having N be larger than the total, otherwise it will return a small likelihood (1e10).

I tried using 1/r as the parameter and the mle2 optimizer was not very happy with it, many times giving the response

Error in .local(fitted, …) :

Hessian is ill-behaved or missing,can't find an initial estimate of std. error(consider specifying std.err in profile call)

Error in .local(object, parm, level, …) :

Problem with profiling: Error in .local(fitted, ...) : 
Hessian is ill-behaved or missing,can't find an initial estimate of std. error(consider specifying std.err in profile call)

Calls: SIR_fit → confint → confint → .local Execution halted

1989-1990 Epidemic 2 SIR model baseline Different Alpha BFGS Jul 9 (worked)

1989-1990 Epidemic 2 SIR model baseline Different Alpha NM Jul 9 (did not work, “Execution halted”)

1989-1990 Epidemic 2 SIR model baseline Same Alpha BFGS Jul 9 (did not work, “Execution halted”)

1989-1990 Epidemic 2 SIR model baseline Same Alpha NM Jul 9 (did not work, “Execution halted”)

1989-1990 Epidemic 3 SIR model baseline Different Alpha BFGS Jul 9 (did not work, “Execution halted”)

1989-1990 Epidemic 3 SIR model baseline Different Alpha NM Jul 9 (did not work, “Execution halted”)

1989-1990 Epidemic 3 SIR model baseline Same Alpha BFGS Jul 9 (did not work, “Execution halted”)

1989-1990 Epidemic 3 SIR model baseline Same Alpha NM Jul 9 (did not work, “Execution halted”)

I think all of these other models did not work because I changed the starting values to reflect the good values that seemed to be working with the

“1989-1990 Epidemic 2 SIR model baseline Different Alpha BFGS Jul 9”.

Additionally, I started the parameter r at for all the others 760, whereas I started r at 2 for the model that worked. I think it's probably better to let the optimizer “find” the value of r, than make it so high initially.

SIR equations written in C code

SIR Equations in C

Test code in R to make sure the SIR equations in C returned the same values as if they were written in R.

Test code for SIR Equations in C

The code is correct as evidenced by the results of the test code.

Results of the Test Code for SIR Equations in C

Method for Even Faster C Code (ode solver also in C using gsl)

SIR Equations and ODE Solver in C

Test code for SIR Equations and ODE Solver in C

Unfortunately, the above test code does not work. The same error is repeated

Error in dyn.load(“SIR_solver.so”) :

unable to load shared library '/home1e/wcr/SIR_solver.so':
/home1e/wcr/SIR_solver.so: undefined symbol: gsl_odeiv_step_rkf45

Execution halted

It appears that other people have had this problem too with linking the gls library from C code to R code.

https://stat.ethz.ch/pipermail/r-devel/2006-December/043927.html

http://www.mail-archive.com/r-help@r-project.org/msg76240.html

The PDF below for linking is good but complicated and in C++ rather than C. http://www.google.ca/#hl=en&source=hp&q=%24+PKG_LIBS%3D%22-lgsl+-lblas%22&btnG=Google+Search&aq=f&aqi=&aql=&oq=%24+PKG_LIBS%3D%22-lgsl+-lblas%22&gs_rfai=&fp=77a98ce8ed84c158

I tried

$ PKG_LIBS = “-L/user/local/lib/libgsl -lgsl” R CMD SHLIB SIR_solver.c

on the command line and nothing happened. I know that gsl is installed on the “euler”.

Thanks for your help! Wes

SIR equations and Jacobian matrix in C

Result with Poisson distribution after r > 10,000:

Result from 1989-1990 epidemic 2 SIR model Baseline Different Alpha BFGS Jul 9

Code for the model with only the Negative Binomial option:

Code for 1989-1990 Epidemic 2 SIR model Baseline Different Alpha BFGS Neg Bin

Result with only Negative Binomial option:

Result 1989-1990 Epidemic 2 SIR model Baseline Different Alpha BFGS Neg Bin

July 16th

Code for the model using C code with Poisson distribution after r > 10,000 (same code as Jul 9th model with new faster C code):

Code for 1989-1990 Epidemic 2 SIR model Baseline Different Alpha BFGS (poi/neg C code)

Result using C code with Poisson distribution after r > 10,000:

Result 1989-1990 epidemic 2 SIR model Baseline Different Alpha BFGS Jul (poi/neg C code)

July 19th (very quick to finish 3 and 1/2 hours versus 3-4 days with July 16th code)

Code for the model using lsoda, Jacobian, and C code with Poisson distribution after r > 1,000,000:

Code for 1989-90 Epidemic 2 SIR model Baseline Different Alpha BFGS (poi/neg fast C code)

Result using very fast C code with Poisson distribution after r > 1,000,000:

Result 1989-1990 epidemic 2 SIR model Baseline Different Alpha BFGS Jul (poi/neg fast C code)

Code for the model using lsoda, Jacobian, and C code with only Neg. Bin. distribution:

Code for 1989-90 Epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result using very fast C code with only Neg. Bin. distribution:

Result 1989-1990 epidemic 2 SIR model Baseline Different Alpha BFGS Jul (neg fast C code)

AICc value 175.1937

Best so far

July 21st

Code for 1989-90 Epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1989-1990 epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 169.8512

July 28th

NA value(s) in the CI

Code for 1989-90 Epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1989-1990 epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 182.0587

Best so far

July 21st

Code for 1990-91 Epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1990-91 epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 190.7815

July 26th

Code for 1990-91 Epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1990-91 epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 210.1084

Best so far

July 21st

Code for 1991-92 Epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1991-92 epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 189.7440

July 28th

NA value(s) in CI

Code for 1991-92 Epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1991-92 epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 197.7814

July 28th

NA value(s) in CI

Code for 1991-92 epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1991-92 epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 196.0226

Best so far

July 22nd

NA value(s) in the CI

Code for 1992-93 Epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1992-93 epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 202.8748

July 22nd

Code for 1992-93 Epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1992-93 epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 207.8083

July 22nd

Code for 1995-96 Epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1995-96 epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 275.9575

July 28th

NA value(s) in the CI

Code for 1995-96 epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1995-96 epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 275.9059

Best so far

July 25th

Code for 1995-96 Epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1995-96 epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 270.8821

July 22nd

Code for 1996-97 Epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1996-97 epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 352.8229

Best so far

July 25th

Code for 1996-97 Epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1996-97 epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 342.7335

July 22nd

Code for 1997-98 Epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1997-98 epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 327.6377

Best

July 28th

NA value(s) in the CI

Code for 1997-98 Epidemic 2 SIR model Baseline Same Alpha SANN/BFGS (neg fast C code)

Result 1997-98 Epidemic 2 SIR model Baseline Same Alpha SANN/BFGS (neg fast C code)

AICc value 325.0415

July 23rd

Code for 1997-98 Epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1997-98 epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 326.4857

July 25th

Code for 1997-98 Epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1997-98 epidemic 3 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 330.13

July 23rd

Code for 1998-99 Epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1998-99 epidemic 2 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 308.5382

Best so far

July 30th

NA value(s) in the CI

Code for 1998-99 Epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

Result 1998-99 epidemic 2 SIR model Baseline Same Alpha BFGS (neg fast C code)

AICc value 307.3100

July 23rd

Code for 1998-99 Epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

Result 1998-99 epidemic 3 SIR model Baseline Different Alpha BFGS (neg fast C code)

AICc value 323.2764