Exposed compartment jumps significantly in value

[bryce-carson]

The possible cause of #24 is this sudden jump in the exposed compartment; in the following table the jump occurs on the Eighth of March.

# A tibble: 10 × 15
   Date           N     S     V     E     I     R     D  newV  newE  newI  newR  newD  cumE  cumI
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-03-03 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 2 2020-03-04 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 3 2020-03-05 38137 38136     0     0     1     0     0     0     1     0     0     0     1     1
 4 2020-03-06 38137 38135     0     1     1     0     0     0     0     0     0     0     1     1
 5 2020-03-07 38137 38135     0     1     1     0     0     0     1     0     0     0     2     1
 6 2020-03-08 38137 38134     0     2     1     0     0     0    -1     1     0     0     2     2
 7 2020-03-09 38137 38134     0     1     2     0     0     0   193     0     0     0   195     2
 8 2020-03-10 38137 37940     0   194     2     0     0     0     0     0     0     0   195     2
 9 2020-03-11 38137 37940     0   194     2     0     0     0     0     0     0     0   195     2
10 2020-03-12 38137 37940     0   194     2     0     0    NA    NA    NA    NA    NA   195     2

SVEIRD.BayesianDataAssimilation(
       ## Parameters
       alpha = 0,
       beta = 0.5,
       gamma = 0.5,
       sigma = 0.001,
       delta = 0,
       lambda = 18,
       ## Model runtime
       n.days = simulation.days,
       ## Model data
       seedData = lieseed,
       neighbourhood.order = 0,
       populationSpatRaster = lierast,
       subregionsSpatVector = lievect,
       aggregationFactor = 0,
       startDate = "2020-03-03",
       countryCodeISO3C = "LIE",
       ## Model options
       simulationIsDeterministic = TRUE,
       dataAssimilationEnabled = FALSE,
       callback = callback)

[bryce-carson]

lievect <- getCountrySubregions.SpatVector("LIE", folder = "/tmp")
lierast <- getCountryPopulation.SpatRaster("LIE")

[bryce-carson]

seedLiechtensteinCOVID19.csv

[bryce-carson]

Transmission likelihoods

This is a counter-example of the edge effect... the transmission likelihoods does not exhibit this.

[bryce-carson]

On the eighth of March, the newExposed is calculated as the following

Browse[2]> newExposed
class       : SpatRaster 
dimensions  : 26, 20, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 9.47375, 9.640417, 47.05792, 47.27458  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
name        : Susceptible 
min value   :    0.000000 
max value   :    1.058809 
Browse[2]> terra::plot(newExposed)
Browse[2]> savePlot("tests/testthat/Liechtenstein/newExposedOverOne.png")
Browse[2]> 

I believe the issue lies with the maximum value being over one, but it could be the values between these (because they're everywhere).

[bryce-carson]

The values of growth range from 0.46 to 1.11. Clearly growth should never be more than

max value   :   0 
Browse[2]> growth
class       : SpatRaster 
dimensions  : 26, 20, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 9.47375, 9.640417, 47.05792, 47.27458  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source(s)   : memory
name        : Susceptible 
min value   :   0.4569843 
max value   :   1.1084694 
Browse[2]> terra::plot(growth)
Browse[2]> savePlot("tests/testthat/Liechtenstein/lieGrowth.png")
Browse[2]> 

The only explanation I can think of is line 1054 (as quoted) leads to the exposed being permitted much beyond where they should be, where the global sum is then inflated (because it's including exposure probabilities less than one) when it should only be a few. It's confusing because $\beta$ is what is actually the probability of exposure, so I'm not sure why the old code had these checks and re-assignment, but that's what I followed in the new code.

Thoughts, @ashokkrish?

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1053 to 1056 in 6dd5f47

	newExposed <- terra::mask(if(simulationIsDeterministic) growth else stats::rpois(1, growth),
	c(proportionSusceptible < 1, transmissionLikelihoods < 1),
	maskvalues = TRUE,
	updatevalue = 0)

[bryce-carson]

Changing the previously quoted lines to the following results in a drastic change in the results.

spatialEpisim.foundation/R/SVEIRD.BayesianDataAssimilation.R

Lines 1053 to 1056 in 6dd5f47

	newExposed <- terra::mask(if(simulationIsDeterministic) growth else stats::rpois(1, growth),
	c(proportionSusceptible < 1, transmissionLikelihoods < 1),
	maskvalues = TRUE,
	updatevalue = 0)

      newExposed <- terra::mask(if(simulationIsDeterministic) growth else stats::rpois(1, growth),
                                proportionSusceptible < 1,
                                maskvalues = TRUE,
                                updatevalue = 0) %>%
        terra::mask(transmissionLikelihoods < 1,
                    maskvalues = TRUE,
                    updatevalue = 0)

> SVEIRD.BayesianDataAssimilation(
       ## Parameters
       alpha = 0,
       beta = 0.5,
       gamma = 0.5,
       sigma = 0.001,
       delta = 0,
       lambda = 18,
       ## Model runtime
       n.days = simulation.days,
       ## Model data
       seedData = lieseed,
       neighbourhood.order = 0,
       populationSpatRaster = lierast,
       subregionsSpatVector = lievect,
       aggregationFactor = 0,
       startDate = "2020-03-03",
       countryCodeISO3C = "LIE",
       ## Model options
       simulationIsDeterministic = TRUE,
       dataAssimilationEnabled = FALSE,
       callback = callback)
$table
# A tibble: 10 × 15
   Date           N     S     V     E     I     R     D  newV  newE  newI  newR  newD  cumE  cumI
   <date>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1 2020-03-03 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 2 2020-03-04 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 3 2020-03-05 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 4 2020-03-06 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 5 2020-03-07 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 6 2020-03-08 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 7 2020-03-09 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 8 2020-03-10 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
 9 2020-03-11 38137 38136     0     0     1     0     0     0     0     0     0     0     0     1
10 2020-03-12 38137 38136     0     0     1     0     0    NA    NA    NA    NA    NA     0     1

[bryce-carson]

@ashokkrish, I've identified that I made a small, understandable mistake when refactoring. I used the proportion susceptible (pSusceptible in the old code, proportionSusceptible in the new code), not the susceptible layer from the previous timestep (i.e. valSusceptible and layers$Susceptible) when I was sub-assigning newExposed.

I was going to ask you "What is the significance of sub-assigning zero wherever the pSusceptible and I_tilda ( $\tilde{I}$) are less than one?," but I have realized this small mistake (in the previous paragraph) and that drastically changes the meaning to simply be the following:

"Wherever the likelihood of transmission from an infected person to a susceptible person is less than one, no exposure occurs, and wherever there are no susceptible people no new exposures can occur there."

That's a drastic change in meaning, and given I've observed the results change for the better due to this and the meaning of the code has been made clearer, I'm going to push these changes once I verify that this issue and issue #24 have been fixed.

    if(deterministic) {
      newExposed <- beta*pSusceptible*I_tilda
    } else {
      rpois(1, beta*pSusceptible*I_tilda)
    }
    newExposed[valSusceptible < 1] <- 0
    newExposed[I_tilda < 1] <- 0

[bryce-carson]

Have you had a chance to review this, or other issues, @ashokkrish ?