Metody pracy z danymi na przykładzie… COVID-19

Ten wpis ma na celu pokazać jak wczytać dane do R z pliku csv, a następnie jak nieco zmodyfikować postać tabeli z danymi i jak utworzyć prosty wykres. Ponieważ mamy obecnie epidemię wirusa SARS-CoV-2, dlatego (wiem, że dla niektórych osób nie jest to wdzięczny temat szkoleniowy), posłużymy się danymi związanymi z liczbą zachorowań na chorobę COVID-19, którą ten koronawirus wywołuje.

Na początek musimy pobrać aktualne (dane użyte w tym wpisie są z dnia 21.03.2020) „surowe” dane statystyczne w postaci pliku csv ze źródła zlokalizowanego pod adresem: https://data.humdata.org/dataset/novel-coronavirus-2019-ncov-cases. Pobrany plik o nazwie time_series_2019-ncov-Confirmed.csv najlepiej zapisać w katalogu domyślnym programu RStudio, wówczas łatwo go znajdziemy i wczytamy.

Teraz uruchamiamy RStudio i w edytorze wpisujemy kod jak poniżej:

packages <- c("utils",
              "dplyr",
              "tidyr",
              "ggplot2"
)
install_packages <- function(n = length(packages)) {
  for (i in 1:n) {
    if(packages[i] %in% rownames(installed.packages()) == FALSE) {
      install.packages(packages[i])
      library(packages[i], character.only = TRUE)
    } else {
      library(packages[i], character.only = TRUE)
    }
  }
}
install_packages()

setwd("~/R")
covid19_data <- read.csv2(file = "time_series_2019-ncov-Confirmed.csv", sep = ",")
covid19_data <- as.data.frame(covid19_data)
head(covid19_data, n = 10)

Linijki 1-16 to prezentowany już wcześniej kod instalowania i ładowania potrzebnych pakietów (zobacz wpis). Funkcja z linii 18 to polecenie ustawienia domyślnego katalogu R jako głównego katalogu zlokalizowanego w systemie Windows w folderze Dokumenty. Nasz pobrany plik musimy wczytać do pamięci programu i utworzyć obiekt o nazwie – dajmy na to, covid19_data. W tym celu użyjemy funkcji read.csv2() z pakietu utils, jak pokazano w kodzie w linii 19. W funkcji tej podajemy nazwę pliku wraz z rozszerzeniem oraz deklarujemy, jaki znak używany jest jako separator w naszym pliku csv, co jest potrzebne do prawidłowego rozdzielenia kolumn z danymi, w tym przypadku jest to przecinek. Przekształcamy nasz obiekt, który jest w tym momencie listą, na ramkę danych -linia 20. Po wykonaniu kodu wraz z kodem z linii 21, dzięki której podejrzymy zawartość wczytanego pliku, powinniśmy ujrzeć w konsoli programu RStudio poniższą zawartość:

> head(covid19_data, n = 10)
     Province.State Country.Region      Lat      Long X1.22.20 X1.23.20 X1.24.20 X1.25.20 X1.26.20 X1.27.20 X1.28.20 X1.29.20 X1.30.20 X1.31.20 X2.1.20
1                         Thailand       15       101        2        3        5        7        8        8       14       14       14       19      19
2                            Japan       36       138        2        1        2        2        4        4        7        7       11       15      20
3                        Singapore   1.2833  103.8333        0        1        3        3        4        5        7        7       10       13      16
4                            Nepal  28.1667     84.25        0        0        0        1        1        1        1        1        1        1       1
5                         Malaysia      2.5     112.5        0        0        0        3        4        4        4        7        8        8       8
6  British Columbia         Canada  49.2827 -123.1207        0        0        0        0        0        0        1        1        1        1       1
7   New South Wales      Australia -33.8688  151.2093        0        0        0        0        3        4        4        4        4        4       4
8          Victoria      Australia -37.8136  144.9631        0        0        0        0        1        1        1        1        2        3       4
9        Queensland      Australia -28.0167     153.4        0        0        0        0        0        0        0        1        3        2       3
10                        Cambodia    11.55  104.9167        0        0        0        0        0        1        1        1        1        1       1
   X2.2.20 X2.3.20 X2.4.20 X2.5.20 X2.6.20 X2.7.20 X2.8.20 X2.9.20 X2.10.20 X2.11.20 X2.12.20 X2.13.20 X2.14.20 X2.15.20 X2.16.20 X2.17.20 X2.18.20
