Na wojnie z koronawirusem najlepszą bronią jest informacja

To był dla mnie szalony rok, stąd brak aktywności na blogu fundacji SmarterPoland w ostatnich miesiącach. Jednak dzisiejszy artykuł na zdrowie.gazeta.pl sprowokował mnie do zabrania głosu nt. skandalicznej (nie)dostępności danych o rozwoju pandemii CoviD. O tym poniżej, ale jeżeli interesuje Cię tytułowy obrazek to więcej informacji o nim znajdziesz na końcu.

Ale po kolei.

Od kilku miesięcy współpracuję z zespołem MOCOS.pl, międzynarodową grupą badaczy zasilaną nieskończoną energią lidera – prof. Tylla Krugera z Politechniki Wrocławskiej. Jakiś czas temu zespół MOCOS przedstawił zbiór 10 rekomendacji w celu powstrzymania dalszego rozwoju epidemii. Poza punktami takimi jak testuj, ograniczaj kontakty, itp, znalazła się na tej liście rekomendacja szczególnie mi bliska, czyli ”Strategia cyfrowego zbierania i publicznego udostępniania danych”. Tytuł tego posta pochodzi właśnie z wspomnianej rekomendacji.

Dlaczego dostępność dobrej jakości danych jest taka ważna?

Brak dobrej jakości danych paraliżuje prace niezliczonej ilości podmiotów. Począwszy od modelarzy, którzy miesiącami opracowują modele predykcyjne dla możliwych scenariuszy pandemii, przez właścicieli małych biznesów, którzy nie wiedząc jak wygląda sytuacja nie mogą efektywnie planować swojej pracy, po szeregowych obywateli, którzy po prostu chcą wiedzieć co się dzieje. Bardzo chcą wiedzieć, więc jeżeli nie będzie oficjalnych rzetelnych danych to tworzy się przestrzeń dla fakenewsów.
Podstawowym wymogiem jest pomiary były porównywalne. W gwarze statystycznej mówi się o unikaniu porównywania jabłek do pomarańczy.
Tymczasem w przypadku danych covidowych raportowanych w Polsce zmieniają się sposoby liczenia podstawowych statystyk, takich jak liczba testów, liczba zakażonych (patrz np. tutaj). Przez to raportowane 10 zdiagnozowanych chorych liczonych we wrześniu liczonych jest w inny sposób niż 10 chorych w listopadzie. Podobnie z wykonanymi testami (chodzi o sposób w jaki wliczane są testy wykonane prywatne). Jeżeli w danych nie mamy informacji jak były liczone poszczególne współczynniki w kolejnych dniach, to dane są słabej jakości.

Równie ważna co jakość danych jest ich dostępność. Z niepojętego dla mnie powodu Ministerstwo Zdrowia jako główny kanał komunikacji o liczbie zakażeń wybrało Twittera. Aby robić jakiekolwiek powtarzalne analizy dane muszą być dostępne! Twitter to nie jest format przechowywania danych ani nawet nie jest dobre medium do udostępniania danych. Przed wakacjami dr. hab. Anna Ochab-Marcinek napisała parser przetwarzający twity Ministerstwa w bardziej dostępny format umieszczony na GitHubie. Ale albo sposób raportowania się zmienił, albo API Twittera się zmieniło i ten sposób pozyskiwania danych przestał działać.
Obecnie chyba wszyscy polscy modelarze korzystają z bazy danych ręcznie uzupełnianej przez Michała Rogalskiego (w wolnym czasie). Baza jest dostępna jako google spreadsheet. Cieszy mnie proaktywna postawa osób gotowych poświęcić dużo własnego czasu aby coś pożytecznego zrobić, ale spodziewałbym się, że w XXI wieku średniej wielkości państwo jest w stanie wystawić w kilka miesięcy oficjalne API z dostępem do kluczowych danych.

Pomijając już sam format danych. Całkowicie niezrozumiały jest dla mnie fakt, że Ministerstwo Zdrowia przestało publikować informacje o wieku, płci i chorobach towarzyszących zmarłych osób. Zasłaniając się tym, że zgonów jest tak dużo, że nie zmieściłyby się w twitach na Twitterze (tak jakby nie można było tych danych umieścić na stronach www MZ).

