[Fix] [S-R trend analysis] Italy shows weird spike in estimated cases

[Inglezos]

Summary

I run the estimator for Italy:

import covsirphy as cs

data_loader = cs.DataLoader(directory="kaggle/input")
jhu_data = data_loader.jhu()
population_data = data_loader.population()
pcr_data = data_loader.pcr()
oxcgrt_data = data_loader.oxcgrt()
vaccine_data = data_loader.vaccine()

ita_scenario = cs.Scenario(jhu_data, population_data, "Italy")
ita_scenario.interactive = True
pcr_data.positive_rate("Italy")
ita_scenario.records(variables=["Confirmed", "Infected", "Fatal", "Recovered"], color_dict={"Confirmed": "blue", "Infected": "orange", "Fatal": "red", "Recovered": "green"})
_ = ita_scenario.trend()
ita_scenario.summary()
ita_scenario.estimate(cs.SIRF)
ita_scenario.summary()
# Add future phase to main scenario
ita_scenario.add(name="Main", days=30)
# fit_dict = ita_scenario.fit(oxcgrt_data=oxcgrt_data, name="Forecast")
# ita_scenario.predict(name="Forecast").summary(name="Forecast")
ita_scenario.fit_predict(oxcgrt_data=oxcgrt_data, name="Forecast").summary(name="Forecast")
# Simulation of the number of cases
ita_scenario.summary()
ita_sim_df = ita_scenario.simulate(variables=["Confirmed", "Infected", "Fatal", "Recovered"],
                                   color_dict={"Confirmed": "blue", "Infected": "orange", "Fatal": "red", "Recovered": "green"},
                                   name="Main")
ita_sim_df = ita_scenario.simulate(variables=["Confirmed", "Infected", "Fatal", "Recovered"],
                                   color_dict={"Confirmed": "blue", "Infected": "orange", "Fatal": "red", "Recovered": "green"},
                                   name="Forecast")
# ita_scenario.describe()
_ = ita_scenario.history(target="Confirmed")
_ = ita_scenario.history(target="Infected")
_ = ita_scenario.history(target="Fatal")
_ = ita_scenario.history(target="Recovered")
_ = ita_scenario.history(target="Rt")

And I get this spike result:

I compared this to an old plot I had from early January of Italy estimation:

Which does not have any such spike. Why do you think this is happening now?

Environment

• CovsirPhy version: v2.17.0-eta
• Python version: 3.8.5
• Installation: conda/pip
• System: Windows

[lisphilar]

Because new records were registered, total length of phase series became loger (old: Feb2020 - Jan2021, new: Feb2020 - Mar2021). With this change, the number of phases decreased from 6 to 2 in Feb2020 - Oct2020). Repeated S-R trend analysis for Feb2020 - Oct2020 may solve this issue.

[Inglezos]

Yes I can confirm this, the following is the trends plot with records up to today for Italy:

and this is for records up to 31Dec20:

We can clearly see that in the second case there are many phases in the period March 2020 - August 2020, whereas in the first one there are only two, which is very wrong.

What seems weird to me is why past phases are affected that much by such future phases? I mean these phases are at least 6-8 months apart (April 20 - Feb 21). How can the change points analysis of the records March 20 - May 20 consider data from 2021 and therefore the phases separation is affected so much negatively?

[lisphilar]

I did not take time to investigate the internal algorithms / papers of ruptures, but ruptures has two approaches, top-down and bottom-up. At this time, we use ruptures.Pelt class in covsirphy._SRChange class. ruptures.Pelt is based on top-down approach. This class reads the all records and then determines the change points. Finding criteria may be modified when new data is available with more drastic changes.

On the other hand, bottom-up approach is as follows.

bottom-up segmentation is generous: it starts with many change points and successively deletes the less significant ones. First, the signal is divided in many sub-signals along a regular grid. Then contiguous segments are successively merged according to a measure of how similar they are.

