COVID-19 transmission dynamics underlying epidemic waves in Kenya

SARS-CoV-2: To have or to have not In June 2021, official records in Kenya showed fewer than 4000 confirmed deaths and 180,000 confirmed cases of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). These data tend to reflect the economically advantaged strata of society who can afford smartphones and have access to medical attention and tests. Brand et al. developed an epidemiological model to estimate the impact of the pandemic in Kenya, the population of which was split into two socioeconomic strata. The authors predicted that 75% of the Kenyan population (about 39 million people) had been exposed to the virus by June 2021. If a fourth wave of infection is observed in the future, it would likely be driven by a variant with enhanced transmissibility or natural immune escape. —CA


Data description
In this paper, we make inferences of the penetration of SARS-CoV-2 transmission into each Kenyan county using a mechanistic transmission model. Joint posterior distributions for the parameters of the transmission model were inferred for each county using a synthesis of three data sources: • Kenyan Ministry of Health National linelist. We were provided linelist information about confirmed cases and swab tests performed. However, the available data differed by date and county. Below we describe which bits of data were available at different times: ○ Confirmed cases. PCR positive swab samples by collection date, ( 1 + ) , denoting the number of positive test results collected on day n for each county,, using the lab confirmation date of the swab test. Negative tests were not recorded in the national linelist. Those confirmed cases who died were also recorded by date of death for each county, ( + ) . The national linelist case data was available over the period we used for parameter inference in this paper (13 th March to 26 th April 2021), and over the period we used for out-of-sample validation of the model's predictive ability (27 th April -21 st May 2021). We screened out tests due to contact tracing or at entry to Kenya. ○ National combined laboratory test results. The combined reported tests, by laboratory collection date, across Kenyan laboratories, ( 3 − ) , , respectively denote the positive and negative test results by collection date in coastal county c on day n that were confirmed at KWTRP. As with the national linelist we screened out tests due to contact tracing or at entry to Kenya. The KWTRP linelist data was available from 13 th March onwards. The counties covered by the KWTRP linelist were Kilifi, Kwale, Lamu, Mombasa, Taita Taveta and Tana River.
• KWTRP serological surveillance programme (rounds 1 and 2). Numbers of seropositive ( + ) , and seronegative ( − ) , blood samples collected from regional centres of the Kenyan National Blood Transfusion Service (KNBTS) on day n originating from county c. Residual blood samples for serology were obtained from regular blood donors attending 4 regional KNBTS centres (Mombasa, Nairobi, Eldoret and Kisumu) in two rounds of collection in May -September 2020. The study methodology is fully described in Uyoga et al (2). We also show data from round 3 of the KWTRP serological surveillance programme (unpublished) in the main manuscript Fig.   2, which was collected January -March 2021, as validation data for the model predictive ability of seroprevalence.
• Google mobility data. Daily estimates of relative human mobility , compared to a baseline of the same date in the previous year (2019) derived from Google mobility trends (32). We assumed that changes in trends in SARS-CoV-2 transmission in Kenya were due to changes in the underlying population mobility. In particular, by changing frequency of indoor congregations. Therefore, we calculated , as the average change in baseline mobility over the "retail and recreation", "grocery and pharmacy", "transit stations", and, "workplaces" settings (Google defined categories), and also over the week prior to day n, in order to average over weekend effects. Due to incomplete data, and the likely bias introduced by using a mobility estimate derived from smartphone users in predicting the mobility of semi-rural populations outside of the major urban conurbations in Kenya, we consider only three areas: Nairobi, Mombasa and the pan-Kenyan aggregate (Fig. S1).
We performed a data cleansing process on the PCR datasets to improve their quality regarding the test date, age of the individual, symptoms and location. The data contained 4 spatial attributes; county, sub-county, ward and village, which we correct for inaccuracies in line with a ward level spatial baseline which accounts for the 47 Kenyan counties, 295 sub-counties and 1450 wards (https://data.humdata.org/dataset/administrative-wards-in-kenya-1450). This process includes spelling checks, string distance calculations, and an automated geo-search of addresses constructed using different combinations of the available spatial data per record. Finally, we run a distributed script that allows collating the number of positive/negative tests per day, county, sub-county, and age.
For making inference about unknown parameters for each county (see below) we used different sources of PCR swab testing data depending on the county and date according to the following rules: • For coastal counties the KWTRP linelist was used.
• For non-coastal counties the national case linelist was used between 13 th March and 6 th July, negative test results are unavailable over this period. • For non-coastal counties the national laboratory test linelist was used after 7 th July and 21 st May 2021. During this period both positive and negative test results were available.
The combined linelist data used in this paper is provided (Data S2), as well as the daily serology samples (Data S3).

Transmission model
We model the spread of SARS-CoV-2 in each of the 47 Kenyan counties as a two group SEIRS transmission process, with waning immunity returning completely protected recovered individuals to a waned immunity partially protected state (W). We attempt to capture the basic features of transmission between households, that is that we focus on the reduction in mixing indoors settings outside the home (termed places-of-interest/POIs) due to changes in behaviour in response to Kenyan government measures and an increased sense of personal risk. This approach is conceptually similar to some modelling studies in high-income countries (e.g. (33)), however, we used a simpler model structure reflecting the lower resolution data available in Kenya compared to most HICs.
We hypothesize that social disaggregation occurred in Kenyan counties; that is that there were broadly two groups of people in each county who differed in their ability to reduce social mobility. These two groups represented the people with lower and higher socio-economic status (SES) in each county. Both the higher and lower SES groups are assumed to have similar epidemiological characteristics except that we assumed that: • The lower and higher SES groups responded to Kenyan government measures by reducing their time spent in POIs at different rates. • The higher SES group reduction in time spent in POIs was well estimated by Google mobility trend data (32).
Alternatively, we can think of the two groups as being defined as those people whose visitation rates to POIs are well represented by Google mobility data (an unknown sized group in each county), and since this requires ownership of a smartphone, estimated as 30% of Kenyan adults in 2017 (34), this definition selects for being in a higher socio-economic status.
We assumed that transmission is sustained by mixing at POIs and, therefore, if the SES groups reduced their mixing in POIs by proportions, respectively, ( ) and ( ), at each time t then the reproductive number for each county was reduced from some baseline reproductive ratio 0 that would have be realised if social mobility had remained unchanged. By assumption, ( ) was fitted to Google mobility data, whereas ( ) was inferred jointly with other model parameters (see below). The per capita force of infection on individuals in the lower and higher SES groups, denoted respectively ( ) and ( ), was Where was the recovery rate, was the coupling factor between the two groups due to mixing in the same POIs (35,36), and, = , where N was the population size of the county, and ∈ [0,1] was the proportion of the population in the lower SES group. The full set of ODEs for each county was, Where is the reduction in susceptibility following a primary infection and is the waning rate of natural immunity.
The S,E,I,R disease compartments correspond to numbers of susceptible/naïve, exposed, infectious and recovered/immune (37), with subscripts denoting lower or higher SES group. The effective R numbers for an infected in either SES group were, Fitting the incubation and recovery rates. The incubation rate and recovery rates were chosen so that 1) The exponential growth rate r to basic reproductive number R0 relationship in this model matches that implied by Ferretti et al (38). The relationship was calculated using the formula 0 = �( ) −1 , where � was the moment generating function (mgf) of the generation distribution for the transmission model (39). Ferretti et al derived a best-fit generation time distribution for SARS-CoV-2 transmission: Weibull(2.8,5.7), and the mgf of the generation time distribution implied by (2) can be solved directly by using Laplace transforms. The importance of matching to r to R0 relationship was that r is the directly observed quantity from case data, whereas R0 scales directly with reduction in mobility. Therefore, in this study, which is both data-driven and estimates the underlying mobility of the lower SES group, matching the r to R0 relationship to an inferred SARS-CoV-2 generation time distribution was important.
2) The mean period between infections and symptom on-set was 5.1 days (40), assuming a mean pre-symptomatic period of 2 days (38, 41).
By combining these two criteria we chose incubation rate as = 1/3.1 days -1 , and the recovery rate = 2.4 days -1 , see Fig. S2 for agreement between the r and R0 relationship for this model and the generation time distribution given by Ferretti et al.
Waning immunity. The W disease compartment corresponds to the number of people who have had waning immunity following a natural infection of SARS-CoV-2. Full immunity is assumed to be lost at a rate =1/180 days -1 , and after loss of full immunity the W-group individuals have a decreased susceptibility to SARS-CoV-2 infection of = 0.16 relative to S-group individuals. There is substantial uncertainty about reinfection, and how transmissible reinfected individuals might be, despite a series of recent studies. We chose the reinfection rate, and subsequent protection from reinfection so that we broadly match these recent studies (Fig. S3). However, we couldn't find any data on whether reinfected individuals are typically as transmissible as people undergoing their first episode of SARS-CoV-2 infection. We defaulted to the maximalist assumption that reinfected individuals are as infectious as the typical person undergoing their primary episode.
Fitting contact rates for visiting POIs for lower and upper SES groups. A key assumption in this paper is that we assumed that there exists a group of people in each Kenyan county who are well described by Google mobility data (32). The Google mobility data tracks time spent by users/visitation rate in various settings relative to immediately before the pandemic. We took a fairly simple approach to using this information, we considered the averaged relativity mobility to indoors settings: "retail and recreation", "grocery and pharmacy", "transit stations", and, "workplaces", which we weighted equally. This gave a relative mobility rate to "risky" settings on each day after 20 th February 2020, which we assumed was data for fitting ( ). We fitted ( ) using a simple piecewise linear functional form: ( ) = 1 for ≤ 15 th March 2020 and ≥ 1 st November 2020, between these dates we fitted a rapid linear decline in ( ) from 15 th March to a minimum point on 15 th April 2020 (15 ℎ ) = 0.56, and then slower linear growth until 1 st November 2020 (Fig S1).
We assumed that the lower SES group followed the same pattern as the Google mobility data: decreasing from 15 th March 2020 until 15 th April 2020 then increasing back to baseline on 1 st November 2020 except that the minimum point of the lower SES mobility was treated as an unknown parameter , . The reason we chose Feb 1 st 2021 as the possible change point for increased transmissibility was based on the increase in frequency of alpha and beta variants in Kenya (20), which went above 50% of reported samples in Feb 2021 after low frequency in December 2020 and January 2020. No variants of concern were detected in Kenya before December 2020 (20).
The cumulative infection processes. The number of people who would test positive, either as PCR positive, or as seropositive, on each day n depended on: 1) the rate of new incidence on each day s < n, and, 2) the probability that someone who was infected on day s is detectable by either PCR or serology respectively = − days later. We calculated the daily incidence rate from the model by including the cumulative infections among the lower and upper SES groups and separating between first infection and all infections (including reinfections) as dummy degrees of freedom in the ODE system Where ( ) ⋅ and ⋅ ( ) were respectively the cumulative first infections, and all infections, by time t in each SES group. The daily numbers of new infections on each day n among each group and differentiating between first and all infections predicted by the transmission model was, , , = ( + 1) − ( ) for each day n.
And similarly for other combinations of SES group and infection episode.