Równie ważna co jakość i dostępność danych jest też odpowiednia prezentacja danych. Od początku pandemii zbieram z różnych mediów sposoby prezentowania danych covidowych. Nigdy wcześniej nie widziałem takiego zainteresowania w mediach wykresami i statystykami. Niestety często przy okazji dobrych chęci wychodzi też na jaw kompletny brak umiejętności analizy i prezentacji danych w sposób czytelny. Dotyczy to zarówno gazet jak i oficjalnych rządowych komunikatów (widać nie jest czytelny, skoro wspomniane zdrowie.gazeta poświęciło cały artykuł na próbę zrozumienia jak czytać prezentowane statystyk śmiertelności).
Przykładowo, w danych jest ewidentny efekt dnia tygodnia wynikający z tego, że przez weekend wykonywanych jest mniej testów. Zanim zacznie się liczyć i pokazywać jakikolwiek trend, absolutnie podstawowym krokiem przetwarzania wstępnego jest uwzględnienie tego efektu. Niewiele jest jednak gazet, które to robią (jedyna, którą znam, pokazuje skumulowaną liczbę przypadków z ostatnich 7 dni, co jest jakąś formą radzenia sobie z tym problemem).

Gdybyśmy mieli dobrej jakości danych, moglibyśmy przygotowywać najróżniejsze raporty i statystyki (tutaj jest kilka statystyk opracowanych w grupie MOCOS).

Bez dobrych danych jesteśmy (społeczeństwo, modelarze, lekarze i też pewnie rząd) ślepi.

Wspomniane rekomendacje grupy MOCOS są dostępne tutaj: https://mocos.pl/pl/recommendations.html.

Korzystając z danych udostępnionych grupie MOCOS z okresu maj – wrzesień przygotowaliśmy kalkulator opisujący ryzyka zgonu lub hospitalizacji. Kalkulator ten jest dostępny na stronie https://crs19.pl/. Dane są dosyć stare, więc trudno powiedzieć na ile te relacje utrzymają się w listopadzie. Ten kalkulator to jedna z wielu rzeczy, które można by zrobić gdyby mieć dobrej jakości dane.

RBioMeSs – R, uczenie maszynowe, statystyka medyczna i bioinformatyka

masterR

TL;DR: 24 listopada, w ramach Spotkań Entuzjastów R, odbędzie się spotkanie poświęcone R, bioinformatyce i statystyce medycznej. Więcej informacji tutaj.

LV:
Rozmawiałem ostatnio ze znajomym o ciekawych wyzwaniach związanych z analizą dużych danych. Zaczęło się od wyników w obszarze sieci konwolucyjnych i deep learningu ale zbaczaliśmy na różne tematy gdzie dane są niemałe a wyzwania być może i większe.
Gdy myśleć o klasyfikacji obrazów o rozmiarach 64×64 piksele (4096 piksle) to o ileż bardziej złożona jest predykcja losów pacjenta na bazie ekspresji dla 20 tysięcy genów czy informacji o stanie mutacji/metylacji dla milionów sond (miliony markerów dla każdego pacjenta! to już jest wysokowymiarowa przestrzeń).

Czytaj dalej RBioMeSs – R, uczenie maszynowe, statystyka medyczna i bioinformatyka

Czy przekroczą 55 milionów?

Już jutro finał 24. Wielkiej Orkiestry Świątecznej Pomocy. Z roku na rok WOŚP zbiera coraz więcej środków, w tym roku na wsparcie oddziałów pediatrycznych i opieki medycznej seniorów. Jak myślicie ile pieniędzy uda się zebrać?
Zobaczmy co na ten temat mają do powiedzenia modele liniowe 😉

Screen Shot 2016-01-09 at 02.34.59

Czytaj dalej Czy przekroczą 55 milionów?

Tabela 1 a pakiet Gmisc

Tworząc raporty często początkowe tabele są do siebie podobne – przedstawiają statystyki opisowe zmiennych. Bardziej złożone statystyki są zazwyczaj później.

W przypadku prac bio-medycznych używa się sformułowania Tabela 1 – czyli pierwsza tabela w artykule, zazwyczaj przedstawiająca statystyki opisowe porównywanych grup (np terapia A/B/C).

Ostatnio odkryłem pakiet Gmisc – fantastyczne wsparcie do szybkiego tworzenia (dobrze wyglądających) tabel prosto z poziomu knitra. Poniższy przykład dotyczy tabeli z podsumowaniem, ale możliwości tego pakietu są znacznie większe.

Przykładowo, taki krótki kawałek kodu (nazwy ważnych zmiennych usunąłem ze zrozumiałych powodów)