1       19      19      25      25      25      25      32      32       32       33       33       33       33       33       34       35       35
2       20      20      22      22      45      25      25      26       26       26       28       28       29       43       59       66       74
3       18      18      24      28      28      30      33      40       45       47       50       58       67       72       75       77       81
4        1       1       1       1       1       1       1       1        1        1        1        1        1        1        1        1        1
5        8       8      10      12      12      12      16      16       18       18       18       19       19       22       22       22       22
6        1       1       1       2       2       4       4       4        4        4        4        4        4        4        4        5        5
7        4       4       4       4       4       4       4       4        4        4        4        4        4        4        4        4        4
8        4       4       4       4       4       4       4       4        4        4        4        4        4        4        4        4        4
9        2       2       3       3       4       5       5       5        5        5        5        5        5        5        5        5        5
10       1       1       1       1       1       1       1       1        1        1        1        1        1        1        1        1        1
   X2.19.20 X2.20.20 X2.21.20 X2.22.20 X2.23.20 X2.24.20 X2.25.20 X2.26.20 X2.27.20 X2.28.20 X2.29.20 X3.1.20 X3.2.20 X3.3.20 X3.4.20 X3.5.20 X3.6.20
1        35       35       35       35       35       35       37       40       40       41       42      42      43      43      43      47      48
2        84       94      105      122      147      159      170      189      214      228      241     256     274     293     331     360     420
3        84       84       85       85       89       89       91       93       93       93      102     106     108     110     110     117     130
4         1        1        1        1        1        1        1        1        1        1        1       1       1       1       1       1       1
5        22       22       22       22       22       22       22       22       23       23       25      29      29      36      50      50      83
6         5        5        6        6        6        6        7        7        7        7        8       8       8       9      12      13      21
7         4        4        4        4        4        4        4        4        4        4        4       6       6      13      22      22      26
8         4        4        4        4        4        4        4        4        4        4        7       7       9       9      10      10      10
9         5        5        5        5        5        5        5        5        5        5        9       9       9      11      11      13      13
10        1        1        1        1        1        1        1        1        1        1        1       1       1       1       1       1       1
   X3.7.20 X3.8.20 X3.9.20 X3.10.20 X3.11.20 X3.12.20 X3.13.20 X3.14.20 X3.15.20 X3.16.20 X3.17.20 X3.18.20 X3.19.20 X3.20.20
1       50      50      50       53       59       70       75       82      114      147      177      212      272      322
2      461     502     511      581      639      639      701      773      839      825      878      889      924      963
3      138     150     150      160      178      178      200      212      226      243      266      313      345      385
4        1       1       1        1        1        1        1        1        1        1        1        1        1        1
5       93      99     117      129      149      149      197      238      428      566      673      790      900     1030
6       21      27      32       32       39       46       64       64       73      103      103      186      231      271
7       28      38      48       55       65       65       92      112      134      171      210      267      307      353
8       11      11      15       18       21       21       36       49       57       71       94      121      121      121
9       13      15      15       18       20       20       35       46       61       68       78       94      144      184
10       1       2       2        2        3        3        5        7        7        7       33       35       37       51

Funkcja head(), z określonym parametrem n, pozwala na wyświetlenie tylko n pierwszych wierszy wczytanej tabeli. Bez wpisywania tego parametru, domyślnie wyświetlanych jest 6 pierwszych wierszy.

Aby zobaczyć listę z nazwami kolumn tabeli używamy polecenia:

colnames(covid19_data)

Otrzymujemy:

> colnames(cov19_data)
 [1] "Province.State" "Country.Region" "Lat"            "Long"           "X1.22.20"       "X1.23.20"       "X1.24.20"       "X1.25.20"      
 [9] "X1.26.20"       "X1.27.20"       "X1.28.20"       "X1.29.20"       "X1.30.20"       "X1.31.20"       "X2.1.20"        "X2.2.20"       
[17] "X2.3.20"        "X2.4.20"        "X2.5.20"        "X2.6.20"        "X2.7.20"        "X2.8.20"        "X2.9.20"        "X2.10.20"      
[25] "X2.11.20"       "X2.12.20"       "X2.13.20"       "X2.14.20"       "X2.15.20"       "X2.16.20"       "X2.17.20"       "X2.18.20"      
[33] "X2.19.20"       "X2.20.20"       "X2.21.20"       "X2.22.20"       "X2.23.20"       "X2.24.20"       "X2.25.20"       "X2.26.20"      
[41] "X2.27.20"       "X2.28.20"       "X2.29.20"       "X3.1.20"        "X3.2.20"        "X3.3.20"        "X3.4.20"        "X3.5.20"       
[49] "X3.6.20"        "X3.7.20"        "X3.8.20"        "X3.9.20"        "X3.10.20"       "X3.11.20"       "X3.12.20"       "X3.13.20"      
[57] "X3.14.20"       "X3.15.20"       "X3.16.20"       "X3.17.20"       "X3.18.20"       "X3.19.20"       "X3.20.20" 