Observation model
The underlying transmission of SARS-CoV-2 is not observed, rather we have access to swab tests and serological tests (positive and negative) aggregated by date and county. Therefore, we developed an observation model that connects unobserved daily transmission rates, which depend on the unknown transmission parameters, to the number of individuals in each county who were detectable as infected. Then we define parametric distributions for the daily pattern of positive and negative swab and serological tests on each day.
There was substantial day-to-day variation in both reported numbers of positive swab tests and percentage of positive tests among all samples confirmed that day (where negative test data is available). The underlying causes of the high level of day-to-day volatility are probable multiple including variation in daily testing rate, as well as in the settings at which swab test were collected, e.g. at the hospital, from a walk-up testing facility etc, as the focus of Kenyan public health teams has shifted over the course of the epidemic. Because of the substantial day-to-day variation, we use the standard robust alternatives to the natural Poisson and Binomial models for count data, the Negative-binomial and Beta-Binomial models respectively (42). Moreover, we assumed that the two SES groups had different parameters with respect to the chance of an infection being detected.
Observable infection status in each county. The probability that an infected individual would be determined as having been infected days after infection if tested by either a PCR swab test or a serology test was denoted, respectively, ( ) and ( ). By combining the underlying infection processes and the delay between infection and observability in our available data sets we find that the number of people who would test positive on each day in each county from either SES group, with a PCR test ( + ) , / , and/or a serology test ( + ) , / , is, where was the midpoint of day n. (6) was the false positive rate for the serology assay (see table S1). Underlying equation (6) is an assumption that the PCR test is 100% specific to SARS-CoV-2. Note that we are assuming that only the first infection contributes to the serological status of individuals, but that reinfections contribute to PCR status equally to first infections. Fitting ( ) . We fitted the sensitivity of a PCR swab test on each day s post-symptoms, ( ), to data on diagnostic accuracy given in Zhou et al (43), the fitted functional form for PCR-detectability more than 5 days after infection was: Where was the tail distribution function of a Gamma distribution with fitted shape parameter � = 18.4 and fitted scale parameter � = 1.1. This aligns with Zhou et al that the median period to become PCR undetectable after symptoms was 20 days with reported interquartile range of 17-24 days (43).
doesn't account for the delay between infection and becoming PCR detectable. To account for this delay we assumed that the distribution of delay between infection and maximum detectability followed the same distribution as the delay between infection and onset of symptoms (among those infected individuals who present with COVID-19 symptoms) as reported by Lauer et at (40), Where ( ) is the probability of developing symptoms on day s after infection. The true maximum sensitivity of the PCR test in a typical Kenyan setting is absorbed into our observation model via detection probability parameters (see below).
is displayed in Fig. S3.
Fitting ( ). The lag between symptoms and maximum detectability by serological assays has been reported as 21 days in a metastudy of reported diagnostic sensitivities (44). We assumed that, given an additional 5 day lag after infection (mean of onset time distribution), the sensitivity of the serological assay increased linearly from 0 at infection to a maximum of 82.5% (2) over 26 days.
is displayed in Fig. S3. Note that this implementation does not include the possibility of waning levels of detectable antibodies, that is seroreversion (23). For comparison in main Fig. 2 we include a forecast of seroprevalence using a simple model of seroreversion, Where is the daily probability of seroreversion if this only occurs on days after the 26 th day post-infection.
Detection of cases and number of swab tests performed each day. Obviously, our observations depend on the number of swab tests being performed on each day. Because of the fluctuations in testing rate in Kenya, and because of difficulty in asserting the reason for each swab test, we fit to the proportion of positive swab tests on each day whenever positive and negative test results are available (see below). However, we also aim to capture the true underlying detection rate over time; that is the % of all infections who are identified as a case. This detection rate changes with time due to the increasing and decreasing availability of tests, however, the number of tests performed on each day was also somewhat dependent on the demand for tests. Therefore, we expected the number of tests on each day to be correlated with the proportion of tests returning positive, because in periods with high PCR positivity there was also likely to be higher demand for tests to be performed.
We found that the daily number of tests performed across Kenya increased almost monotonically from the beginning of March 2020 until reaching 4000 tests nationwide per day in early July 2020, and afterwards was correlated with the proportion positive over all the tests, as expected (Fig. S4). We interpreted this finding as the testing regime being composed of two main periods: 1) March -June 2020, under-capacity of testing in Kenya when infections were more likely to be unidentified due to simply lacking available swabs, and, 2) July 2020 onwards, Kenyan testing is at maximum capacity and the testing rate responds to demand. We give a simple piece-wise linear form for the relative detection rate (RDR) in Kenya on day n, Negative-binomial model for number of daily positive swab tests. The mean detection rate per PCR-detectable individual per day by swab testing for each SES group ( , and , ), and the clustering factor 1 of the daily detections ( ), which we point estimate as = 0.5. The observed number of positive swab tests in each county on day n was connected to the underlying detection rate using a negative binomial distribution, Where is the mean number of PCR positives expected by the model accounting for the detection rates by SES group, relative detection rate (which reduces the chance of detecting a case before July 2020) Beta-binomial model for proportion of daily positive swab tests. We didn't expect the daily swab test results to reflect an unbiased sample of the underlying PCR-detectable fraction of the population. Therefore, we include a relative bias parameter for a PCR-detectable individual being tested relative to a PCR-undetectable individual, for each SES group ( , ), and an effective sample size parameter ( ; (42)). We used a point estimate of = 30. On days where negative swab tests were available, we connect the observable status of the epidemic to the data thus, Where , is the total number of PCR swab samples collected on day n in the county being fitted and is the proportion of tests performed returning positive expected by the model, accounting for bias in the sampling regime and the possibility that tests occur among the upper SES group with probability independently of the underlying numbers of PCR detectable individuals. Beta-binomial model for proportion of daily positive serological tests. The reported uncertainty in the maximum sensitivity of serology assay was fairly high: the posterior mean sensitivity was 82.5% (credible interval 69.6-91.2%; (2). The posterior uncertainty in the serological sensitivity influenced the confidence the inference method placed on the serological sample data; if the test sensitivity was known to high precision we would treat each day's serological samples as a binomial draw from an underlying proportion of seropositive individuals given by equation (6). Given that the sensitivity of the serological assay was itself an uncertain factor we fitted the posterior uncertainty in the testing sensitivity to a beta distribution: ∼ ( � = 33.6,̂= 7.13). This implied that the appropriate observation model for the number of positive serological samples on day n (( + ) ), out of the total number of serological samples being collected on day n, , was a Beta-binomial distribution, Given an underlying realization of the transmission process the mean number of positive serological samples on day n is , ( + ) . The "total-count" parameter (42), = � + = 40.73, allowed for greater dispersion in the observed seropositive count data than would be allowed by a Binomial model.