Czytaj dalej Tabela 1 a pakiet Gmisc

Myszy, testy post hoc i diffogramy


Pracowałem ostatnio z ciekawym problemem.

Mamy dwa rodzaje myszy. Z każdego rodzaju wybieramy trzy osobniki. Jesteśmy ciekawi jak pewien sposób traktowania komórek nerwowych wpływa na kolce dendrytyczne, czyli małe wypustki na neuronach. Z każdej myszy pobieramy ileś skrawków mózgu, każdy skrawek traktujemy na dwa badane sposoby. Każdemu skrawkowi robimy kilkadziesiąt zdjęć. Na każdym zdjęciu oznaczamy setki kolców dendrytycznych.

Mierzymy parametry kolców (np. długość), ale ich pomiary nie są niezależne. Musimy uwzględnić efekt zdjęcia, zagnieżdżony w efekcie skrawka, zagnieżdżony w efekcie myszy, zagnieżdżony w efekcie typu myszy. A interesuje nas porównanie pomiędzy sobą czterech warunków eksperymentalnych (skrzyżowane efekty typu myszy i sposobu traktowania).

Czym to modelować?

Czytaj dalej Myszy, testy post hoc i diffogramy

Jak często w Polsce występują depresje?

Temat depresji ostatnio często gości na łamach gazet, głównie za sprawą Justyny Kowalczyk. W różnych gazetach różni eksperci wskazują grupy narażone na wyższe ryzyko depresji. Skłoniło mnie to do poszperania za statystykami depresji.

Nie było łatwo/ Znalazłem statystyki z GUS z roku 2009 z raportu ,,Stan zdrowia ludności Polski w 2009 r.”. GUS ma też starsze raporty, ale tam depresja jest raportowana razem z nerwicami i są też raporty nowsze, ale te nie zawierają podziału na grupy wiekowe.

Graficzne podsumowanie poniżej.



Jak zaplanować płeć dziecka? Co wyjdzie z fuzji Big Data, advanced analytics i PISA 2012?

Istnieje wiele naturalnych metod planowania płci dziecka (oparta o kalendarz, pole magnetyczne, dietę itp), jednak większość z nich ma zerową (a w najlepszym przypadku kilkuprocentową) skuteczność.

Okazuje się jednak, że istnieją metody znacznie skuteczniejsze, choć niespodziewane. Udało się je odkryć dzięki analizom Big Data danych z międzynarodowego badania PISA 2012. Poniżej przedstawiamy wyniki dla Polski, ale podobne otrzymuje się dla praktycznie każdego kraju. To nie może być przypadek!

Czytaj dalej Jak zaplanować płeć dziecka? Co wyjdzie z fuzji Big Data, advanced analytics i PISA 2012?

Dostępność lekarzy specjalistów, tutorial

Dwa dni temu, w tym wpisie przedstawialiśmy grafiki przygotowane przez Michała Kurtysa w ramach wakacyjnego projektu.

Pan Michał przygotował też krótki tutorial wyjaśniający jak samodzielnie zrobić takie wykresy w R. Poniżej wklejamy ten tutorial. Jest on moim zdaniem bardzo ciekawy i porusza wiele ciekawych technicznych problemów z analizą danych przestrzennych.

Dane
Informacje o kolejkach zapisane są w plikach Excelowych o rozszerzeniu xls.
Na każde województwo przypada takich plików kilka. Nas teraz będzie interesował tylko plik dotyczący świadczeń specjalistycznych, zawierający w nazwie skrót AOS.
W sumie jest więc 16 plików (po jednym na województwo), które wyglądają mniej-więcej w ten sposób.
01_AOS_31072013.xls
02_AOS_31072013.xls
...
16_AOS_31072013.xls

Województwa są ułożone alfabetycznie – najmniejszy numer odpowiada województwu dolnośląskiemu, a największy zachodniopomorskiemu. Aby powiązać numer z nazwą województwa przygotowałem plik województwa.csv, w którym są one wypisane w należytej kolejności.

wojwództwa.csv
WOJ. DOLNOŚLĄSKIE
WOJ. KUJAWSKO-POMORSKIE
WOJ. LUBELSKIE
WOJ. LUBUSKIE
WOJ. ŁÓDZKIE
WOJ. MAŁOPOLSKIE
WOJ. MAZOWIECKIE
WOJ. OPOLSKIE
WOJ. PODKARPACKIE
WOJ. PODLASKIE
WOJ. POMORSKIE
WOJ. ŚLĄSKIE
WOJ. ŚWIĘTOKRZYSKIE
WOJ. WARMIŃSKO-MAZURSKIE
WOJ. WIELKOPOLSKIE
WOJ. ZACHODNIOPOMORSKIE

