-Door dr. Erik Giltay, psychiater en epidemioloog in het LUMC-
Dit is een blog voor studenten, co-assistenten, arts-assistenten en anderen met interesse voor onderzoek in de psychiatrie. Ik richt me deels op vraagstellingen die met open source databestanden benaderd kunnen worden. Ik zal proberen om de verschillende stappen van het ‘data processing’ proces te laten zien, gebruik makend van R code.
Data opzoeken
Ik probeer als onderzoeksvraag te formuleren: “Hebben dag-tot-dag veranderingen in het weer een invloed op depressieve symptomen bij mensen met psychische problemen? Ik hypothetiseer (zonder achtergrondskennis over weersinvloeden op de stemming) dat de windsnelheid, windstoten, de maximum tempertuur, de minimum temperatuur, de duur van de zonneschijn, de duur van de neerslag, de duur van de bewolking, en de luchtvochtigheid de stemming zou kunnen beïnvloeden.
Ik ga op zoek naar meteorologische data, en het liefst van dag tot dag. Via google kom ik uiteindelijk uit op de knmi site, waar je informatie van 30 verschillende weerstations in Nederland kan downloaden. Ik ga de data van tussen 2001 en 2011 downloaden, omdat we ook informatie hebben van Routine Outcome Monitoring (ROM) data in die periode. Ik wil proberen om die twee bestanden dan aan elkaar te koppelen met join om zo te kijken of en welke relaties er te vinden zijn.
Omdat Leiden geen weerstation blijkt te bevatten, kies ik wat weerstations in de buurt, en ik ben van plan om die data dan te middelen. Maar eerst even kijken of dat wel een reële afweging is. Dan hebben we wel een beetje een beeld van de weersomstandigheiden in en rond Leiden waar metingen van de stemming zijn gedaan.
De data in het textbestand KNMI is goed gedocumenteerd; de eerste regels bevatten allemaal informatie. Die zal ik niet in een dataframe moeten inladen, dus die moet ik negeren bij het inladen. Deel van deze informatie is:
BRON: KONINKLIJK NEDERLANDS METEOROLOGISCH INSTITUUT (KNMI)
STN LON(east) LAT(north) ALT(m) NAME
210: 4.430 52.171 -0.20 VALKENBURG
215: 4.437 52.141 -1.10 VOORSCHOTEN
225: 4.555 52.463 4.40 IJMUIDEN
240: 4.790 52.318 -3.30 SCHIPHOL
257: 4.603 52.506 8.50 WIJK AAN ZEE
330: 4.122 51.992 11.90 HOEK VAN HOLLAND
343: 4.313 51.893 3.50 R’DAM-GEULHAVEN
344: 4.447 51.962 -4.30 ROTTERDAM
YYYYMMDD = Datum (YYYY=jaar MM=maand DD=dag)
FHX = Hoogste uurgemiddelde windsnelheid (in 0.1 m/s)
FXX = Hoogste windstoot (in 0.1 m/s)
TN = Minimum temperatuur (in 0.1 graden Celsius)
TX = Maximum temperatuur (in 0.1 graden Celsius)
SQ = Zonneschijnduur (in 0.1 uur) berekend uit de globale straling (-1 voor <0.05 uur)
DR = Duur van de neerslag (in 0.1 uur)
NG = Etmaalgemiddelde bewolking (bedekkingsgraad van de bovenlucht in achtsten, 9=bovenlucht onzichtbaar)
UG = Etmaalgemiddelde relatieve vochtigheid (in procenten)
Inladen van de data
De eerste regels zal ik moeten overslaan. Dat kan met de skip commando, maar nog handiger is om aan te geven achter de comment statement dat alles achter een “#” genegeerd moet worden. Alleen kom ik erachter dat ik deze ‘hashtag’ voor de kolom-namen moet weghalen in een teksteditor, zodat deze dan ook direct worden ingelezen.
KNMI <- read_csv("KNMI_20131215.txt.tsv", comment = "#")
KNMI %>% glimpse()## Observations: 33,004
## Variables: 10
## $ STN <int> 210, 210, 210, 210, 210, 210, 210, 210, 210, 210, 210...
## $ YYYYMMDD <int> 20010101, 20010102, 20010103, 20010104, 20010105, 200...
## $ FHX <int> 90, 90, 80, 100, 110, 70, 60, 60, 40, 80, 80, 40, 50,...
## $ FXX <int> 150, 160, 140, 170, 190, 110, 110, 140, 70, 120, 130,...
## $ TN <int> 11, 69, 54, 67, 67, 49, 40, 25, 28, 25, -22, -36, -18...
## $ TX <int> 87, 110, 93, 91, 104, 78, 74, 68, 66, 44, 33, 52, 24,...
## $ SQ <int> 0, 25, 48, 0, 0, 36, 15, 2, 8, 2, 64, 7, 65, 10, 71, ...
## $ DR <int> 74, 2, 44, 47, 98, 30, 0, 21, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ NG <int> 7, 5, 5, 8, 8, 5, 5, 5, 6, 8, 3, 4, 1, 4, 0, 0, 1, 6,...
## $ UG <int> 91, 84, 87, 91, 94, 89, 90, 94, 94, 94, 79, 89, 85, 9...Dat is goed ingeladen, alleen heeft de datum nog niet het goede format als datum. Dat is wel vaker het geval met datums, daar moet altijd goed op gelet worden (ook bij excel). Dat passen we aan met een functie van het lubridate package, die alle datums kan omzetten. Ik verwijder daarna de kolom met de naam “YYYYMMDD”. Ik hercodeer ook alle codes van de verschillende meetstations.
Locatie van de meetstations
In de file zitten ook de locaties van de meetstations. Ik wil de meetstations in de buurt van Leiden. Die locaties zou ik kunnen visualiseren in R via een package die google maps gebruikt, namelijk ggmap. Maar eerst wil ik die data inlezen. Dat is wel even een precisiewerkje, maar wel mooi als ik het niet hoef over te typen, dat gaat uiteindelijk in de volgende stappen:
- lees de hele file nog eens in
- gebruik alleen rijen 6 t/m 13 waar de locatie data zich bevindt
- knip de data in aparte velden (op plek 1 van links, dan 15 van links, enz.)
- eerste kolom ‘even’ kan eruit worden gehaald
- definieer de longitude en latitude als getallen
meetplekken <- read.csv("KNMI_20131215.txt.tsv", header = FALSE) %>%
slice(6:13) %>%
separate(V1, sep = c(1, 15, 27, 39, 46),
c("even", "meetstation", "long", "lat", "alt", "naam")) %>%
select(-even) %>%
mutate (long = as.double(long),
lat = as.double(lat))
library(ggmap)
ned <- get_map(location = 'netherlands', zoom = 8)
ggmap(ned) + geom_point(aes(long, lat), data = meetplekken,
color = "darkblue", size = 7, alpha = 0.5)
De meetstations liggen dus mooi rondom Leiden.
Visualiseren
Laten we proberen er een figuur van te maken van alle dagelijkse data van het weer. Ik begin maar met de variabele FHX. Dat is de hoogste uurgemiddelde windsnelheid (in 0.1 m/s).
KNMI_opschoon %>%
ggplot(aes(date, FHX, color = STN)) +
geom_smooth(size=1.4, method = "loess", span = 0.1) +
ylab ("Hoogste uurgemiddelde windsnelheid (in 0.1 m/s)") 
De windsnelheid is flink sterker aan de kust (en in de winter) dan wat verder in het binnenland. Nu zou ik nog 7 andere figuren kunnen maken voor de andere variabelen. Maar ik kan ook een long dataformat maken, en dan in 1 instructie alle figuren maken.
KNMI_long <- KNMI_opschoon %>%
gather(meting, waarde, FHX:UG) %>%
glimpse()## Observations: 192,528
## Variables: 4
## $ STN <fct> Valkenburg, Valkenburg, Valkenburg, Valkenburg, Valkenb...
## $ date <date> 2002-08-02, 2002-08-03, 2002-08-04, 2002-08-05, 2002-0...
## $ meting <chr> "FHX", "FHX", "FHX", "FHX", "FHX", "FHX", "FHX", "FHX",...
## $ waarde <dbl> 50, 40, 40, 50, 70, 60, 50, 40, 50, 50, 70, 60, 40, 40,...Er zijn nu dus nog maar 4 variabelen. Daar maak ik dan in 1 instructie een figuur van met 8 facets.
Dat geeft veel informatie. Behalve windstoten en windsnelheid vertonen de 7 verschillende meetcentra ongeveer dezelfde meetwaardes. Ik zal nu de verschillende meetcentra gaan grouperen om de gemiddeldes in een database te krijgen. Die zal ik dan relateren aan de BDI data.
KNMI_gemiddeld <- KNMI_opschoon %>%
group_by(date) %>%
summarise(windsnelheid = mean(FHX, na.rm = TRUE),
windstoot = mean(FXX, na.rm = TRUE),
temp_min = mean(TN, na.rm = TRUE),
temp_max = mean(TX, na.rm = TRUE),
zonneschijn = mean(SQ, na.rm = TRUE),
neerslag_duur = mean(DR, na.rm = TRUE),
bewolking = mean(NG, na.rm = TRUE),
vochtigheid = mean(UG, na.rm = TRUE)) %>%
ungroup()Ik toon een figuur van de hoogste windsnelheid gemiddeld over de 8 meetstations rondom Leiden, waarbij ik ook de dag tot dag variatie toon. Er is dus voldoende variatie over iedere dag heen, en zulke variantie is belangrijk om een verband aan te kunnen tonen en dagen met elkaar te kunnen vergelijken.