Parameter inference
As described in Data description section we had access to data on daily swab tests (positive and negative) and the KWTRP serological survey daily samples. We grouped the parameters for inference into transmission model parameters ( ) and observation model parameters ( ) and used Bayesian likelihood-based inference to infer parameters. A challenge with using the linelist data in Kenya for inference of transmission was that the metadata concerning the reason for receiving a swab test, the levels of symptoms of people who tested positive, and their healthcare outcomes were often missing. Overall, more than 90% of the people who tested positive in Kenya, and for whom we have a description of their symptoms, reported no symptoms (asymptomatic). Therefore, unlike model-based inference for COVID-19 transmission in high-income countries we didn't use severe outcomes such as hospitalisation or death as data sources for inference, e.g. (7), because this data was unreliable. Instead, we concentrated on fitting to the proportion positive of daily swabs test and serological tests jointly with detection rate of cases.
We use the Bayesian inference to infer a joint posterior distribution for the unknown parameters (both transmission-based parameters and observation-based parameters) for each county. We describe the three main ingredients for our Bayesian approach below: 1) the log-likelihood function for the data given a set of parameters, 2) the county-specific hierarchy of prior distributions for the parameters, and, 3) the Markov-chain Monte Carlo method used to draw parameter sets from the posterior distribution.
Log-likelihood function. The observation model gives the following log-likelihood function for the unknown parameters = ( , ) given the sampling data for a county, Where and are, respectively, the probability mass functions for the negative binomial and beta-binomial distributions, as described in the observation model subsection. The index set of days in covered days where PCR negative tests were available in the county being fitted. The first day where samples were included in the log-likelihood calculation was 12th April 2020, due to testing being even more irregular before that date, and the last day was 27 th April 2021.
County-specific hierarchy of priors. As described in the main text, when designing priors for Bayesian inference we, first, grouped the counties by socio-economic and epidemiological factors using an unsupervised machine learning technique. The socio-economic/epidemiological factors are described and were generated in Macharia et al (28), and included factors such as county rates of obesity, informal settlement habitation, population density and access to major cities. We used ordered leaf clustering by centred and standardized 1 (Manhattan) distance over all factors to group the 47 counties and identified four groups of similar counties (Fig S5). Nairobi and Mombasa were sufficiently distinct from other counties that they formed their own singlet groupings. There were two further county groupings evident from the unsupervised learning, which corresponded to predominantly remote and rural counties and to counties with partial urbanization and greater connectivity to the main Kenyan cities. The full county list of the identified groups was: 1) The capital city Nairobi (also a county) and, 2) The second city Mombasa (also a county) 3) Semi-urban/semi-rural counties that either contain significant sized cities and/or neighbour Nairobi county: We also distinguished the counties by the number of serological tests in the KWTRP serological surveillance trial rounds 1 and 2 that had been performed on people living in them. Serology test data was critical for identifying parameters involved in the detection of infections in each county by giving an estimate of the true population exposure at different time points. 11 counties had notably more serological tests than the 36 other counties: Embu, Kilifi, Kisii, Kisumu, Kwale, Mombasa, Nakuru, Nairobi, Nyeri, Siaya, Uasin Gishu (Fig S5).
For the 11 counties with higher numbers of serological tests we assigned priors that encoded our a priori beliefs about the epidemic in Kenya. The prior distributions for ϵ, , , , were the same for each county reflecting our a priori beliefs. We had 1) strong confidence that mixing between SES groups would be assortative, 2) weak confidence that the lower SES group would be able to reduce their mobility less than the upper SES group, 3) moderate confidence that the daily detected proportion PCR tests being positive would be higher than the actual proportion of the population PCR detectable, with this bias being higher among the lower SES group who were believed would be less likely to seek a test when asymptomatic compared to individuals in the upper SES group, and, 4) moderate confidence that new variants dominating transmission from February 2021 onwards were about 50% more transmissible: The county group specific priors for the 11 counties with higher numbers of serological tests were: • Nairobi and Mombasa.   The county group specific priors were based on the view that although the most a priori likely possibility was that most counties had very few infected individuals on 21 st February, and that we might expect the unknown numbers of exposed people to be concentrated in cities. We had moderate a priori confidence that the overall detection rate among the upper SES would be around 5% of infections, and about 1% among the lower SES group. We had strong confidence that the inferred community size of the upper SES group would be higher in cities rather than other counties.
For the other 36 counties with less serological testing, we used the same county-group specific priors except that the PCR positivity bias parameters ( , ) and daily detection rates of positives ( , , , ) had priors derived from the inferred posteriors of counties among the 11 with higher numbers of serology tests. This was an approximation to a hierarchical Bayesian formulation of a joint likelihood for all 47 Kenyan counties; we fitted gamma distributions to the pooled posterior draws for these parameters over the semi-urban/semi-rural counties with larger numbers of serological tests (Embu, Kisii, Kisumu, Nakuru, Nyeri, Siaya, Uasin Gishu) and the rural counties with larger numbers of serological tests (Kilifi and Kwale). These fitted priors were used for parameter inference.
for each county using the NUTS-HMC sampler implemented by the Julia language package dynamicHMC.jl. Solving the likelihood function for a proposed value of involved solving the ODE system (2), we used the highly performant DifferentialEquations.jl package for ODE solutions (45). The HMC method required a log-likelihood gradient, , which, for our use-case of a small ODE system with a low number of parameters, was most efficiently supplied by forward-mode automatic differentiation implemented by the package ForwardDiff.jl.
The MCMC chain converged for each county (all MCMC chains and MCMC diagnostics can be accessed through the linked open code repository). The posterior mean (and 95% CIs) for each parameter can be found in supplementary data: Data S1