Dane o granicach administracyjnych Polski pobrałem z Geoportalu. Ze względów licencyjnych początkowo chciałem wykorzystać OpenStreetMap. Gotowe pliki shapefile oferuje m.in. geofabrik/cloudmade. Niestety, jeżeli chodzi o Polskę nie miałem tam czego szukać – wszystkie mapy były zdeformowane.
Mapy z Geoportalu nie są dostępne w formacie shapefile. Na szczęście dzięki instrukcjom we wpisie Pawła Wiechuckiego przedstawianego na tym blogu ich własnoręczne stworzenie nie było trudne.

Biblioteki
Potrzebujemy następujących bibliotek:
** maptools – funkcja readShapePoly pozwala na wczytanie pliku ShapeFile.
** sp – do manipulacji i wyświetlania danych kartograficznych.
** rgeos, rgdal – funkcja spTransform, pozwala zmienić układ współrzędnych
** FNN – funkcja get.knnx – do regresji
** SmarterPoland – funkcję getGoogleMapsAddress ułatwi pobranie współrzędnych punktu o określonym adresie
** Cairo – do produkcji wykresów
** XLConnect – obsługa plików excela

Kod R

Zacznijmy od wczytania bibliotek.

1
2
3
4
5
6
7
8
library(maptools)
library(sp)
library(rgeos)
library(SmarterPoland)
library(FNN)
library(rgdal)
library(Cairo)
library(XLConnect)

Następnie zmiana katalogu roboczego. Wczytanie listy województw.

Z każdego pliku „aos” odczytujemy dane, które zaczynają się w 3 wierszu(razem z nagłówkiem) i mają 9 kolumn.
Postawnowiłem zmienić nazwy kolumn, gdyż oryginalne były niesłychanie długie.
Dodatkowo tworzymy dodatkowe kolumny – ID i nazwę województwa.

9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
setwd("C:\\Users\\mic\\Documents\\eurostat\\nfzqueue\\new")
wojewodztwa = read.csv("wojewodztwa.csv", header=FALSE)
wojewodztwa$V1 = as.character(wojewodztwa$V1)
 
poradnie = data.frame()
for( i in 1:16) {
	filename = sprintf("%02d_AOS_31072013.xls", i)
	print(filename)
	wb = loadWorkbook( filename, create = FALSE)
	dane = readWorksheet(wb,sheet="Zestawienie", startRow = 3, startCol=0, endCol=8)
	dane$IDWojewodztwa = i
	dane$Wojewodztwo = wojewodztwa$V1[i]
	poradnie = rbind(poradnie, dane)
}
 
#zmiana nazwy kolumn
names(poradnie) = c(
"Nazwa.komorki.realizujacej",
"Kategoria",
"Nazwa",
"Nazwa.komorki",
"Adres",
"Liczba.oczekujacych",
"Liczba.skreslonych",
"Sredni.czas",
"ID.wojewodztwa",
"Wojewodztwo")

Pacjenci są podzieleni przez NFZ na dwie grupy – „przypadek stabilny” i „przypadek ostry”.
Graficznie przedstawiać będziemy wyłącznie dane o przypadkach stabilnych.

Dalej dzielimy adres na części składowe – nazwa miejscowości i ulica razem z numerem numer lokalu.

36
37
38
39
poradnie = poradnie[which(poradnie$Kategoria == "przypadek stabilny"), ]
poradnie$Miejscowosc = sapply(strsplit(poradnie$Adres, "\n"), function(x) x[1] )
poradnie$Ulica = sapply(strsplit(poradnie$Adres, "\n"), function(x) x[2] )
poradnie$Sredni.czas = as.numeric(poradnie$Sredni.czas)

Pobieramy współrzędne wyłącznie miejscowości w której znajduje się poradnia.
Jest to znacznie szybsze, a z pewnością nie potrzebujemy większej dokładności.
Po za tym Google ogranicza darmowy dostęp do usługi do 2500 zapytań dziennie.

Konstruujemy więc tabelę zawierającą nazwę miejscowości w jednej kolumnie, a w drugiej województwo w którym się ona znajduje.
Dodanie województwa pomaga uściślić zapytanie – istnieją przecież miejscowości, które noszą tę samą nazwę.

