[fb-survey] synthetic dataset and unit tests

[krivard]

• spec out what kind of data is needed
• generate data or adapt existing data to be delphi-safe
• complete tests

[capnrefsmmat]

See https://github.com/cmu-delphi/covid-19/pull/73

I'll leave this open for @statsmaths to note any additions or changes that have to be done to the synthetic data, and to track the cleanup and testing of the Facebook code.

[krivard]

Collecting here a list of cases the code habitually needs to handle that won't necessarily get caught by the synthetic data:

• Should accept multiple --full arguments
• Should only keep the earliest record for a token when sorted by StartDate
• Should merge --full arguments with a savefile when --load is specified, keeping the earliest record for a token when sorted by StartDate, and excluding that record from recoded.cleaned.recorded.partA iff it was originally excluded (debugging this was no fun at all)

Integrity checks can't be run with any kind of synthetic data. They include:

• Make sure the local files in git are up-to-date with what's currently in the API
• Make sure the updates to local files (we get ~4 days of backfill with a small number of respondents) are small
• Make sure there are no big jumps in API columns vs what's currently in the API
• Make sure there are no big jumps in effective_sample_size vs what's currently in git

We used to need to handle files that were encoded in Mountain time instead of Pacific, but with the Qualtrics API that's not necessary anymore. It should stay in, just in case anyone ever needs to load an old file (generated before 20 April), but it's not a testing priority.

[taylor-arnold]

I have had a chance now to starting looking at the synthetic data and existing code. There is a lot of code in the current pipeline that has been copied and slightly changed (with little to no comments). This makes it nearly impossible to try to reorganise the way I did for JHU and GHT, namely refactoring the code into smaller functions and cleaning up as I went along. It seems like a better approach will be starting over with a new code base, and borrowing existing code chunks as possible. In order to do this, I need to understand what the code is supposed to be doing. Note that I am (at least currently) only attacking what's covered in the R code inside of the covid-19/facebook directory and not consider the pull of data from the directory covid-19/qualtrics-cron.

[taylor-arnold]

Here is my current understanding of the code and data:

1. The input dataset processed by the R code has the structure of the synthetic data constructed here. This is a csv file with additional information pushed into the header. Each record has a field token that we use to match with the weights from Facebook and the field A3 captures the self-reported zip code. Currently we are interested in the responses A1_1 through A1_5, which are of the form "Does someone in your household have this symptom" and the data are recorded as 1="No", 2="Yes", and NA="Did not answer".
2. Our main variables of interest are CLI (does someone have a fever and one of: cough, shortness of breath, difficulty breathing) and ILI (does someone have a fever and one of: cough, soar throat).
3. There are two additional files that give "weights" for each of the households.
4. The output variables are given at four geographic levels (hrr, msa, county, state), given as a raw and smoothed version, and given as both weighted and unweighted forms. That yields a total of 4x2x2x2 = 32 total indicators.

Here are the outstanding questions I have and was not able to find in the code after looking at it for the better part of the day:

1. Is there something incorrect in my high-level outline above?
2. What method do we use to go from a pair of vectors (responses of CLI or ILI, weights) for a single geographic region to an estimated value, standard error, and sample size? Is the un-weighted version the same as the weighted version with uniform weights, or is there
3. There are several date fields in the dataset. Which one is used to assign a survey to a specific date in time?
4. What is the difference between the two sets of weight files? When do we use one and when do we use the other?
5. How is the smoothed signal computed? Does it literally smooth the raw signal, or does it require going back to the raw data?
6. I see the metadata files in "geographical_scope" and lots of code, but am lost in the weeds. What's the general strategy for aggregating zip code up to the other geographic levels? Summing raw counts? Population weights? A particular pre-computed transition matrix?

Apologies for the all the questions. I've done my best going through the code, but I am still having trouble with some of the details. If there is some other documentation no in the main GH that would answer some of these, please point we to them as well!

[taylor-arnold]

Working on creating a single script that replicates a cleaned-up version of main logic of the existing code. Currently committed on the fb branch with: 4fc709b. See the file here for the current code: test_script.R

[krivard]

Monumentous!