Jak widać, w tabeli z liczbą potwierdzonych przypadków zachorowań, mamy kolumny o nazwie prowincji/stanu, kraju/regionu, z danymi o szerokości i długości geograficznej oraz kolumny przedstawiające kolejne daty począwszy od 22.01.2020, skończywszy na 20.03.2020.

Pakiet dplyr jest jednym z najczęściej używanych pakietów R, służy do manipulowania i organizowania danych w tabelach – do filtrowania, agregacji, podsumowań w nowych kolumnach, itp. Pełen opis tego pakietu znajdziecie pod tym linkiem: https://cran.r-project.org/web/packages/dplyr/dplyr.pdf. Spróbujmy teraz zastosować pakiet o nazwie dplyr do filtrowania danych, a konkretnie do sprawdzenia przypadków zachorowań na terenie Polski. Dodatkowo wybierzemy do wyświetlenia tylko przydatne kolumny:

covid19_data1 <- covid19_data %>%
  dplyr::filter(Country.Region == "Poland") %>%
  dplyr::select(-Province.State, -Lat, -Long)

poleceniem z linii 2 filtrujemy dane związane z przypadkami w naszym kraju, natomiast poleceniem z linii 3 usuwamy z wyświetlanej tabeli kolumny o podanych nazwach, używamy w tym przypadku znaku minus przed nazwą każdej kolumny. Dane zapisujemy w nowym obiekcie covid19_data1. Jak zauważyliście, w powyższym kodzie użyty został operator %>%, który służy do przenoszenia danych do kolejnej wykonywanej funkcji – jest to tzw. przetwarzanie potokowe. Tabela, jaką w rezultacie otrzymujemy, wygląda jak poniżej:

> covid19_data1
  Country.Region X1.22.20 X1.23.20 X1.24.20 X1.25.20 X1.26.20 X1.27.20 X1.28.20 X1.29.20 X1.30.20 X1.31.20 X2.1.20 X2.2.20 X2.3.20 X2.4.20 X2.5.20 X2.6.20
1         Poland        0        0        0        0        0        0        0        0        0        0       0       0       0       0       0       0
  X2.7.20 X2.8.20 X2.9.20 X2.10.20 X2.11.20 X2.12.20 X2.13.20 X2.14.20 X2.15.20 X2.16.20 X2.17.20 X2.18.20 X2.19.20 X2.20.20 X2.21.20 X2.22.20 X2.23.20
1       0       0       0        0        0        0        0        0        0        0        0        0        0        0        0        0        0
  X2.24.20 X2.25.20 X2.26.20 X2.27.20 X2.28.20 X2.29.20 X3.1.20 X3.2.20 X3.3.20 X3.4.20 X3.5.20 X3.6.20 X3.7.20 X3.8.20 X3.9.20 X3.10.20 X3.11.20 X3.12.20
1        0        0        0        0        0        0       0       0       0       1       1       5       5      11      16       22       31       49
  X3.13.20 X3.14.20 X3.15.20 X3.16.20 X3.17.20 X3.18.20 X3.19.20 X3.20.20
1       68      103      119      177      238      251      355      425

Dane ułożone są w sposób niezbyt wygodny i czytelny, ponieważ chcielibyśmy, aby daty były wypisane w kolejnych wierszach tabeli, a liczba przypadków podana była w oddzielnej kolumnie. Taki jak obecnie układ tabeli nazywany jest szerokim i bywa rzadziej stosowany niż, interesujący nas, układ wąski. Możemy w R przekształcić postać tabeli na układ wąski, w tym celu posłużymy się funkcją gather() z pakietu tidyr, jak poniżej w linii 1 oraz od razu usunąć z daty symbol X (linie 4-6):

covid19_data2 <- tidyr::gather(covid19_data1, "Country.Region")
colnames(covid19_data2) <- c("data", "przypadki")

for (i in 1:nrow(covid19_data2)) {
  covid19_data2[i, "data"] <- substr(covid19_data2[i, "data"], 2, 10)
}

Od razu warto zmienić nazwy kolumn, by były czytelniejsze – linia 2, mamy wówczas:

> covid19_data2
      data przypadki
1  1.22.20         0
2  1.23.20         0
3  1.24.20         0
4  1.25.20         0
5  1.26.20         0
6  1.27.20         0
7  1.28.20         0
8  1.29.20         0
9  1.30.20         0
10 1.31.20         0
11  2.1.20         0
12  2.2.20         0
13  2.3.20         0
14  2.4.20         0
15  2.5.20         0
16  2.6.20         0
17  2.7.20         0
18  2.8.20         0
19  2.9.20         0
20 2.10.20         0
21 2.11.20         0
22 2.12.20         0
23 2.13.20         0
24 2.14.20         0
25 2.15.20         0
26 2.16.20         0
27 2.17.20         0
28 2.18.20         0
29 2.19.20         0
30 2.20.20         0
31 2.21.20         0
32 2.22.20         0
33 2.23.20         0
34 2.24.20         0
35 2.25.20         0
36 2.26.20         0
37 2.27.20         0
38 2.28.20         0
39 2.29.20         0
40  3.1.20         0
41  3.2.20         0
42  3.3.20         0
43  3.4.20         1
44  3.5.20         1
45  3.6.20         5
46  3.7.20         5
47  3.8.20        11
48  3.9.20        16
49 3.10.20        22
50 3.11.20        31
51 3.12.20        49
52 3.13.20        68
53 3.14.20       103
54 3.15.20       119
55 3.16.20       177
56 3.17.20       238
57 3.18.20       251
58 3.19.20       355
59 3.20.20       425

Jak wiemy, pierwszy potwierdzony przypadek zachorowania na COVID-19 miał miejsce w dniu 4 marca (linia 45), zatem możemy zmniejszyć liczbę wierszy, tak by pierwszy wiersz wypadał 3 marca – począwszy od liczby przypadków równej zero – a przy tym jednocześnie poprawmy daty, by miały odpowiedni format.

covid19_data3 <- covid19_data2[42:59, ]
covid19_data3$data <- as.Date(paste0(3:20,"-3-2020"), format = "%d-%m-%Y")
rownames(covid19_data3) <- c(1:nrow(covid19_data3))
covid19_data <- cbind(covid19_data3, data.frame(numer_porzadkowy = c(1:nrow(covid19_data3))))

W linii 1 przypisujemy wiersze o numerach od 42 do 45 do nowego obiektu. W kolejnej linii do kolumny zawierającej daty wprowadzamy te same daty, ale tym razem w formacie %d-%m-%Y. Do utworzenia wektora dat użyta została funkcja paste0(), gdzie wyrażenie 3:20 pozwala na wpisanie ciągu dni. W linii 3 wpisujemy kolejne liczby jako nazwy wierszy. Efekt końcowy naszych działań zapisujemy do ramki danych covid19_data, którą tworzymy przez połączenie dotychczasowej tabeli z ramką danych zawierającą numery porządkowe. Efekt poniżej:

> covid19_data
         data przypadki numer_porzadkowy
1  2020-03-03         0                1
2  2020-03-04         1                2
3  2020-03-05         1                3
4  2020-03-06         5                4
5  2020-03-07         5                5
6  2020-03-08        11                6
7  2020-03-09        16                7
8  2020-03-10        22                8
9  2020-03-11        31                9
10 2020-03-12        49               10
11 2020-03-13        68               11
12 2020-03-14       103               12
13 2020-03-15       119               13
14 2020-03-16       177               14
15 2020-03-17       238               15
16 2020-03-18       251               16
17 2020-03-19       355               17
18 2020-03-20       425               18

Spróbujmy narysować prosty wykres pokazujący zmianę liczby potwierdzonych przypadków zachorowań w naszym kraju w kolejnych dniach marca. Wykorzystajmy funkcję ggplot() z pakietu ggplot2.

ggplot(covid19_data, aes(x = data, y = przypadki)) + 
  geom_point() + 
  scale_x_date(limits = as.Date(c("2020-03-03","2020-03-20")), 
               labels = date_format("%d"),
               breaks = date_breaks("days"))
Rys 1 Wzrost liczby zachorowań na COVID-19 w Polsce wg kolejnych dni marca

Widzimy jaki jest wzrost liczby przypadków, ale spróbujmy teraz pójść krok dalej i wyznaczmy prognozę przyrostu przypadków tej choroby w Polsce. Oczywiście nasze założenia co do prognozy opierać się będą wyłącznie na przypadkach historycznych i ich trendzie. Model wykładniczy, jaki zastosujemy, nie będzie uwzględniał żadnych działań mających na celu przeciwdziałać rozprzestrzenianiu się wirusa. Można powiedzieć, że model będzie naiwny, jak czasem określa się właśnie bardzo proste modele. Ale nie przejmujmy się tym, ponieważ chodzi głównie o demonstrację metody.