Funkcja unique eliminuje zduplikowane wiersze.
W pętli tworzymy nową zmienną region, w której zapisujemy jedynie nazwę województwa bez skrótu „Woj.” na początku.
Google Geocoding raczej nie trawi tego przedrostka.

Województwo łączymy z nazwą miejscowości i przekazujemy do argumentu city funkcji getGoogleMapAddress.
Może to wyglądać to dziwnie, ale zapytanie i tak zostanie skonstruowane poprawnie.

40
41
42
43
44
45
46
47
48
49
50
51
52
places_unique = data.frame( poradnie$Miejscowosc, poradnie$Wojewodztwo, check.names=FALSE)
places_unique = unique( places_unique )
names(places_unique) = c("Miejscowosc", "Wojewodztwo")
places_unique$Wojewodztwo = as.character(places_unique$Wojewodztwo)
places_unique$Miejscowosc = as.character(places_unique$Miejscowosc)
 
for( i in 1:nrow(places_unique) ) {
	region = strsplit( places_unique$Wojewodztwo[i], " " )[[1]][2]
	temp_coords = getGoogleMapsAddress(street="",city=paste(places_unique$Miejscowosc[i],",",region))
 
	places_unique$lat_wgs84[i] = temp_coords[1]
	places_unique$lng_wgs84[i] = temp_coords[2]
}

Zawczasu warto przygotować sobie współrzędne w układzie, który jest bardziej przyjazny do zadania:
Tworzymy zmienną typu SpatialPoints i określamy ich układ współrzędnych.
Współrzędne, które pobraliśmy są w układzie odniesienia WGS 84, znanym także jako EPSG:4326

Zmieniamy układ współrzędnych przy pomocy funkcji spTransform.

Ostatecznie łączymy dwie tabele przy pomocy funkcji merge.
Wreszcie mamy wszystkie interesujące nas wartości.

53
54
55
56
57
58
59
60
places_sp = SpatialPoints( cbind(places_unique$lng_wgs84, places_unique$lat_wgs84) )
proj4string(places_sp) = CRS("+init=epsg:4326")
places_sp = spTransform( places_sp, CRS("+init=epsg:2180"))
 
places_unique$lng = coordinates(places_sp)[,1]
places_unique$lat = coordinates(places_sp)[,2]
 
poradnie = merge(poradnie, places_unique)

Wczytanie pliku shapefile, który zawiera granice administracyjne województw.
Jest już w prawidłowym układzie współrzędnych – wystarczy tylko go określić.

61
62
wojewodztwa.shp = readShapePoly("geoportal\\woj.shp")
proj4string(wojewodztwa.shp)=CRS("+init=epsg:2180")

Funkcja generate_grid posłuży do wygenerowania kraty, która będzie pokrywać powierzchnię określonego pola.
Na obrazku widać punkty kraty. W istocie są to punkty, które będą służyły jako punkty startowe linii.

Pierwszy argument funkcji to ilość komórek kraty w osi X.
Komórka kraty ma być kwadratem, więc ich ilość w osi Y nie będzie przekazywana do funkcji.

Następnie wyliczamy topologię kraty. Przekazujemy:
** współrzędne dolnego, lewego rogu
** wymiary komórek
** ilość komórek
I na koniec zwracamy obiekt klasy SpatialGrid.

63
64
65
66
67
68
69
70
71
72
73
generate_grid = function( cell_n_x=1000, sp_polygons) {
	bbox_min = bbox(sp_polygons)[,1]
	bbox_max = bbox(sp_polygons)[,2]
	cell_width=((bbox_max-bbox_min)/cell_n_x)[1]
	cell_height=cell_width
	cell_n_y = ceiling(((bbox_max-bbox_min)/cell_height)[2])
 
	gt = GridTopology( bbox_min, c(cell_width,cell_height), c(cell_n_x,cell_n_y)  )
	sg = SpatialGrid(gt, CRS("+init=epsg:2180"))
	return(sg)
}

Jak widać na obrazku, część punktów kraty znajduje się poza granicami kraju.
Funkcja which_cells_inside zwraca indeksy komórek kraty, które znajdują się wewnątrz obiektu SpatialPolygons.