Some additional weird logic: Currently I call run.R three times (once for each of three different output sets), and while we always start at 00:00:00 and end at 23:59:59, each output needs a different date range (it's gross, i know). Taking the run from 2020-05-21 as an example:

--output	--start	--end	--load?
cids	2020-05-20	2020-05-20	no
individual	2020-05-14	2020-05-20	no
covidalert; archive; community	2020-04-06	2020-05-20	yes

Why this is necessary:

1. The CID lists have to bracket a single day of data
2. The individual response sets need to receive small amounts of backfill for 4 days but then should level off; we use the extra days to double check there are no more changes
3. The smoothed aggregations have to run from the beginning of the study or else the figures diverge from what's in the API

We can probably merge (1) and (2) together by just saying only output CIDs for the most recent day in the range. Then if we solve (3), we can put everything on the same 7-day range. If not, we can probably do a similar thing with "only output individual responses for the last seven days".

We have occasionally needed to do mass re-generations of cumulative CID lists for debugging purposes, so whatever we do should permit control through params.json even if there is a standard configuration for the normal automated behavior.

---

Small nit: you probably meant name instead of "full" here:

write_cid <- function(data, name)
{
  fname <- sprintf(
    "cvid_cids_%s_response_%s_-_%s.csv",
    "full",
    format(start_time, "%H_%M_%Y_%m_%d", tz = tz_to),
    format(end_time, "%H_%M_%Y_%m_%d", tz = tz_to)
  )

[taylor-arnold]

Thanks for the additional information and catching the bug in write_cid. I was confused about where the past data was being merged in for the smoothing step; good to know that I was not missing anything but it simply has not been implemented. I think we can build in some logic to avoid having to run the entire thing three times, while still making it flexible for backfill. A brute-force solution would be to include three different start dates, but as you suggest perhaps something more elegant can be done as well.

Likely I am misunderstanding something in the smoothing logic, but I would have thought that was only a need to go back 7 days given that it seems that smoothing only occurs within a 7-day window.

I pushed a revised version of my single script that now handles the smoothing of the signal, grabbing the community data, and storing the archive. From there, it should be close to ready to wrap up into documented and tested functions. It will need some very close code review, as it is not as easy for me to directly verify that it produces output consistent with the existing code.

[taylor-arnold]

I pushed a commit today that has a packaged version of the refacted/rewritten facebook survey code. You can see the current structure here.

Generally, it's well covered with both unit and integration tests. For integration tests, I have no real data to compare against, so I made another very small set with only 20 responses put into just 2 zip codes to try to sanity check that the results are reasonable. While its possible that there are some edge-case bugs, I believe that the current code is mostly accurate for the archive, weights, and individual outputs. For the community and covidcast outputs, I believe that the raw and weighted versions are also correct.

There are now two blocking things that I am not entirely understanding. Namely:

1. I understand the standard binary point estimate and standard errors produced for the community survey. However, I am not familiar with the methods being used for the household counts (or at least, do not recognise the method from the code). I believe my code's logic follows what's in the current script, but hard to feel very confident without knowing what the output is supposed to do. This is particularly true when introducing weights and smoothing. The current implementation has a lot of intermediate variables that are created—some do not even seem to be used anywhere—that make parsing the method even more complicated.
2. For smoothing, what is supposed to happen when there are missing values (no data in a region for a particular day). Similarly, how is smoothing handled when the smoothing window goes beyond the time-frame of the existing data. This is even more confusing when there are sample weights and I'm not sure what is supposed to be done.

[krivard]

The household CLI and ILI is computed using importance sampling; more information may be found here:

• Master COVIDcast methods doc
• Facebook methods doc (older)
• https://statweb.stanford.edu/~owen/mc/Ch-var-is.pdf; particularly 9.3: Importance sampling diagnostics, equation 9.13

For smoothing, in theory we're taking the weighted average of all data collected in the last seven days. If there is not at least one day of raw data in the last 7, we don't report a smoothed value for that day.

[taylor-arnold]

Thanks. Trying to read through these docs, but I'm still somewhat confused though (sorry!). I actually find the documentation you sent fairly straightforward, but I'm not sure how it completely explains the code.

Specifically, there is a lot of code that seems to be focused on something related to "mixing". For example, in aggregations-setup.R on lines 161-187. I can't figure out where this switch statement is described in the documentation or what I should expect it to do.

Searching through the three documents you linked did not surface anything with a similar name to "Mixing Coefficient".

Thanks again for dealing with all my questions!

[krivard]

Ah right, mixing -- sorry about that. The mixing code ensures that if some household is responsible for more than 1/100 of the score for the whole dayXregion bucket for whatever reason, then for small dayXregions (n<=100) we replace the normalized weight with 1/n, and for larger ones we use a number in between 1/n and the normalized weight, where "in between" depends on exactly how large the dayXregion is versus the maximum normalized weight in the dayXregion.

This should 100% be added to the signals explanation in the PDF -- @ryantibs @brookslogan I know documentation is no one's priority but this is the 3rd time I've had to explain this part of the facebook process (to umd twice, to lester, and now to taylor) and my stats knowledge is decidedly underpowered for the job. Can one of you either update the signals doc, or coordinate with me to schedule a call and I'll take dictation?

[taylor-arnold]

Great, glad I was not going crazy.

Currently this is not implemented for the Facebook Community surveys. That makes sense: we were previously not creating weighted versions and it's just a binary response, so no household can be more than 1/100 of the score if there are more than 100 households (the current reporting cut-off). I was asked to build in a weighted version of the community survey. Should mixing be implemented there as well?

[krivard]

Oof, probably yes, though I'm not sure if the change in statistic means any of the mixing thresholds need to change.

[ryantibs]

I wasn't even aware of this! @brookslogan It'd be ideal if Logan could just directly put a new part in signal_descriptions.tex to explain this.

[taylor-arnold]

I just pushed edits to the refactored facebook code. I believe that it now covers everything other than megacounties. Though as mentioned in my other notes, these will need to be carefully checked and possible debugged by someone with access to the full dataset.

@ryantibs @brookslogan — Going through the documentation, there are actually a number of differences between / details missing in the current documentation for the facebook signal. For example:

• The documentation describes a standard frequentist estimator of the community ILI signal, but the code implements a Bayesian estimator with a Beta(0.5, 0.5) prior. The standard error is, however, still in terms of the frequentist estimator.
• The code uses Jeffrey's correction to the standard error of the HH ILI and CLI signals, but this is not documented anywhere, or justified in any code comments.
• The documentation doesn't mention the "weight_in_geography" weights and how they are treated. This is confusing because at some points in the code there can be 4 or more different sets of weights all mixed together.
• The documentation does not document how sample size is computed. With all the different types of weight floating around, this get quite complicated.
• As mentioned above, there's no description of the mixing logic to handle large inverse probability weights.