dynmod: Dynamic models

##' Estimate immunity from MMR and case data
##'
##' @param countries the countries to investigate. If set to NULL
##' (default), will estimate immunity for all countries
##' @param vaccine.efficacy vaccine efficacy
##' @param estimation.year year for which to estimate the profile
##' @param age.dist.cases age distribution of cases
##' @param reporting proportion of cases that get reported
##' @param average.last average over the last X years to get further
##' years of cases (0 if no averaging is to be done)
##' @import reshape2
##' @return A data table with immunity by age in 2014
##' @author Sebastian Funk
##' @export
estimateMMRImmunity <- function(countries = NULL,
                                vaccine.efficacy = 0.95,
                                estimation.year = 2014,
                                age.dist.cases = NULL,
                                reporting = 1,
                                average.last = 0,
                                coverage.estimate = FALSE,
                                use.sia = FALSE) {

    # load data
    data(europe_mmr_timings)
    if (coverage.estimate)
    {
        data(mcv_world_estimate)
        mcv.world <- mcv.world.estimate
    } else
    {
        data(mcv_world)
    }

    # if countries aren't specified, estimate immunity for all
    # countries for which there is data
    if (is.null(countries)) {
        countries <- intersect(mcv.world[, country], mmr.timings[, country])
    }

    # merge vaccination and timing data
    setkey(mcv.world, country, year)
    setkey(mmr.timings, country)
    mmr.timings[, vaccine := NULL]
    mmr.immunity <-
        merge(mcv.world[country %in% countries],
              mmr.timings[country %in% countries],
              by = "country")

    # set variables
    mmr.immunity <- mmr.immunity[, year := as.integer(year)]
    mmr.immunity <- mmr.immunity[vaccine == "MCV", mcv.min := mcv1.min]
    mmr.immunity <- mmr.immunity[vaccine == "MCV", mcv.max := mcv1.max]
    mmr.immunity <- mmr.immunity[vaccine == "MCV2", mcv.min := mcv2.min]
    mmr.immunity <- mmr.immunity[vaccine == "MCV2", mcv.max := mcv2.max]
    mmr.immunity <- mmr.immunity[year <= estimation.year]
    max.year <- mmr.immunity[, ceiling(mcv.max)]
    # age relative to estimation year
    mmr.immunity <-
        mmr.immunity[, age := estimation.year - year + max.year]
    mmr.immunity <- mmr.immunity[age >= 0]

    # estimate immunity that is generated by each vaccine dose
    mmr.immunity$immunity <- apply(mmr.immunity, 1, function(x) {

        age <- as.numeric(x[["age"]])
        mcv.min <- as.numeric(x[["mcv.min"]])
        mcv.max <- as.numeric(x[["mcv.max"]])
        uptake <- as.numeric(x[["uptake"]])

        min(1, 1 - min(0, age - mcv.min) /
                ifelse(mcv.max > max(age + 1),
                       mcv.max - mcv.min, 1)) * uptake
    })

    # wide table (for combining 1st and 2nd dose)
    mmrd <- data.table(dcast(mmr.immunity, country + age ~ vaccine,
                             value.var = "immunity"))
    mmrd[, age := as.integer(as.character(age))]
    if ("MCV2" %in% colnames(mmrd))
    {
        mmrd <- mmrd[is.na(MCV2), MCV2 := 0]
    } else
    {
        mmrd <- mmrd[, MCV2 := 0]
    }
    # combine first and second dose for vaccinal.immunity
    mmrd <- mmrd[, vaccinal.immunity := MCV * (vaccine.efficacy +
                                                   (1 - vaccine.efficacy) *
                                                       MCV2/100 * vaccine.efficacy)]
    setkey(mmrd, country)

    if (use.sia || !is.null(age.dist.cases))
    {
        data(pop_world_age)
        data(ms_world)

        # world population only available in 5-year age groups, so
        # let's group our immunity data in this way

        max.age <- ceiling(max(mmrd[, age]) / 5) * 5
        mmrd[, agegroup := cut(x = mmrd[, age],
                        breaks = c(0, seq(5, max.age + 5, 5)), right = FALSE)]
        mmrd.agegroups <- mmrd[, list(lowest.age = min(age),
                                      highest.age = max(age),
                                      vaccinal.immunity = mean(vaccinal.immunity),
                                      MCV1 = mean(MCV),
                                      MCV2 = mean(MCV2),
                                      age.part = (max(age) - min(age) + 1) / 5),
                               by = list(country, agegroup)]
        mmrd.agegroups[, lower.age.limit := lowest.age %/% 5 * 5]

        mmrd.agegroups <-
            merge(mmrd.agegroups,
                  pop.world.age[, list(country, lower.age.limit,
                                       both.sexes.population)],
                  by = c("country", "lower.age.limit"))
        setnames(mmrd.agegroups, "both.sexes.population", "population")
        mmrd.agegroups[, immunity := vaccinal.immunity]

        # do we want to take into account historical case data?
        if (!is.null(age.dist.cases))
        {
            # now deal with cases
            mmrd.cases <- apply(mmrd.agegroups, 1, function(x)
            {
                country.cases <- ms.world[country == x[["country"]]]
                country.cases <- country.cases[!is.na(cases)]
                country.cases[, years.passed := estimation.year - year]
                # fill cases with average of last few years
                if (average.last > 0)
                {
                    max.passed <- max(country.cases[, years.passed])
                    if (max.passed < as.integer(x[["highest.age"]]))
                    {
                        last.year <- min(country.cases[, year])
                        avg.cases <-
                            mean(sapply(seq(last.year, last.year + average.last - 1),
                                        function(current.year)
                                        {
                                            if (current.year %in% country.cases[, year])
                                            {
                                                country.cases[year == current.year, cases]
                                            } else {
                                                0
                                            }
                                        }), na.rm = T)
                        for (add.year in seq(max.passed, as.integer(x[["highest.age"]])))
                        {
                            # copy last row
                            country.cases <- rbind(country.cases,
                                                   country.cases[nrow(country.cases)])
                            # change year and cases in new last row
                            country.cases[nrow(country.cases), year := year + 1]
                            country.cases[nrow(country.cases), cases := avg.cases]
                            country.cases[nrow(country.cases),
                                          years.passed := years.passed + 1]
                        }
                    }
                }

                country.cases[, true.cases := cases / reporting]
                ages.in.group <- seq(as.integer(x[["lowest.age"]]),
                                     as.integer(x[["highest.age"]]))
                country.cases <-
                    country.cases[years.passed <= as.integer(x[["highest.age"]])]

                if (nrow(country.cases) > 0)
                {
                    agegroup.cases <- sum(apply(country.cases, 1, function(y)
                    {
                        age.group.cases <-
                            ages.in.group - as.integer(y[["years.passed"]]) + 1
                        age.group.cases <- age.group.cases[age.group.cases >= 1]
                        round(sum(age.dist.cases[age.group.cases]) *
                                  as.integer(y[["true.cases"]]))
                    }))
                    agegroup.cases / as.integer(x[["population"]]) * 100
                } else {
                    0
                }
            })

            mmrd.agegroups[, natural.immunity := mmrd.cases]
            mmrd.agegroups[, immunity := min(100, immunity + natural.immunity),
                           by = 1:nrow(mmrd.agegroups)]
        }

        # do we want to take into account SIAs?
        if (use.sia)
        {
            ## initialise
            mmrd.agegroups[, sia.immunity := 0]
            mmrd[, sia.immunity := 0]
            data(sia)

            ## set an oldest age (so the highest age group has an appropriate upper bound
            oldest.age <- 123
            ## only take into account SIAs which have an upper and lower age bound
            sia <- sia[!(is.na(age.lower) | is.na(age.upper))]

            sia <- sia[country %in% countries]

            ## loop over SIAs
            for (row in seq_len(nrow(sia)))
            {
                ## get country's 2013 age-based population figures
                ## country.pop <-
                ##     pop.world.age[as.character(country) == sia[row, country]]
                ## ## set upper limits of age groups in population figures
                ## for (country_row in seq_len(nrow(country.pop) - 1))
                ## {
                ##     country.pop <-
                ##         country.pop[country_row,
                ##                     upper.age.limit :=
                ##                         country.pop[country_row + 1, lower.age.limit] - 1]
                ## }
                ## ## set upper limit of oldest ageg roup
                ## country.pop[nrow(country.pop), upper.age.limit := oldest.age]

                ## get age range (at the time of estimation) of targeted population
                age.range <- unlist(sia[row, list(age.lower, age.upper)]) +
                    estimation.year - sia[row, year]
                ## identify age groups in population data that correspond
                ## to age range of targeted population
                youngest.vacc.age <-
                    mmrd.agegroups[as.character(country) == sia[row, country],
                                   min(lowest.age)]
                oldest.vacc.age <-
                    mmrd.agegroups[as.character(country) == sia[row, country],
                                   lower.vacc.group <-
                                       max(highest.age)]

                if (youngest.vacc.age <= age.range[1])
                {
                    lower.vacc.group <-
                        mmrd.agegroups[as.character(country) ==
                                           sia[row, country] &
                                               lowest.age <= age.range[1],
                                       max(lowest.age)]
                } else
                {
                    lower.vacc.group <- youngest.vacc.age
                }
                if (oldest.vacc.age >= age.range[2])
                {
                    higher.vacc.group <-
                        mmrd.agegroups[as.character(country) ==
                                           sia[row, country] &
                                               highest.age >= age.range[2],
                                       min(highest.age)]
                } else
                {
                    higher.vacc.group <- oldest.vacc.age
                }

                ## lower.group <-
                ##     country.pop[lower.age.limit <= age.range[1], max(lower.age.limit)]
                ## higher.group <-
                ##     country.pop[lower.age.limit >= age.range[2], min(upper.age.limit)]

                ## get population size of target group
                ## pop <- country.pop[lower.age.limit >= lower.group &
                ##                        upper.age.limit <= higher.group,
                ##                    sum(both.sexes.population)]

                pop <- mmrd.agegroups[as.character(country) == sia[row, country] &
                                          lowest.age >= lower.vacc.group &
                                              highest.age <= higher.vacc.group,
                                      sum(population)]

                ## get proportion of the population in the target age
                ## groups that has been targeted
                proportion.targeted <- min(sia[row, target] / pop, 1, na.rm = TRUE)

                ## get proportion of the population targeted that has been reached
                proportion.reached <- sia[row, pcent.reached / 100]
                if (is.na(proportion.reached))
                {
                    ## if there is no information, assume that no one has been reached
                    proportion.reached <- max(sia[row, reached] / pop, 0, na.rm = TRUE)
                }

                mmrd.agegroups[as.character(country) == sia[row, country] &
                                   lowest.age >= lower.vacc.group &
                                       highest.age <= higher.vacc.group,
                               sia.immunity := sia.immunity +
                                   (100 - immunity - sia.immunity) *
                                       proportion.targeted *
                                           proportion.reached]
            }
            mmrd.agegroups[, immunity := immunity + sia.immunity]
        }

        mmrd.agegroups[, susceptible.population :=
                           round((100 - immunity) / 100 * population *
                                     (highest.age - lowest.age + 1) / 5)]
        return(mmrd.agegroups)
    } else {
        return(mmrd)
    }
}