Wykorzystamy do tego funkcję over. Dla argumentów SpatialPoints i SpatialPolygons funkcja zwraca wektor o długości równej ilości punktów.
Wartości wektora mówią o tym, wewnątrz którego wielokąta znajduje się punkt.
Tych mamy 16 – odpowiadają województwom.

Jeżeli jest poza jakimkolwiek wielokątem składowym to danemu punktowi odpowiadać będzie wartość NA.

Szczegółowo wielokąt składowy to obiekt klasy Polygons:

74
75
76
77
78
79
80
81
82
83
84
85
> class(wojewodztwa.shp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> class(wojewodztwa.shp@polygons)
[1] "list"
> class(wojewodztwa.shp@polygons[[1]])
[1] "Polygons"
attr(,"package")
[1] "sp"
> length(wojewodztwa.shp@polygons)
[1] 16
86
87
88
89
90
91
92
93
which_cells_inside = function( sp_grid, sp_polygons) {
	grid_coords = coordinates(sp_grid)
	grid_points = SpatialPoints(grid_coords, CRS( proj4string(sp_grid)) )
	inside = over( grid_points, sp_polygons)
	proj4string(grid_points) = proj4string(wojewodztwa.shp)
	inside = which( inside$ID1 > 0 )
	return(inside)
}

Generowanie linii. W argumentach znajduje się zmienna cell_distance, która mówi co ile komórek kraty będzie znajdować się początek linii.
Gdy pisałem kod uznałem, że być może takie rozwiązanie się przyda. Teraz jestem pewien, że bardziej elegancko byłoby generować mniejszą kratę.
Co z resztą każdy zauważy.

Zwracamy obiekt typu SpatialLines. Strukura obiektów (Spatial)Line(s) jest podobna do Polygons.

Przypominam, że linie rysujemy od punktu kraty do najbliższej przychodni.
Musimy więc wiedzieć, która jest najbliższa.
Przekazujemy więc macierz knn_indices w której będzie to zapisane.

94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# sp_points to obiekt typu SpatialPoints, zawierający listę przychodni.
# cell_subset ogranicza kratę, tak by linie były prowadzone tylko z wewnątrz kraju. 
 
generate_lines = function( knn_indices, sp_grid, sp_points, cell_distance=1, cell_subset=NULL ) {
	cell_n_x = sg@grid@cells.dim[1]
	cell_n_y = sg@grid@cells.dim[2]
	#macierz ktora zawiera indeksy komorek w ktorych zaczynaja sie linie
	line_anchors = matrix(nrow=floor(cell_n_y/cell_distance), ncol=floor(cell_n_x/cell_distance))
	for( i in 1:nrow(line_anchors)) {
		line_anchors[i,] = seq(cell_distance, cell_n_x, by= cell_distance) +
						  rep(cell_n_x*cell_distance*i,floor(cell_n_x/cell_distance))
	}
	#tylko te linie ktorych poczatek jest w wybranych komorkach
	#na przyklad wewnatrz kraju
	line_anchors = intersect(line_anchors, cell_subset)
 
	#konstrukcja obiektow biblioteki sp
	line_starts = coordinates(sp_grid)[line_anchors,]
	closest_medic_idx = knn_indices[line_anchors,1]
 
	line_stops = coordinates(sp_points)[closest_medic_idx,]
	line_list = sapply( 1:nrow(line_starts), function(x) Line( rbind( line_starts[x,], line_stops[x,])) )
	line_obj = Lines(line_list, ID="a")
	sp_lines_obj = SpatialLines( list(line_obj) )
 
	return(sp_lines_obj)
}

Odpalamy funkcje.
W pętli przejdziemy po kolei po ka¿dym typie poradni.
Wewnątrz niej wybierzemy interesujące nas wiersze w tabeli poradnie.
Dzięki bibliotece Cairo łatwo zapiszemy wykres w wysokiej rozdzielczości.
get.knnx jest funkcją obsługującą regresję k-sąsiadów.
W pierwszym argumencie jest zbiór danych, w drugim zbiór zapytań.
W tej sytuacji interesuje nas wyłącznie najbliższy sąsiad, więc k wynosi 1.
Zwraca dwie macierze – w pierwszej zapisuje indeksy najbliższych sąsiadów, a w drugiej odległości.

123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
#wczytanie danych, etc..
#kod
sg = generate_grid( cell_n_x=100, wojewodztwa.shp)
inside_poland = which_cells_inside( sg, wojewodztwa.shp)
 
for( specialist_type in unique( (poradnie$Nazwa.komorki.realizujacej) ) ) {
	print(specialist_type)
	flush.console()
	specialist = poradnie[which(poradnie$Nazwa.komorki.realizujacej == specialist_type),]
	specialist_sp = SpatialPoints( cbind( specialist$lng, specialist$lat) ,CRS("+init=epsg:2180") )
 
	knn_res = get.knnx( specialist_sp@coords, coordinates(sg), k = 1 )
 
	sp_lines_obj = generate_lines( knn_res[[1]], sg, specialist_sp, cell_distance=1, cell_subset=inside_poland )
 
	Cairo(900, 700, file=paste(specialist_type,".png",sep=""), type="png", bg="white")
	plot(wojewodztwa.shp)
	plot(sp_lines_obj,add=T)
	#plot(specialist_sp,col="blue",lwd=2, add=T)
	title(main=specialist_type)
	dev.off()
}

Hidden track. Tworzenie map „kolorowych”. Tworzymy gęstszą kratę.
Następnie uzyskujemy punkty kraty i zmieniamy ich układ współrzędnych do WGS 84.

Jest to niezbędne – dla układu współrzędnych epsg:2180 funkcja spDistsN1 nie zwracała odległości w kilometrach.

145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
sg = generate_grid( cell_n_x=500, wojewodztwa.shp)
inside_poland = which_cells_inside( sg, wojewodztwa.shp)
 
grid_points_wgs84 = SpatialPoints( coordinates(sg), CRS(proj4string(sg)))
grid_points_wgs84 = spTransform(grid_points_wgs84, CRS("+init=epsg:4326"))
gp_coords = coordinates(grid_points_wgs84)
 
for( specialist_type in unique( (poradnie$Nazwa.komorki.realizujacej) ) ) {
	print(specialist_type)
	flush.console()
	specialist = poradnie[which(poradnie$Nazwa.komorki.realizujacej == specialist_type),]
	specialist_sp = SpatialPoints( cbind( specialist$lng, specialist$lat) ,CRS("+init=epsg:2180") )
	specialist_coords_wgs84 = cbind( specialist$lng_wgs84, specialist$lat_wgs84)
 
	#knn_res = get.knnx( specialist_sp@coords, coordinates(sg), k = ifelse(nrow(specialist)>10, 10, nrow(specialist) ))
	knn_res = get.knnx( specialist_sp@coords, coordinates(sg), k = 1)
 
	knn_indices = knn_res[[1]]
	km_distances =  matrix( nrow=nrow( knn_indices ), ncol=ncol(knn_indices))
	closest_coords = apply( knn_indices, c(1,2), function(x) specialist_coords_wgs84[x,])
	closest_coords = matrix( closest_coords, ncol = 2, byrow=TRUE ) #2*ncol(knn_indices)
 
	for (rowID in 1:nrow( knn_indices ) ) {
		pts = matrix(closest_coords[rowID,],ncol=2)
 
		p = gp_coords[rowID,]
		km_distances[rowID,] = spDistsN1(pts, p, TRUE)
	}
	closest_distance = km_distances[,1]
	sg_df = SpatialGridDataFrame( sg, data.frame(closest_distance))
 
	#obszar położony poza granicami kraju ma kolor biały
	sg_df@data$closest_distance[-inside_poland]=0
	#80 jest górną granicą wartości na wykresie
	sg_df@data$closest_distance[ sg_df@data$closest_distance>80] = 80
 
	Cairo(900, 700, file=paste(specialist_type,"_c",".png",sep=""), type="png", bg="white")
	print(spplot( sg_df, panel = function(x,y,...){
	panel.gridplot(x,y,...,  )
	sp.polygons( wojewodztwa.shp )
	sp.points( specialist_sp, pch=".", col="black" )
	}, at=c(1:80), col.regions = rev(heat.colors(100)), main=specialist_type))
	dev.off()
 
}

Było to mój pierwszy ,,większy” projekt w R. Z pewnością wiele rozwiązań nie jest najlepszych, ale mimo wszystko mam nadzieję, że ten wpis będzie pomocny.

Dostępność lekarzy specjalistów

Ten wpis jest przygotowany przez Michała Kurtysa podczas wakacyjnego projektu wykonanego w ramach wolontariatu dla naszej fundacji.

Pustynie dostępności do lekarzy specjalistów, czyli gdzie w Polsce najdalej trzeba jechać do specjalisty.
Michał Kurtys

Zainteresowany wakacyjnym wolontariatem dla Fundacji SmarterPoland zgłosiłem się do Pana Przemysława Biecka z pytaniem o interesujący temat.
Wziąwszy pod uwagę moje zainteresowania, interesującym tematem okazało się wizualizacja danych o czasach czekania w kolejkach do specjalisty. Dane te są publikowane przez NFZ i na tym blogu pojawiły się już dwa wpisy na ten temat: http://smarterpoland.pl/index.php/2013/07/jak-dlugo-czeka-sie-w-kolejce-do-lekarza/, http://smarterpoland.pl/index.php/2013/08/gdzie-w-polsce-czeka-sie-najdluzej-w-kolejce-do-lekarza-specjalisty/

Bardziej niż czas czekania do specjalisty zainteresował mnie problem odległości, w jakiej dostępny jest pożądany specjalista.
Jak daleko musi pojechać pacjent, żeby dostać się do specjalisty, do którego nie musi czekać w wielomiesięcznej kolejce.
Lub podobne pytanie, jak wielu ,,dostępnych” specjalistów znajduje się w umiarkowanej odległości od pacjenta.

Czytaj dalej Dostępność lekarzy specjalistów

Statystyki dotyczące nowotworu piersi

Jakiś czas temu agencja ExpertPR poprosiła naszą Fundację o zaprezentowanie statystyk dotyczących nowotworu piersi w Polsce. Temat bardzo mnie zainteresował (spędziłem w przeszłości kilka lat na współpracy z onkologami w analizie danych dotyczących nowotworów). Przygotowaliśmy więc prezentację, a Pani Prezes ją przedstawiła.

Dziś chciałbym w skrócie tę prezentację zaprezentować na blogu.

Link do prezentacji znajduje się tutaj. Poniżej prezentowane wykresy są ,,klikalne” i przekierowują na wskazaną stronę prezentacji. Cała prezentacja została wykonana w czystym R [o technikaliach napiszę więcej w czwartek].

W prezentacji przedstawione są wybrane informacje o epidemiologii choroby, o postępach w walce z tą chorobą i o postrzeganiu tej choroby [wybrane, temat jest zdecydowanie szerszy niż pozwalała na jego omówienie półgodzinna prezentacja]. Poniżej przedstawiam cztery wybrane slajdy.

Rośnie liczba zdiagnozowanych przypadków raka piersi

Jak widzimy na poniższym wykresie liczba zdiagnozowanych przypadków szybko rośnie. Ta informacja ma pozytywny wydźwięk. Dlaczego? Rosnąca liczba zdiagnozowanych przypadków jest związana z coraz wcześniejszym wykrywaniem nowotworu. A to z kolei znacząco zwiększa szansę na powodzenie leczenia.

Rośnie liczba artykułów na temat raka piersi

I artykułów naukowych (jest to bardzo gorąca dziedzina badań) i artykułów ,,popularnonaukowych”. To również pozytywne zjawisko. Jednym z problemów pacjentów jest trudność w dotarciu do informacji o chorobie, niepewność, jak leczenie wygląda w Polsce. Większa dostępność informacji o chorobie pomaga chorym przynajmniej w tym aspekcie.

Media społecznościowe a rak piersi

Okazuje się, że media społecznościowe spełniają ważną rolę grupy wsparcia dla osób chorych. Nie mówią tu o akcjach typu klikamy ,lubię tu', ale o miejscu gdzie osoby chore mogą się spotkać i podzielić obawami/przeżyciami/pytaniami.

Zaskoczony byłem, że np. na Twitterze funkcjonuje kanał #BCSM (Breast Cancer Social Media), na którym organizowane są wideo konferencje, odczyty, przemowy zaproszonych osób. W trakcie takiej konferencji na wskazanym kanale słuchacze dzielą się komentarzami.

Podobne inicjatywy można spotkać na innych mediach.

Przeżywalność 5-letnia a różne czynniki

Nikt nie wie jak potoczą się losy konkretnego chorego. Ale patrząc na średnie można odczytywać jak szanse na przeżycie np. pięciu lub dziesięciu lat zmieniają się w zależności od wielkości i złośliwości odkrytego guza. Zupełnie zgodne z oczekiwaniami jest to, że guzy wcześnie wykryte znacznie lepiej rokują, jeżeli chodzi o szanse leczenia. Ale o ile lepiej i jakie są ,,wyjściowe” szanse?

Więcej informacji można znaleźć w cytowanej prezentacji.