Back-calculation and forecasting for Kenyan counties
The parameter draws from the posterior distribution defined the uncertainty in our model nowcasts and forecasts for each county, since the underlying transmission model was deterministic. Therefore, posterior distributions for epidemic quantities were created by simulating the epidemic for each ( ) , = 1,2, … posterior draws. We used this technique for both back-calculating estimates of quantities that were unobserved during the data collection period, such as the daily infection rate in each SES group and the seropositivity rates after round 2 of the serological survey finished, and, forecasting ahead of the data collection period by assuming that the parameters remained unchanged in May-July 2021.
When forecasting the number of positive tests that would occur on days in the future, we combined a forecast of the proportion positive expected on each day, for each posterior draw of a parameter set, with a forecast of the number of tests that would be performed on each day. We noted that the number of tests performed was not independent of the proportion positive, and for each county fitted a simple linear model For days that cover the most recent 60 days of available test data ( are iid normal errors).

Inference of observed mortality rate per infection
The commonly used infection fatality ratio (IFR) by age estimates from Verity et al (46), weighted by the Kenyan population distribution given by the 2019 census, implied a basic IFR prediction in Kenya of = 0.264%; that is 264 deaths per 10,000 infections. This assumed a uniform attack rate across age-groups in Kenya.
We used the posterior predictions of the underlying daily infections in Kenya counties to infer a crude infection fatality ratio ( ) for each Kenyan county. The lag between infection and death, for those infected individuals who die, was defined as the convolution of three time duration distributions: 1. The duration of time between infection and symptoms (days), which we assumed was distributed ( � = 1.64, � = 0.36) (40). 2. The duration of time between initial symptoms and severe symptoms (days), sufficient to seek hospitalisation, which we assumed was distributed (1,5) (47). 3. The duration of time between severe symptoms and death estimated from UK hospital data (7). This was an empirical distribution with mass function .
We discretized the first two distributions to give probability mass functions for the number of days between infection and symptoms, and for the number of days between symptom onset and severe symptom onset. The probability mass function for the (discrete) number of days between infection and death, for those who died, , was given as a discrete convolution over these probability mass functions: ( ) = * * ( ) for the probability that death occurs days after infection.
We didn't use mortality data for inference of parameters; therefore, we estimate mortality detection rates from back-calculation of the underlying daily incidence rate. Because of the possibility of substantial under-reporting of mortality due to COVID-19 in Kenya, which may differ across socio-economic groups, we jointly fit a mortality-detection rate to both SES groups in each county = ( , ). We used maximum likelihood estimation to infer mortalitydetection rates: Where ( ) were the number of deaths attributed to COVID-19 on each day t, � and ��� were the posterior mean back-calculations for the daily rate of new infections among the lower and upper SES groups, which were convolved with a probability kernel for the lag in days between infection and death.
It should be noted that because the under-reporting rate of deaths in Kenya is unknown, the county-and-SES-group specific mortality parameters estimated should be interpreted as the IFR scaled by the chance of the death being reported.