Dopasujmy funkcję wykładniczą postaci f(x)=ae^{{x}/{b}} do naszych punktów, która zwykle opisuje rozprzestrzenianie się epidemii. Do wykonania regresji nieliniowej stosujemy polecenie nls(), w której podajemy jaką nasza funkcja powinna mieć postać:

covid19_data_fit <- nls(formula = przypadki ~ a*exp(numer_porzadkowy/b), start = list(a = 0.1, b = 1), data = covid19_data)
summary(covid19_data_fit)
predict(covid19_data_fit)

Funkcją summary() wyświetlamy szczegóły o dopasowaniu funkcji modelowej do punktów doświadczalnych, natomiast polecenie predict() pozwala na wyznaczenie wartości funkcji modelowej dla założonych argumentów:

> summary(covid19_data_fit)

Formula: przypadki ~ a * exp(numer_porzadkowy/b)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   4.5259     0.7697    5.88 2.33e-05 ***
b   3.9302     0.1582   24.85 3.29e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.85 on 16 degrees of freedom

Number of iterations to convergence: 29 
Achieved convergence tolerance: 5.987e-06

> predict(covid19_data_fit)
 [1]   5.837233   7.528495   9.709778  12.523059  16.151452  20.831125  26.866672  34.650940  44.690597  57.639115  74.339296  95.878138 123.657577 159.485747 205.694662 265.292007 342.156905 441.292405

Zobaczmy, efekt dopasowania również na wykresie:

ggplot(covid19_data2, aes(x = numer_porzadkowy, y = przypadki)) + 
  geom_point() +
  geom_smooth(method = "loess", formula = y ~ exp(x/3.9302))
Rys 2 Dopasowanie funkcji modelowej do puntów doświadczalnych

Teraz porównajmy liczbę przypadków zachorowań w kolejnych dniach z liczbą przypadków wyznaczoną w oparciu o model oraz policzmy współczynnik R^{2} dla danych modelowych i historycznych:

fitting_points <- data.frame(data = covid19_data$data, przypadki = covid19_data$przypadki, prognoza_przypadki = predict(covid19_data_fit))
R2 <- cor(fitting_points$prognoza_przypadki, fitting_points$przypadki)^2
> fitting_points
         data przypadki prognoza_przypadki
1  2020-03-03         0           5.837233
2  2020-03-04         1           7.528495
3  2020-03-05         1           9.709778
4  2020-03-06         5          12.523059
5  2020-03-07         5          16.151452
6  2020-03-08        11          20.831125
7  2020-03-09        16          26.866672
8  2020-03-10        22          34.650940
9  2020-03-11        31          44.690597
10 2020-03-12        49          57.639115
11 2020-03-13        68          74.339296
12 2020-03-14       103          95.878138
13 2020-03-15       119         123.657577
14 2020-03-16       177         159.485747
15 2020-03-17       238         205.694662
16 2020-03-18       251         265.292007
17 2020-03-19       355         342.156905
18 2020-03-20       425         441.292405

> R2
[1] 0.9911105

Jak widzimy, model jest bardzo dobrze dopasowany do punktów doświadczalnych. Policzmy w takim razie, ile wg modelu powinno być potwierdzonych zachorowań w dniu 31 marca. Dzień ten będzie miał przyporządkowany numer 29. w nomenklaturze naszego numeru porządkowego.

a <- summary(covid19_data_fit)$coefficient[1]
b <- summary(covid19_data_fit)$coefficient[2]

a*exp(29/b)
> a*exp(29/b)
[1] 7248.348

Wg tego bardzo prostego modelu na dzień 31 marca możemy spodziewać się liczby zachorowań na COVID-19 na poziomie około 7248 przypadków.

W tym wpisie chciałem nakreślić sposoby pracy z danymi i jest to przegląd metod, jednak może być dobrym wstępem do dalszego ich zgłębiania. Oczywiście mam nadzieję, że przedstawiona tutaj prognoza nie sprawdzi się w naszym kraju, a tempo zachorowań na COVID-19 będzie wyhamowywać.

Dodaj komentarz

Twój adres email nie zostanie opublikowany. Pola, których wypełnienie jest wymagane, są oznaczone symbolem *