I en artikel i SvD för några veckor påståds det att vi efter Corona-pandemins värsta dagar nu har en underdödlighet i Sverige. Det här kan enligt artikelförfattaren ses som en effekt av att många personer i riskgrupper gått bort under våren.
Jag fick en fråga av min kompis Benjamin som undrade om det här stämmer, givet data från SCB, och om det finns ett sätt för oss att slå fast huruvida vi haft en underdödlighet eller om det egentligen sett likadant ut.
Det här liknar i mångt och mycket ett klassiskt hypotestest, där vi har en nollhypotes \(H0\) att ingenting har hänt och \(H1\) att nollhypotesen förkastas, d.v.s. vi har haft en under eller överdödlighet.
För att besvara den här frågan har vi en Excel-fil med antal rapporterade dödsfall som SCB fått in under ett antal år. Vi börjar med att läsa in den.
library(tidyverse)
library(readxl)
data <- read_excel("scb-doda.xlsx", sheet = "Tabell 1", range = "A7:G373")
data## # A tibble: 366 x 7
## DagMånad `2015` `2016` `2017` `2018` `2019` `2020`
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 januari 295 245 339 275 300 254
## 2 2 januari 272 272 317 312 276 300
## 3 3 januari 281 293 310 295 271 256
## 4 4 januari 256 260 305 285 282 272
## 5 5 januari 268 260 304 291 275 265
## 6 6 januari 256 246 304 268 281 285
## 7 7 januari 309 270 313 272 280 246
## 8 8 januari 294 292 296 260 274 292
## 9 9 januari 272 276 288 274 239 277
## 10 10 januari 285 256 304 287 254 253
## # … with 356 more rowsFör att kunna arbeta med den här datan behöver vi samla alla observationer i en kolumn och skapa en datumkolumn.
library(lubridate)
library(tsibble)
tidy_data <- data %>%
pivot_longer(cols = `2015`:`2020`, names_to = "år", values_to = "antal_döda") %>%
mutate(datum = paste0(år,
str_extract(DagMånad, "[:alpha:]+"),
str_extract(DagMånad, "[:digit:]+")) %>%
as.Date(format = "%Y%B%d")) %>%
mutate(dag = day(datum), månad = month(datum)) %>%
select(datum, antal_döda) %>%
filter(datum < "2020-08-18") %>%
na.omit() %>%
as_tsibble(index = datum)
tidy_data## # A tsibble: 2,056 x 2 [1D]
## datum antal_döda
## <date> <dbl>
## 1 2015-01-01 295
## 2 2015-01-02 272
## 3 2015-01-03 281
## 4 2015-01-04 256
## 5 2015-01-05 268
## 6 2015-01-06 256
## 7 2015-01-07 309
## 8 2015-01-08 294
## 9 2015-01-09 272
## 10 2015-01-10 285
## # … with 2,046 more rowsVi kan nu visualisera det här över tid:
p <- tidy_data %>%
ggplot(aes(x = datum, y = antal_döda)) +
geom_line() +
theme_minimal() +
labs(title = "Antal rapporterade dödsfall över tid",
caption = "Källa: SCB")
plotly::ggplotly(p)För jämförelsens skull kan vi visualisera varje år som en separat linje.
library(feasts)
library(gghighlight)
tidy_data %>%
filter(datum < "2020-08-17") %>%
as_tsibble() %>%
gg_season() +
theme_minimal() +
labs(title = "Antal rapporterade dödsfall över tid",
caption = "Källa: SCB") +
gghighlight(max(antal_döda) > 350, label_key = år)
Vi kan undersöka distributionen under dessa perioder.
tidy_data %>%
filter(month(datum) %in% c(07)) %>%
mutate(år = as.factor(year(datum))) %>%
ggplot(aes(antal_döda, color = år)) +
geom_density() +
labs(
title = "Fördelning antal döda under juli"
) +
theme_minimal()
Vi kan testa att göra en STL dekomposition.
Hur exakt man ska testa hypotesen för det här inte självklart. En väg framåt skulle vara någon form av bootstrapping. Men det är ett projekt i sig att sätta upp en sån struktur och kräver nog ganska gedigen kunskap i att analysera tidsserier.
Ett enkelt alternativ är att använda en linjär modell endast på juli. Men det är en väldigt enkel modell, kanske man skulle behöva dubbelkolla och tänka vidare på.
linjär_modell <- tidy_data %>%
mutate(år = if_else(datum >= "2020-01-01", "2020", "2015-2019")) %>%
filter(month(datum) == 7) %>%
lm(antal_döda ~ år, data = .)
summary(linjär_modell)##
## Call:
## lm(formula = antal_döda ~ år, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.032 -12.032 -0.903 10.968 49.968
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 225.032 1.393 161.536 <2e-16 ***
## år2020 -5.129 3.412 -1.503 0.135
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.34 on 184 degrees of freedom
## Multiple R-squared: 0.01213, Adjusted R-squared: 0.006761
## F-statistic: 2.259 on 1 and 184 DF, p-value: 0.1345Vi har ett stort P-värde och även om estimatet är negativt, d.v.s. dödligheten är lägre under 2020, så är inte skillnaden signifikant.
Inspektera residualerna.