Further evidence to support the two social group model
Comparison to previous modelling of the first wave in Kenya. A previous model for the first wave in Kenya used a similar methodology, but rather than using two SES groups each county was modelled as forming a single transmission group with an effective population size to account for heterogeneity in transmission and with R(t) fitted daily rather than assuming the simple parametric form used in this paper (6). The one-group model was successful in capturing the dynamics of the first wave in Kenya, and is close agreement with the model presented in this study over the course of the first wave in Kenya, for example, the one-group model predicted 43.3% (CI 35.3%-49.5%) by 30 th September (6) whereas the posterior mean estimate for overall exposure in Nairobi for the two group model was 40.0%. However, the one-group model failed to explain the second wave in Kenya without large increasing in R(t) relative to the February 2020 estimates, with these increases starting months before VOCs like B.1.1.7 were first detected in Kenya (30, Fig S6). These large increases in R(t) for the one-group model derived from the model having a single detection rate for all infecteds, whereas the two-group model favours an explanation where infections in the lower SES group are detected at lower rates than infections in the upper SES group.
Randomised seroprevalence survey in Nairobi. The model that we present in our current paper proposes separate SES groups with infections in the upper SES group being both delayed and more densely sampled compared to the lower SES group. An obvious alternative explanation is that the serology data used to fit the model was an over-estimate of exposure, possibly due to the underlying blood donor serology data being heavily biased in favour of the detection of seropositive individuals. However, the overall model prediction of seropositivity in Nairobi was in reasonable agreement with the overall seroprevalence random selection of the Nairobi population performed in the first half of November 2020 (48). The agreement between the model prediction and the randomized serological surveys was even better for the central age groups (ages 10-60), whereas the overall seropositivity (including <10 year olds and > 60 year olds) was lower. This may reflect model bias in that we calibrate to blood donor serology, which limits to donors in the 15-65 age interval. However, first, the predictive error was not substantial (< 10% seropositivity prediction error), and, second, it is not clear if this bias will persist after general reopening of schools in Kenya (Fig S7).
The randomised seroprevalence survey also reported a breakdown by subcounty (48). For comparison to model predictions we compared some exemplar sub-counties to predictions for lower and SES group model predictions of seroprevalence. The Kibera sub-county, containing the Kibera slum area, ranked highest out of 17 Nairobi sub-counties for both poverty and density of informal settlement and we used the seroprevalence of this county as an exemplar of the lower SES group. No Nairobi sub-county was lowest ranked for both poverty and informal settlement density so as the upper SES group exemplar we used the seroprevalence over three sub-counties: Dagoretti South (12 th for poverty, 15 th -17 th for informal settlement), Embakasi East (15 th -17 th for poverty, 12 th for informal settlement) and Roysambu (14 th for poverty, 13 th for informal settlement). Poverty and informal settlement density data was sourced from Macharia et al (27). The agreement between the exemplar sub-county unadjusted seroprevalences and the model predictions were also reasonable, which suggested that segregating the population by SES status was a reasonable approach to capturing the wide disparity in seroprevalence between Nairobi sub-counties reported in (48).
As well as agreeing on the sero-prevalence, the estimated mortality-detection rate for Nairobi (c.f. Inference of observed mortality rate per infection, and Data S1) was 0.03%, fitted on reported deaths over the entire period. This is in reasonable agreement with the IFR reported by the November seroprevalence survey, 0.04% (48).
Data on public and private hospital admissions. We had imperfect data on hospital admissions in Nairobi, derived from the Kenyan national linelist, and did not use admissions to fit the model. However, in the Kenyan capital Nairobi we had a subset of cases with PCR swab test confirmation and with metadata on admission to a named health facility. During the first wave we note that most of the identified admissions occurred in public hospitals, which would be selected by people in the lower SES group, whereas in the second wave private health facilities were more frequently reporting admissions (Fig S8).

