FMSN30/MASM22

Use the official data from Folkhälsomyndigheten to investigate the change in the daily number of new confirmed cases over time in different Swedish Regions.

You can pick one region to study, or a few that you study separately. Or use a joint model for several regions at the same time. That will require interaction terms between region and whatever time variables you have in the model.

Discussion forum; Project 3A Covid-19 help

You should always

plot the data with the date on the x-axis together with the estimated mean, $LaTeX: \hat{\mu}_i$ , and its confidence interval. If you want to, you can add the bootstrap prediction interval as well.
If you want to compare regions using cases by 100000 inhabitants, you should still fit the model using the raw number of cases, with the population size as an offset variable, see, below. You should calculate and predict everything for the raw number of cases, before you transform it to cases by 100000 inhabitants.
report the rate ratio, $LaTeX: e^{\beta_1}$ , for the rate of change over time, together with its confidence interval.
plot the standardized deviance residuals against time, together with suitable reference lines and a comment on their behaviour.

Suitable model: Negative binomial regression

Since the number of positive tests depends on the size of the population you have to use an "offset" variable. This means that the expected value is directly proportional to the offset variable:

$LaTeX: \mu_i = \text{population}_i \cdot e^{\mathbf{x}_i\beta}$

and

$LaTeX: \ln \mu_i = \ln(\text{population}_i) + \mathbf{x}_i\beta$

The offset variable is specified with

glm.nb(new_cases ~ x1 + x2, offset(log(population)), data = fhm.data).

Note: strictly speaking it should ,of course, also depend on the number of tests that were performed, but we, unfortunately, don't have that information.

Data

The data frame, fhm.data, is located in fhm_data.RData.

It consists of the daily reported new confirmed positive cases, separately for each of the 21 Swedish regions, from Folkhälsomyndigheten.

The last day in the data set is 8 May, using the updated information from 11 May.

Note: Jönköping had the first confirmed case in Sweden in early February. The person stayed in quarantine and no new cases were found until three weeks later. I have deleted this first observation so the data set starts 24 February, two days before the second confirmed case.

The days are numbered consecutively, starting 24/2. Since different regions had their first case on different dates I have also numbered the days within each region, starting two days before the first case in the region. This starting date varies from 24/2 to 11/3.

obs_date	The date (in R date format). Can be used on the x-axis in plots
day_nbr	Consecutive day number 1 = 24/2-20
day_nbr_region	Consecutive day number within region. 1 = 2 days before the first case in the region
region	The name of the region. Use factor(region) to turn it from text into a factor variable if you want to relevel.
population	The population size
new_cases	: the number of new cenfirmed cases

Some ideas for what you could do

Choose any of the following problems. You should be able to talk for 10-15 minutes about what you have done.

1. The rate of the exponential growth

Model the rate of the initial exponential growth in the number of new cases in some region(s). In some of the regions the rate has stopped growing or is even declining. In that case, only use the first part of the data.

Some things you could chose to focus on:

(a) Late or early start

You could compare the increase rate in regions that had an early start with regions that had a later start to see if there is any systematic difference. Divide the regions into the categories "early" and "late" using the date of their first observation.

(b) Poisson of Negative binomial

Since we are counting the number of cases we might expect that a Poisson distribution would work. However, the data varies much more than a Poisson distribution allows, making a Negative binomial distribution fit much better. Fit both types of models and illustrate how the bad fit of the Poisson, and better fit of the Negative binomial, shows itself is, e.g., standardised deviance residuals and prediction intervals.

(c) The early tracing

In the early stages, the testing was mostly related to tracing the infection to/from a confirmed case. Later, this approach was abandoned and instead the testing was focused on patients admitted to hospital. This means that a number of the regions have a short peak in the cases in the first week or so, that then drops. Investigate how these outliers influence the estimated rate.

(d) The weekend effect

We might expect that the number of tests performed, and thus the number of new cases, varies over the week. Extract the day of the week by

data$dayofthweek <- weekdays(data$obs_date)

There were also a few days that you might want to treat as if they were Saturdays or Sundays: 10/4 (Good Friday), 13/4 (Easter Monday), 30/4 (Last of April), 1/5 (First of May).

Add the day of the week as a categorical variable to examine whether there are any systematic relative differences in the number of new cases reported on different days.

2. Rise and fall in the number of new cases

In the regions that had an early start, the number of new positive cases, new_cases, has started to decline again. Investigate how this affects the standardized deviance residuals if you fit a model that does not take this into account, i.e.,using only day_nbr (and interaction with region of you want to study several regions in the same model).

You could then try to find some way to model this rise and fall. Note: you do not have to succeed! Interesting failures are also allowed. Things you could try include, but are not limited to:

(a) Polynomial regression: add a square term or other polynomial (day_nbr + day_nbr^2 etc)

(b) Piece wise regression: Identify the day of the peak, day_top, (by hand and some help from geom_smooth()) and fit a piece wise model "Fall ~ day_rise + day_fall" with the new x-variables day_rise and day_fall calculated using:
data$after_fall <- as.numeric(data$day_nbr > day_top)
data$day_rise <- data$day_nbr*(1 - data$day_after) + day_top*data$day_after
data$day_fall <- data$day_nbr*data$day_after + day_top*(1 - data$day_after)

You can investigate the effect of chosing different day_top-values with a combination of residual plots and AIC. If you know about cross-validation, feel free to use it!

If you want to model several regions at once, you will need a separate day_top-value for each region, as well as interaction between region and each of day_rise and day_fall to allow for different slopes in different regions.

In order to better catch the rise-plateau-fall behaviour you might want to identify two or more time points, first when the rise stops and then when the decline starts. This requires the addition of a day_plateau for the days between these time points. Feel free to experiment.

In some regions it is unclear whether the number of cases has started to fall or not. Try to determine this by fitting both a linear and a non-linear model and comparing them,