Repository: cmu-delphi/covidcast-indicators
Issue: changehc indicator outputs no data for hhs, nation when data is present at lower geo levels

[[Slack] context](https://delphi-org.slack.com/archives/C01LZ3A2UMU/p1614698588017700?thread_ts=1614663716.017100&cid=C01LZ3A2UMU)

In the changehc indicator, it seems to be possible for county/state/hrr data to pass minimum denominator thresholds when hhs/nation data does not. Our current hypothesis is that this is due to smoothing numerators and denominators separately and after geo aggregation: the smoother is being overly enthusiastic about curve fitting, such that some counties/state/hrrs may be concave-up when the majority of their siblings are concave-down, causing the overall curve for hhs or nation to fall below threshold.

There are several suspicious things going on here:
* Smoothing after geo aggregation (though doctor-visits also does this)
* Smoothing the numerator and denominator separately before computing the rate
* Backfill causes the counts to be artificially low for the most recent period of data, and the smoother is probably not a good tool for modeling this effect
* ~Zero-filling (iirc the smoother expects NaNs, not zeros)~ (this turns out to be true zeros -- by definition CHNG does not report if the count is 0)
* No corner-case tests, only completion tests
* ~[how did this ever work](https://github.com/cmu-delphi/covidcast-indicators/blob/ca43cda9ffd3e321369ad7911c5d8f44a1feb70a/changehc/delphi_changehc/sensor.py#L152)~ (turns out pandas assigns NaNs for the remaining indices, which is fine)

To fix:
* Write a test that fails when the relationship between having-data is not hierarchical
* Consider smoothing ratios instead of numerators/denominators
* Switch to using raw denominators for thresholding
* Consider whether backfill is adding 2x smoothing
* Consider reducing the bandwidth? though this would have other ramifications

[Sample data](http://rinkitink.ml.cmu.edu/delphi/misc/bad_spiky_business.csv)

Changed files:
changehc/tests/test_update_sensor.py, changehc/delphi_changehc/sensor.py

diff --git a/changehc/delphi_changehc/sensor.py b/changehc/delphi_changehc/sensor.py
index 32c23f66..f4a8934a 100644
--- a/changehc/delphi_changehc/sensor.py
+++ b/changehc/delphi_changehc/sensor.py
@@ -26,19 +26,6 @@ class CHCSensor:
                         poly_fit_degree=1,
                         gaussian_bandwidth=Config.SMOOTHER_BANDWIDTH)
 
-    @staticmethod
-    def gauss_smooth(count,total):
-        """Smooth using the left_gauss_linear.
-
-        Args:
-            count, total: array
-        """
-        count_smooth = CHCSensor.smoother.smooth(count)
-        total_smooth = CHCSensor.smoother.smooth(total)
-        total_clip = np.clip(total_smooth, 0, None)
-        count_clip = np.clip(count_smooth, 0, total_clip)
-        return count_clip, total_clip
-
     @staticmethod
     def backfill(
             num,
@@ -120,30 +107,13 @@ class CHCSensor:
 
         # calculate smoothed counts and jeffreys rate
         # the left_gauss_linear smoother is not guaranteed to return values greater than 0
-
-        smoothed_total_counts, smoothed_total_visits = CHCSensor.gauss_smooth(
-            total_counts.flatten(), total_visits
-        )
-
-        # in smoothing, the numerator may have become more than the denominator
-        # simple fix is to clip the max values elementwise to the denominator (note that
-        # this has only been observed in synthetic data)
-        # smoothed_total_counts = np.clip(smoothed_total_counts, 0, smoothed_total_visits)
-
-        smoothed_total_rates = (
-                (smoothed_total_counts + 0.5) / (smoothed_total_visits + 1)
-        )
-
-        # checks - due to the smoother, the first value will be NA
-        assert (
-                np.sum(np.isnan(smoothed_total_rates[1:])) == 0
-        ), "NAs in rate calculation"
-        assert (
-                np.sum(smoothed_total_rates[1:] <= 0) == 0
-        ), f"0 or negative value, {geo_id}"
+        rates = total_counts.flatten() / total_visits
+        smoothed_rate = CHCSensor.smoother.smooth(rates)
+        clipped_smoothed_rate = np.clip(smoothed_rate, 0, 1)
+        jeffreys_rate = (clipped_smoothed_rate * total_visits + 0.5) / (total_visits + 1)
 
         # cut off at sensor indexes
-        rate_data = pd.DataFrame({'rate':smoothed_total_rates, 'den': smoothed_total_visits},
+        rate_data = pd.DataFrame({'rate': jeffreys_rate, 'den': total_visits},
                                  index=y_data.index)
         rate_data = rate_data[first_sensor_date:]
         include = rate_data['den'] >= Config.MIN_DEN
diff --git a/changehc/tests/test_update_sensor.py b/changehc/tests/test_update_sensor.py
index dd8ca5dd..3ec2dccf 100644
--- a/changehc/tests/test_update_sensor.py
+++ b/changehc/tests/test_update_sensor.py
@@ -81,18 +81,24 @@ class TestCHCSensorUpdator:
                 ""
             )
             su_inst.shift_dates()
-            data_frame = su_inst.geo_reindex(self.small_test_data.reset_index())
+            test_data = pd.DataFrame({
+                "num": [0, 100, 200, 300, 400, 500, 600, 100, 200, 300, 400, 500, 600],
+                "fips": ['01001'] * 7 + ['04007'] * 6,
+                "den": [1000] * 7 + [2000] * 6,
+                "date": [pd.Timestamp(f'03-{i}-2020') for i in range(1, 14)]})
+            data_frame = su_inst.geo_reindex(test_data)
             assert data_frame.shape[0] == multiple*len(su_inst.fit_dates)
             assert (data_frame.sum() == (4200,19000)).all()
 
     def test_update_sensor(self):
         """Tests that the sensors are properly updated."""
-        for geo in ["state","hrr"]:
+        outputs = {}
+        for geo in ["county", "state", "hhs", "nation"]:
             td = TemporaryDirectory()
             su_inst = CHCSensorUpdator(
-                "02-01-2020",
-                "06-01-2020",
-                "06-12-2020",
+                "03-01-2020",
+                "03-22-2020",
+                "03-27-2020",
                 geo,
                 self.parallel,
                 self.weekday,
@@ -100,10 +106,27 @@ class TestCHCSensorUpdator:
                 self.se,
                 ""
             )
-            su_inst.update_sensor(self.small_test_data,  td.name)
+            # As of 3/3/21 (40c258a), this set of data has county outputting data, state and hhs not
+            # outputting data, and nation outputting data, which is undesirable. Ideal behaviour
+            # should be all output or a subregion only outputting if its parent has output,
+            # which is what is being tested here.
+            small_test_data = pd.DataFrame({
+                "num": [0, 100, 200, 300, 400, 500, 600, 100, 200, 300, 400, 500, 600] * 2,
+                "fips": ["01001"] * 13 + ["42003"] * 13,
+                "den": [30, 50, 50, 10, 1, 5, 5, 50, 50, 50, 0, 0, 0] * 2,
+                "date": list(pd.date_range("20200301", "20200313")) * 2}).set_index(
+                ["fips", "date"])
+            su_inst.update_sensor(small_test_data,  td.name)
+            for f in os.listdir(td.name):
+                outputs[f] = pd.read_csv(os.path.join(td.name, f))
             assert len(os.listdir(td.name)) == len(su_inst.sensor_dates),\\
                 f"failed {geo} update sensor test"
             td.cleanup()
+        assert outputs["20200319_county_smoothed_outpatient_covid.csv"].empty
+        assert outputs["20200319_state_smoothed_outpatient_covid.csv"].empty
+        assert outputs["20200319_hhs_smoothed_outpatient_covid.csv"].empty
+        assert outputs["20200319_nation_smoothed_outpatient_covid.csv"].empty
+
 
 class TestWriteToCsv:
     """Tests for writing output files to CSV."""