Formal model selection and sensitivity analysis
Model selection using Deviance Information Criterion. We considered two alternative transmission models: • A one-group transmission model, i.e. homogeneous mixing within the county, with an effective contact rate (relative to pre-pandemic baseline) ( ) which was fixed as ( ) = 1 when < 15 th March 2020, and then allowed to vary daily. Otherwise, all model dynamics were like equations (1)(2)(3)(4)(5), albeit with one group and, therefore, less parameters (apart from the daily ( ) values) to infer both for the transmission model and the observation model. We did not infer a school closure effect for the one-group model, since this is absorbed into estimation of ( ). • A three-group transmission model with fixed population sizes representing the bottom 25% of income, middle 50% of income and top 25% of income in the population. The effective contact rates for each group were assumed to have the same trend as Google mobility data (decreasing from 15 th March 2020 to 15 th April 2020, increasing subsequently until return to pre-pandemic baseline on 1 st November 2020 Fig. S1), however, we did not fix any group to have a 45% decrease in mobility (as per the main model). Otherwise, all model dynamics were like equations (1)(2)(3)(4)(5), albeit with three groups and, therefore, more parameters to infer. As in the main model, the fraction of contacts made within group was defined as a parameter ( ), with outside group contacts made pro-rata according to the outside group size.
For both alternate models we compared to the two-group model using data from Nairobi, due to the higher rates of testing in the Kenyan capital. The three-group transmission model, after adjusting the log-likelihood function to accept three sets of detection rate parameters, had parameter inference using Hamiltonian MCMC (HMC) in a similar manner to above. The onegroup model required inference on each daily ( ) value as well as inferring the other model parameters, e.g. 0 , (0) etc. To infer the ( ) values for the one-group model we used the expectation-minimisation (EM) algorithm, alternatively using HMC to infer posterior distribution model parameters with ( ) fixed (E-step), and converging all daily ( ) values to their maximum value with the posterior distribution of other model parameters fixed (M-step); see for (6) details. The one-group model here is similar, but not identical, to the model used to fit the first wave in Kenya (6) and presented above (Fig. S6). The main difference is that the model used in (6) had effectively two groups, an effective population size of individuals at risk of infection (less than the actual population size) was inferred. Moreover, we now have a longer time series of serology data for use in parameter inference.
The three models were compared using the Deviance information criterion (DIC), as formulated in (42), Where, (⋅) and (⋅), were respectively the posterior mean and posterior variance over the sets of parameters drawn from the HMC process. DIC is known to have an asymptotically (large sample size limit) chi-squared distribution (42) therefore, even away from the large sample size regime a difference in DIC of greater than 15, Δ > 15, should be considered very substantial statistical evidence towards favouring the model with lower DIC.
In both cases, a comparison between the two-group model and the one-group model (Δ = 5375.4) and between the two-group model and the three-group model (Δ = 66.9), found substantial statistical evidence to favour the model used in the main text for Nairobi.
In particular, the one-group model was unable to form a coherent explanation of both the PCR test data and the serological test data in Nairobi. The essential problem for the one-group model is that the seropositivity rate grows differently to the case detection rate; the most coherent solution found by the one-group model was very rapid spread before measures were implemented around 15 th March 2020, and before there were many PCR tests in Kenya, then variable transmission rates over 2020. This agrees with the overall trend of observed PCR cases in Nairobi but fails to capture the increasing seropositivity rate. On the other hand, the two-group model can coherently fit both data streams by allowing PCR-positive individuals in the lower SES group to be detected at a lower rate than higher SES group individuals (Fig S9). Note that as well as the two-group model being heavily favoured by DIC, visually it predicts the round 3 KNBTS serology data in Nairobi (which was kept out-of-sample and not used in inference) much better than the one-group model (Fig S9). Sensitivity analysis. In addition to formal model comparison, which was performed under a baseline set of fixed assumptions that: (i) the period of complete immunity after first infection episode was on average 180 days, (ii) the long-term susceptibility after loss of complete immunity was 16% compared to naïve individuals, and (iii) the relative infectiousness of infected individuals during subsequent infection episodes compared to their first infection episode was 100% compared to first-time infected individuals, we also considered a range of sensitivity scenarios (see Table S2 for the full list of sensitivity scenarios).
For each scenario, we inferred joint posterior distributions for all free parameters using HMC as described above. For scenarios where we assumed that individuals in subsequent infection episodes were less infectious than their first, we modified the equations (1-2) so that the waned immunity state (W) flowed into new states exposed-after-waning (EW) and infectious-afterwaning (IW) for each SES group. Infected individuals in the IW state contributed less to the force of infection (1) according to their assumed relative infectiousness.
After parameter inference we found that in Nairobi inferring two substantial effective social groups, characterised by a strong preference for within-group mixing, and with one group (lower SES group) reducing their contact rates substantially less than the other (higher SES group) early in the pandemic was common across all sensitivity scenarios. The posterior mean estimators for assortative mixing ( ) was between 0.925-0.99 over all scenarios (Fig S10), breakdown of size of lower and higher SES groups was between 54%/46% -63%/37% over all scenarios (Fig S11), and the minimum (15 th April 2020) contact rates relative to pre-pandemic baseline for the lower SES group was between 82-88% over all scenarios (Fig S12).