(From ruptures documentation: https://centre-borelli.github.io/ruptures-docs/user-guide/detection/bottomup/)

Note: After version 2.17.0, I updated codes for refactoring and trying to new algorithms. call tree: Scenario.trend() -> ParamTracker.trend() -> trend.trend_detector.TrendDetector.sr() -> trend.sr_change._SRChange.run()

For screening, I tried some ruptures classes in covsirphy._SRchange.run(). (Class parameters are not adjusted.) ruptures.BottomUp(model="normal") showed quite lower RMSLE score, 0.085 with my local PC and "examples/scenario_analysis_py", Italy. Baseline Pelt(model="rbf") showed 0.381.

BottomUp may be helpful, but this change also raises some issues. 64 phases will be produced. We need to evaluate the impact carefully.

Additionally, Binseg(model="normal") showed RMSLE score 0.065, producing 46 phases. This is a top-down approach, but criteria may not be changed because change point analysis will be repeated automatically. https://centre-borelli.github.io/ruptures-docs/user-guide/detection/binseg/

All results:

setting	setting	trend()	.estimate()	.score()
class	model	phases	Time [sec]	RMSLE
Pelt	l1	0th-1st	NA	NA
Pelt	l2	0th	NA	NA
Pelt	rdf	0th-13th	294	0.381
Binseg	l1	0th-1st	NA	NA
Binseg	l2	0th	NA	NA
Binseg	rbf	0th-13th	202	0.334
Binseg	linear	Error	NA	NA
Binseg	normal	0th-45th	284	0.065
Binseg	ar	0th	NA	NA
Window	l1	0th-1st	NA	NA
Window	l2	0th	NA	NA
Window	rbf	0th-1st	NA	NA
Window	linear	Error	NA	NA
Window	normal	0th-2nd	NA	NA
Window	ar	0th	NA	NA
BottomUp	l1	0th-1st	NA	NA
BottomUp	l2	0th	NA	NA
BottomUp	rbf	0th-15th	213	0.339
BottomUp	linear	Error	NA	NA
BottomUp	normal	0th-63th	365	0.085
BottomUp	ar	0th	NA	NA

[lisphilar]

With #678, we can select algorithms of change detection.

snl = cs.Scenario(country="Italy")
snl.register(jhu_data, population_data)
snl.trend(algo="Binseg-normal")

Argument algo will be selected from "Pelt-rbf", "Binseg-rbf", "Binseg-normal", "BottomUp-rbf" and "BottomUp-normal". At this time, the default value is "Pelt-rbf" and this means we use the same algoritm as that in the previous version.

snl.trend(algo="Binseg-normal") shows the next figure and title and the other settings should be changed later.

[Inglezos]

I tried some combinations myself and the Binseg-normal seems to be indeed the better method for the phases separation (trade-off between accuracy of phases separation, accuracy of estimation and number of phases).

I have also another idea for this to work, perhaps we should change the min_size as well. In my opinion, it would be more logical that each phase has a minimum duration of a couple of weeks, in order the effects of measures or infectivity to be represented in the recorded cases (as we assume after all in the forecast mechanism), but also for the ODE system to have enough data to work on for each phase. So in my tests I tried min_size = 14 along with Binseg-normal and the results were better while keeping acceptable number of phases, around 20. Please try and confirm this whenever you have some time.

In order to reach to a final solution, we have to agree on a combination of these three factors (i.e. Binseg, normal, 14). To take this decision we have to check the behavior for multiple countries, both in trends analysis and in estimation and accept the trade-off. Of course an alternative would be to keep the current configuration and apply repeatedly trends analysis on very long sub-phases, but this could be an improvement for the future instead.

[lisphilar]

Thank you for your confirmation. "Binseg-normal" could be the default value for algo. I agree with the idea we should adjust min_size, but is is a difficult task to specify the default value of min_size for all countries. How about checking total MSE value on S-R plane (without parameter estimation) with some min_size values to determine appropriate value of min_size for a country?

With the next codes and results, min_size=15 (sharp change from 15 to 16) or min_size=7 seems appropriate value for Italy. (I calculated sum of RMSLE scores in phases as a tentative method here, but RMSLE is not suitable for production because RMSLE has root.)

import covsirphy as cs
data_loader = cs.DataLoader()
jhu_data = data_loader.jhu()
population_data = data_loader.population()

country = "Italy"
record_df, _ = jhu_data.records(country=country, population=population_data.value(country))
for i in range(3, 31):
    detector = cs.TrendDetector(record_df, min_size=i)
    _ = detector.sr(algo="Binseg-normal")
    df = detector.summary()
    rmsle_scores = df["RMSLE_S-R"]
    print(f"RMSLE with min_size={i:2}: {rmsle_scores.sum():.6f}")

RMSLE with min_size= 3: 0.001635 RMSLE with min_size= 4: 0.001591 RMSLE with min_size= 5: 0.001564 RMSLE with min_size= 6: 0.001560 RMSLE with min_size= 7: 0.001507 RMSLE with min_size= 8: 0.001625 RMSLE with min_size= 9: 0.001687 RMSLE with min_size=10: 0.001662 RMSLE with min_size=11: 0.001662 RMSLE with min_size=12: 0.001662 RMSLE with min_size=13: 0.001662 RMSLE with min_size=14: 0.001658 RMSLE with min_size=15: 0.001752 RMSLE with min_size=16: 0.002479 RMSLE with min_size=17: 0.002476 RMSLE with min_size=18: 0.002587 RMSLE with min_size=19: 0.002587 RMSLE with min_size=20: 0.002465 RMSLE with min_size=21: 0.002482 RMSLE with min_size=22: 0.002483 RMSLE with min_size=23: 0.002491 RMSLE with min_size=24: 0.002549 RMSLE with min_size=25: 0.002599 RMSLE with min_size=26: 0.002599 RMSLE with min_size=27: 0.002599 RMSLE with min_size=28: 0.002599 RMSLE with min_size=29: 0.002600 RMSLE with min_size=30: 0.002600

Or, if we need to keep the total number of phases, "dynamic programming" algorithm may be better. https://centre-borelli.github.io/ruptures-docs/user-guide/detection/dynp/

[Inglezos]

Or, if we need to keep the total number of phases, "dynamic programming" algorithm may be better.

We don't know beforehand the number of phases. We just need to keep an eye out for this number, because it seems wrong for them to be, let's say, greater than 30 (for now, in the future it would be normal as months are passing by), while the actual could be just 14. It wouldn't be right to impose this restriction to our analysis beforehand. Dynamic programming has as a precondition that this number is known right? Could it be combined instead (for unknown number) with a penalized method?

I agree strongly with the idea that min_size should vary for each country by running a parametric analysis iteratively to select the best value, but this should also consider the number of phases as mentioned above. A big foul for me is that for example for Italy, each phase has a duration of almost the min_size value, around 7-10 days, and it is like the algorithm uses this value blindly without freedom in the phases selection. Like the algorithm is biased to this input value of min_size. While if we use larger min_size value, the algorithm is forced in a way to actually work on the phases separation and use different intervals and not necessarily similar to min_size.

[lisphilar]

Just to confirm, in your idea, the default value is Scenario.trend(algo="Binseg-normal", limits=(3, 31))? limit is the minimum and maximum value of min_size. (If you have any ideas for this new argument, please let me know.)

This method automatically calculate total MSE scores for min_size=i in the range of i=3, 4, 5,....31 within the limits(=(3, 31) as default) and select the best min_size, then determine the change points.

Is this correct?

[Inglezos]

Yes exactly that's what I have in mind, this is correct.

The new factor that we have to somehow include as well is how much the results are unbiased from the min_size and the phases duration are not all very close to that value (since it seems improbable to be actually true indeed!).

We have to achieve the sweet spot between number of phases (and therefore unbiased phases separation in order the ruptures algorithm to actually work on the data) and separation MSE, keeping in mind that the subsequent estimation should lead to good results with acceptable deviation.

Therefore, the upper limit for now, intuitively can be around 25 or even 30, but in the future this will certainly increase inevitably as weeks are passing by, so we need to calculate this upper limit more strictly somehow, associating it with the aforementioned unbiased phases separation dilemma I presented.

For the lower limit, I feel it should be at least one week, because I think that the ODE system uses these days data to reach to a solution. And if we have too few data, which means only few days, then the ODE solution will not be stable/robust, I have this concern. In the ODE system for the SIRF model we have 5 unknown variables/parameters, the alpha1, alpha2, beta, gamma and tau, right? So for each phase we need at least 5 sets of cases data points for the system to converge to a solution, correct? I am not sure if this means we need data from at least 5 days, since if I remember correctly tau does not correspond 1-1 to days, but you get my point I hope.

[lisphilar]

Thank you for your confirmation and the new idea regarding phase duration.

For lower limit, I agree with the idea phase duration should be over the number of variables or parameters (including tau).

For upper limit, just for discussion, with #680, I updated TrendDetector.summary(), a development tool for trend analysis. (This is not linked with Scenario.trend() and Scnerio.summary().)

In Italy case, phase duration determined by ruptures seems not free from min_size, but the relationship is non-linear, if my interpretation is correct.

import covsirphy as cs
data_loader = cs.DataLoader()
jhu_data = data_loader.jhu()
population_data = data_loader.population()
country = "Italy"
record_df, _ = jhu_data.records(country=country, population=population_data.value(country))

for i in range(3, 80):
    detector = cs.TrendDetector(record_df, min_size=i)
    _ = detector.sr(algo="Binseg-normal")
    df = detector.summary(metrics="MSE")
    score = df["MSE_S-R"].sum()
    duration_diff = df["Duration"].mean() - i
    print(f"min_size={i:2}: MSE={score:.2g}, [mean of (duration - min_size)]={int(duration_diff):2}")

min_size= 3: MSE=4.1e+08, [mean of (duration - min_size)]= 0 min_size= 4: MSE=4.8e+08, [mean of (duration - min_size)]= 2 min_size= 5: MSE=5.1e+08, [mean of (duration - min_size)]= 2 min_size= 6: MSE=5.1e+08, [mean of (duration - min_size)]= 1 min_size= 7: MSE=4.7e+08, [mean of (duration - min_size)]= 2 min_size= 8: MSE=9.4e+08, [mean of (duration - min_size)]= 2 min_size= 9: MSE=1e+09, [mean of (duration - min_size)]= 3 min_size=10: MSE=1.1e+09, [mean of (duration - min_size)]= 5 min_size=11: MSE=1.2e+09, [mean of (duration - min_size)]= 5 min_size=12: MSE=1.2e+09, [mean of (duration - min_size)]= 4 min_size=13: MSE=1.2e+09, [mean of (duration - min_size)]= 3 min_size=14: MSE=1.2e+09, [mean of (duration - min_size)]= 3 min_size=15: MSE=1.3e+09, [mean of (duration - min_size)]= 4 min_size=16: MSE=6.8e+09, [mean of (duration - min_size)]= 4 min_size=17: MSE=6.8e+09, [mean of (duration - min_size)]= 5 min_size=18: MSE=7e+09, [mean of (duration - min_size)]= 5 min_size=19: MSE=7e+09, [mean of (duration - min_size)]= 8 min_size=20: MSE=7e+09, [mean of (duration - min_size)]= 9 min_size=21: MSE=7e+09, [mean of (duration - min_size)]=10 min_size=22: MSE=7e+09, [mean of (duration - min_size)]= 9 min_size=23: MSE=7e+09, [mean of (duration - min_size)]= 8 min_size=24: MSE=7.1e+09, [mean of (duration - min_size)]=10 min_size=25: MSE=7.1e+09, [mean of (duration - min_size)]=13 min_size=26: MSE=7.1e+09, [mean of (duration - min_size)]=12 min_size=27: MSE=7.1e+09, [mean of (duration - min_size)]=11 min_size=28: MSE=7.1e+09, [mean of (duration - min_size)]=10 min_size=29: MSE=7.1e+09, [mean of (duration - min_size)]= 9 min_size=30: MSE=7.1e+09, [mean of (duration - min_size)]= 8 min_size=31: MSE=7.1e+09, [mean of (duration - min_size)]=11 min_size=32: MSE=7.1e+09, [mean of (duration - min_size)]=10 min_size=33: MSE=8.1e+09, [mean of (duration - min_size)]= 9 min_size=34: MSE=8.9e+09, [mean of (duration - min_size)]= 8 min_size=35: MSE=1e+10, [mean of (duration - min_size)]= 7 min_size=36: MSE=2.5e+10, [mean of (duration - min_size)]=12 min_size=37: MSE=2.5e+10, [mean of (duration - min_size)]=11 min_size=38: MSE=2.9e+10, [mean of (duration - min_size)]=17 min_size=39: MSE=2.9e+10, [mean of (duration - min_size)]=16 min_size=40: MSE=2.9e+10, [mean of (duration - min_size)]=15 min_size=41: MSE=2.9e+10, [mean of (duration - min_size)]=14 min_size=42: MSE=2.9e+10, [mean of (duration - min_size)]=13 min_size=43: MSE=2.9e+10, [mean of (duration - min_size)]=12 min_size=44: MSE=2.9e+10, [mean of (duration - min_size)]=11 min_size=45: MSE=2.9e+10, [mean of (duration - min_size)]=10 min_size=46: MSE=2.9e+10, [mean of (duration - min_size)]= 9 min_size=47: MSE=2.9e+10, [mean of (duration - min_size)]=17 min_size=48: MSE=2.9e+10, [mean of (duration - min_size)]=16 min_size=49: MSE=2.9e+10, [mean of (duration - min_size)]=15 min_size=50: MSE=2.9e+10, [mean of (duration - min_size)]=14 min_size=51: MSE=2.9e+10, [mean of (duration - min_size)]=13 min_size=52: MSE=2.9e+10, [mean of (duration - min_size)]=12 min_size=53: MSE=2.9e+10, [mean of (duration - min_size)]=11 min_size=54: MSE=2.9e+10, [mean of (duration - min_size)]=10 min_size=55: MSE=2.9e+10, [mean of (duration - min_size)]= 9 min_size=56: MSE=2.9e+10, [mean of (duration - min_size)]= 8 min_size=57: MSE=2.9e+10, [mean of (duration - min_size)]= 7 min_size=58: MSE=2.9e+10, [mean of (duration - min_size)]= 6 min_size=59: MSE=2.9e+10, [mean of (duration - min_size)]= 5 min_size=60: MSE=2.9e+10, [mean of (duration - min_size)]=17 min_size=61: MSE=2.9e+10, [mean of (duration - min_size)]=16 min_size=62: MSE=2.9e+10, [mean of (duration - min_size)]=15 min_size=63: MSE=2.9e+10, [mean of (duration - min_size)]=14 min_size=64: MSE=2.9e+10, [mean of (duration - min_size)]=13 min_size=65: MSE=2.9e+10, [mean of (duration - min_size)]=12 min_size=66: MSE=2.9e+10, [mean of (duration - min_size)]=11 min_size=67: MSE=2.9e+10, [mean of (duration - min_size)]=10 min_size=68: MSE=2.9e+10, [mean of (duration - min_size)]= 9 min_size=69: MSE=2.9e+10, [mean of (duration - min_size)]= 8 min_size=70: MSE=2.9e+10, [mean of (duration - min_size)]= 7 min_size=71: MSE=2.9e+10, [mean of (duration - min_size)]= 6 min_size=72: MSE=2.9e+10, [mean of (duration - min_size)]= 5 min_size=73: MSE=2.9e+10, [mean of (duration - min_size)]= 4 min_size=74: MSE=1.8e+10, [mean of (duration - min_size)]=23 min_size=75: MSE=1.8e+10, [mean of (duration - min_size)]=22 min_size=76: MSE=1.8e+10, [mean of (duration - min_size)]=21 min_size=77: MSE=1.8e+10, [mean of (duration - min_size)]=20 min_size=78: MSE=1.8e+10, [mean of (duration - min_size)]=19 min_size=79: MSE=1.8e+10, [mean of (duration - min_size)]=18

[lisphilar]

Revised my comment: Phase duration is NOT free from min_size, but the relationship is non-linear.

[lisphilar]

New idea:

1. We set Scenario.trend(algo="Binseg-normal", min_size=None) as default.
2. When min_size is None, output of Scenario.estimate_delay() will be used.
3. Detect changes with Binseg(model="normal", min_size=3): this min_size is not the same as Scenario.trend(min_size). Always 3.
4. Phases with < min_size (argument of Scenario.trend()) days will be merged with the previous phases.

[Inglezos]

When I try to execute Scenario.estimate_delay() (for any country) it throws the error:

KeyError: "@indicator and target must be a sub-list of Date, Confirmed, Infected, Fatal, Recovered, Susceptible, but ['Stringency_index', 'Confirmed'] was applied."

Is this normal or create a new issue? I run:

import covsirphy as cs
data_loader = cs.DataLoader()
jhu_data = data_loader.jhu()
population_data = data_loader.population()
ita_scenario = cs.Scenario(jhu_data, population_data, "Italy")
ita_scenario.estimate_delay()

Am I doing something wrong?

[Inglezos]

I like the new idea yes. This fixed 3 value where is exactly used in the code? Due to the problem mentioned in my previous comment, I cannot see the estimated delay for some example countries in order to check whether this value is good enough to be a min_size value.

[lisphilar]

The KeyError was raised because OxCGRTData was not registered in Scenario instance. Please run ita_scenario.register(extras=[oxcgrt_data]) in advance.

https://lisphilar.github.io/covid19-sir/usage_quick.html#Start-scenario-analysis

At the current version, the minimum value of min_size is 3 because we cannot perform curve fitting (to show the .trend() figure) with 1 or 2 points.

[Inglezos]

Yes you're right, I thought that this dataset was pre-registered by default. Now back to the idea.. I think it is a very good idea, but for example for Italy that period is 21 days. If we use in the trends min_size = 21 though, we have a not so good phases separation. Also that applies for France as some first examples I tried. For USA, Greece it was a very good approach though so I guess it depends. For Italy and France this leads to: One idea would be to subtract few days (1 week) from the delay period or try another indicator-target pair?

[Inglezos]

Another idea would be to use 7 days as minimum phase duration and apply the following logic:

snl_scenario = cs.Scenario(jhu_data, population_data, "Italy")
snl_scenario.register(extras=[oxcgrt_data])
snl_delay = snl_scenario.estimate_delay(indicator="Stringency_index", target="Confirmed")

if snl_delay[1]["Period Length"].mode()[0] >= 7:
    phase_period = snl_delay[1]["Period Length"].mode()[0]
elif snl_delay[0] >=7:
    phase_period = snl_delay[0]
else:
    phase_period = 7

_ = snl_scenario.trend(algo="Binseg-normal", min_size=phase_period)

or something like this.

[lisphilar]

Or, can we set Scenario.trend(algo="Binseg-normal", min_size=7) as default?

[lisphilar]

Thank you. Your idea is to use mode value, not mean/median? Should we also change the default output of .estimate_delay()[0] from mean to median?

Italy: mean: 22 median: 21 mode: 14

The folowing figures (histogram/kde of delay period length) suggest that "mode" value will be a better choice.

from matplotlib import pyplot as plt
import covsirphy as cs
data_loader = cs.DataLoader()
jhu_data = data_loader.jhu()
population_data = data_loader.population()
oxcgrt_data = data_loader.oxcgrt()
country = "Italy"
snl_scenario = cs.Scenario(jhu_data, population_data, country)
snl_scenario.register(extras=[oxcgrt_data])
try:
    _, delay_df = snl_scenario.estimate_delay(indicator="Stringency_index", target="Confirmed")
    series = delay_df["Period Length"]
    series.plot.hist()
    plt.show()
    series.plot.kde()
    plt.show()
    print(f"mean: {int(series.mean())}")
    print(f"median: {int(series.median())}")
    print(f"mode: {int(series.mode()[0])}")
except KeyError:
    pass

[Inglezos]

Yes I feel that the "mode" value would be more representative to use.

[lisphilar]

With #685, I update codes.

• TrendDetector(..., min_size=7) as default
• TrendDetector(algo="Binseg-normal") as default
• When min_size=None with Scenario.trend(), min_size for S-R trend anlysis will be set as the delay period calculated with Scenario.estimate_delay()

• 686: use mode value of delay period

[lisphilar]

With #687, the title will be updated to simplify the figure. I checked the performance of new Scenario.trend() with the following script.

import covsirphy as cs
data_loader = cs.DataLoader()
jhu_data = data_loader.jhu()
population_data = data_loader.population()
oxcgrt_data = data_loader.oxcgrt()
countries = [
    "Italy", "India", "USA", "Greece", "Russia",
    "Brazil", "France", "Spain", "United Kingdom", "New Zealand", "Germany", "Japan",
]
filer = cs.Filer(directory=f"example/output/trend", prefix="trend", numbering="01")
for country in countries:
    snl = cs.Scenario(country=country, province=None)
    snl.register(jhu_data, population_data, extras=[oxcgrt_data])
    snl.interactive = False
    snl.trend(**filer.png(country.replace(" ", "_")))

[Inglezos]

Thank you, I think the results are very good! My test countries are the following by the way:

countries_regression = [
    "Italy", "Japan", "China", "India", "USA",
    "Greece", "Brazil", "Russia", "Spain", "France",
    "United Kingdom", "Germany", "Netherlands", "Czech Republic", "Poland",
    "Sweden", "Austria", "Australia", "Switzerland", "Hungary", "South Africa"
]

[lisphilar]

Thank you for your ideas and testing with many countries! I will close this issue.

[Inglezos]

As I now realize, I had run the trends tests with an older version and not having git pulled the updated codes, I am very sorry. The results I get now for many of the countries I mentioned above are not good (if you have some time check this yourself too). Probably we should revisit this issue.

[lisphilar]

I retried testing with version 2.17.0-omicron.

import covsirphy as cs
data_loader = cs.DataLoader()
jhu_data = data_loader.jhu()
population_data = data_loader.population()
oxcgrt_data = data_loader.oxcgrt()
countries = [
    "Italy", "Japan", "China", "India", "USA",
    "Greece", "Brazil", "Russia", "Spain", "France",
    "United Kingdom", "Germany", "Netherlands", "Czech Republic", "Poland",
    "Sweden", "Austria", "Australia", "Switzerland", "Hungary", "South Africa"
]
filer = cs.Filer(directory=f"example/output/trend", prefix="trend", numbering="01")
for country in countries:
    snl = cs.Scenario(country=country, province=None)
    snl.register(jhu_data, population_data, extras=[oxcgrt_data])
    snl.interactive = False
    snl.trend(**filer.png(country.replace(" ", "_")))
    print(country, snl.estimate_delay()[0])

I could not find un-expected results. What countries have un-expected results?

China: The 2nd phase is not fitted. However, to fix this issue, the phase duration will be < 5.

Spain: Accuracy of the 14th phase is lower due to un-updated Recovered data. Complement system could be updated?

Netherlands: The 10th phase could be divided, but I could not find the cause...

Hungary: Abnormal detection due to delay period = 53.

[lisphilar]

We may need to revise Scenario.estimate_delay() to get many data for "Period Length" calculation. We have only 11 records for Hungary as an example.

snl = cs.Scenario(country="Hungary", province=None)
snl.register(jhu_data, population_data, extras=[oxcgrt_data])
# 53, dataframe
delay, df = snl.estimate_delay()

Index	Confirmed	Stringency_index	Period Length
0	48757	40.74	NaN
1	32	50	12
2	32	50	36
3	32	50	29
4	4183	54.63	17
5	94916	57.41	53
6	3990	61.11	31
7	3741	66.67	53
8	300	71.3	2
9	366279	72.22	5
10	2443	76.85	118
11	614612	79.63	19

[lisphilar]

This issue #670 is closed to release new stable version 2.18.0. The issue for Hungary will be discussed in #693 and fixed in version 2.19.0 or 2.18.1.

Note: When OxCGRTData is not registered, delay period will not be used and min_size will be 7.

snl = cs.Scenario(country="Hungary", province=None)
snl.register(jhu_data, population_data)
snl.trend()

[Inglezos]

My findings are:

1. Spain phases separation faces issues -> 13th, 14th
2. Netherlands phases separation faces issues -> 9th, 10th
3. Czech Republic phases separation faces issues -> 2nd, 12th
4. Australia phases separation faces issues -> 16th
5. Hungary phases separation is totally problematic
6. South Africa phases separation faces issues -> 14th

Also the major one is: 0a. Japan, China, India, USA, Greece, Rrussia, UK, Germany, Poland, Sweden, Austria, Australia, Switzerland have too many phases, I would expect normally around 25 max. This is related to the dilemma I described, about phases duration dependence on small min_size. I think that in these countries, the min_size has the minimum value 7 and there are many short phases. 0b. Is this necessarily wrong though? Can the ODE system converge and give a stable solution with only 7 timepoints having 5 unknown variables), is it enough? Because if the 7 points are enough and lead to good results, both in phases separation and estimation terms, then we don't need to worry about the high number of phases.

[lisphilar]

We should revise .estimate_delay() to use appropriate min_size values because .estimate_delay() also uses ruptures package and Pelt class is used there. This will be discussed in #693.

I think 7 timepoints are enough for ODE system, but this could be tested with ModelValidator class. (In ModelValidator class, the number of timepoints = step_n is selected randomly at this time. We need to update codes of analysis/model_validator.py line 102 etc.)

Additionally, many phases sometimes lead too long estimation time. This is not the first priority, but this should be also considered.

[Inglezos]

This discussion will continue in #693 with .estimate_delay() reworks.