Inferred patterns in transmission dynamics across Kenya
Most counties were inferred to either have comparatively low assortative mixing between groups ( < 0.9) and/or an effectively small higher SES group (<15% of the county population size), c.f. Data S1. For majority of counties, we infer an initial number of exposed individuals < 1 on February 21 st , 2020, which we interpret as the epidemic having not arrived in those counties by that date, but rather arriving later spread from neighbouring Counties. These counties were predominantly rural counties, containing most of the Kenyan population, and with R(t) similar to Mandera county, which we used as an exemplar in the main manuscript (figure 3 bottom right). For these counties, R(t) was low but persistently above 1 (R(t) ~ 1-1.5) until November 2020, when cases began to decline due to the depletion of susceptibles over that period. The epidemic was comparatively slow moving in the rural counties due to a combination of delayed epidemic arrival, inferring a basic R0 < 2 for that county, and/or inferring significant reduction in R(t) due to school closures.
The main cities in Kenya, Nairobi and Mombasa, and most counties surrounding Nairobi, were inferred to have a larger proportion in the higher SES group than most rural counties (>15% of the county population size), and to have highly assortative mixing ( > 0.9; Fig. S13). Additionally, Meru county and Garissa county, which contains the Dadaab refugee camp, were in this group of counties (Fig S13). Counties with an inferred higher proportion of higher SES group and a high degree of assortative mixing had a distinct first peak in June-August 2020 due to depletion of susceptibles in the lower SES group whilst individuals in the higher SES group were still at substantially lower mobility rates compared to pre-pandemic. Although only 28% of the Kenyan population live in urban areas, most tests have been performed in urban areas. Consequently, the two-wave pattern observed in urban areas dominates the overall Kenyan epidemic trajectory of detected cases. In addition to Fig S13, the county-specific data and model predictions for PCR positives, deaths, seropositivity, and R(t) can be found in the Github repository associated with this paper https://github.com/SamuelBrand1/kenya-covid-three-waves .

Supplementary Text
Notation for distributions used in this study In this study, we have used a number of parameter symbols that are also the most commonly used symbols for various common parametric distributions. Moreover, these parametric distributions are used in the underlying analysis frequently with their distribution parameters defined as functions of underlying transmission model states. To reduce misunderstanding reserve symbols with "hats" as referring to the parameters of a parametric distribution and use "=" to refer to the value of the parameter. Find below the choice of parametrization for the parametric distributions used in the study: • Exponential distribution.
(̂= ), with mean . to Google defined setting categories, "retail and recreation", "grocery and pharmacy", "transit stations", and, "workplaces". Also, shown is the 7-day moving average of the mean over these settings (black curve) and the approximation for the contact rate used for the upper SES group in this paper (red curve).          Mean latent (infected but not infectious) period 1/ 3.1 days. The mean incubation period (40) was reduced by two days of pre-symptomatic transmission (41) to give a latency period.
Mean infectious period 1/ 2.4 days. Chosen to fit the exponential growth rate r to reproductive number R relationship of a fitted SARS-CoV-2 generation duration (38) Baseline reproductive number, 0 .

Inferred from data
Population proportion in lower SES group, .

Inferred from data
Mobility of upper SES group, ( ).

Mobility of lower SES group, ( ).
Piecewise linear. ( ) = 1 for ≤ 15 th March 2020 and ≥ 1 st November 2020. Linear trend decreasing to, and then increasing from, ( ) which was inferred from data.
Proportion of secondary infections within same SES group, .
Relative susceptibility compared to naïve susceptibles after loss of complete protection after first infection episode, .
Relative effect on transmission of having schools closed, ℎ .

Inferred from data
Relative effect on transmission post 1 st February 2021 due to introduction of variants of concern, . Inferred from data

County-specific observation model parameters and data
Number of people who would test PCR positive on day n in lower/upper SES group, ( / + ) .
Dynamic Number of people who tested PCR positive on day n, ( + ) .

Data
Number of people who would test as seropositive on day n in lower/upper SES group, ( / + ) .
Dynamic Number of people who tested as seropositive on day n, ( + ) .

Data
Probability that an infected individual would test PCR positive on day after infection, ( ) ( ) = * ( ), where was the tail function of a gamma distribution fitted to data given in (44), and is the probability function of onset of symptoms post-infection (40).
Probability that an infected individual would be detectably seropositive on day after infection, ( ) ( ) is linearly increasing over 26 days to saturate at 82.5% sensitivity, based on report delay in seroconversion (44) and maximum sensitivity of serological assay (2).
Daily rate of PCR-positive individual in lower SES group receiving a swab test, , .

Inferred from data
Daily rate of PCR-positive individual in upper SES group receiving a swab test, , .

Inferred from data
Clustering coefficient of daily PCR tests, . Point estimate: = 0.5.
Relative bias in favour of selecting a PCR positive individual for swab testing among lower SES group (used when negative tests are available) ( )

Inferred from data
Relative bias in favour of selecting a PCR positive individual for swab testing among upper SES group (used when negative tests are available) ( )

Inferred from data
Proportion of daily swab tests among upper SES group (used when negative tests are available) ( )

Inferred from data
Effective sample size/clustering coefficient of daily PCR tests (used when negative tests are available), ( ). Point estimate: = 30.

County-specific detected fatality rate parameters
Chance of fatality and being determined as due to COVID-19 disease per infection in the lower and upper SES groups, / .

Inferred from data
Probability mass function of delay lag between infection and death for those who die, ( ).
Derived as a convolution over the lag from infection to symptom onset (40), the lag from symptoms to hospitalisation/severe symptoms (47), and the lag between severe symptoms and death (7). Table S1: Dynamic and observational model variables and parameters. "Dynamic", means that the variable was an output of the transmission and observation model for the county. Table S2.