Data science blog #5: wordt de stemming beïnvloed door het weer?

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:

1. lees de hele file nog eens in
2. gebruik alleen rijen 6 t/m 13 waar de locatie data zich bevindt
3. knip de data in aparte velden (op plek 1 van links, dan 15 van links, enz.)
4. eerste kolom ‘even’ kan eruit worden gehaald
5. 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.

Ik zal nu het seizoen en de maand toevoegen, gebruik makend van een mooie functie die iemand op internet heeft gezet. Dan kan ik de modellen later corrigeren voor het seizoen of maand, want ik ben niet benieuwd naar de winterdepressie, maar meer naar de korte-termein effecten van weersinvloeden op de verschillende symptoomclusters. Uiteindelijk lijkt de maand correctie het beste te zijn. We kijken dus niet naar de seizoenseffecten op de stemming, maar naar de dag tot dag variaties van de verschillende meteorologische parameters.

geef_seizoen <- function(input.date){
  numeric.date <- 100*month(input.date)+day(input.date)
  cuts <- base::cut(numeric.date, breaks = c(0,319,620,921,1220,1231)) 
  levels(cuts) <- c("Winter","Lente","Zomer","Herfst","Winter")
  return(cuts)
}

KNMI_gemiddeld <- KNMI_gemiddeld %>% 
  mutate(seizoen = geef_seizoen(date) %>% as.factor(),
         maand = month(date) %>% as.factor())

Nu zal ik de geanonimiseerde Routine Outcome Monitoring data inlezen van de Beck Depression Scale (BDI), met datum. Deze zal ik mergen (via join) met de KNMI_gemiddeld file.

BDI <- read_sav("BDIerik.sav") %>% 
  select(-ID, -c(BDI01:BDI21)) %>% 
  rename(date = DATE) %>% 
  inner_join(KNMI_gemiddeld)

## Joining, by = "date"

Dat geeft een dataframe met 22972 rijen en 15 kolommen. Nu kunnen eens onderzoeken welke associaties er zijn tussen de 3 symptoomclusters en weersinvloeden.

Dit is een analyse met 1.280.888 observaties (in een long data format). Ik heb de top 8 verbanden laten zien, met de sterkste effecten van de 32 (= 8 weet-metingen * 4 schalen) analyses die ik heb gedaan. Ik heb steeds voor de maand geadjusteerd (12 maanden is 11 vrijheidsgraden).

Het lijkt gelukt om weerdata te downloaden en te linken aan een ander databestand. Enkele effecten zijn statistisch significant. De eerste bevindingen lijken er wel op te wijzen dat hogere (minimum en maximum) temperaturen op de dag met meer depressieve symptomen gerelateerd zijn. Daartegenover is een regenachtige dag juist geassocieerd met iets minder klachten. Zouden we dus moeten concluderen dat men zich wat beter voelt bij een druilerige dag met typisch Hollands weer?

Maar de statistische significante associaties zijn zeer zwak (verklaarde variantie van kleiner dan 0.1%; de beta-estimates zijn de gestandaardseerde coëfficienten, en het kwadraat ervan is een maat voor de verklaarde variantie). Dit kan erop wijzen dat er geen relaties van betekenis zijn tussen dag-tot-dag veranderingen in het weer en depressieve symptomen bij mensen met psychische problemen. Aan de andere kant kunnen deze associaties ook een flinke onderschatting zijn, want eigenlijk zou je binnen één persoon sequentiële metingen van de stemming moeten doen (van opeenvolgende dagen), om veel sterkere statistische mogelijkheden te hebben om dit te analyseren.

Eerdere analyses vonden geen effecten in relatie tot meteorologische data, maar die analyse was kleiner dan de bovenstaande (met 556 en 206 deelnemers, en 8.316 observaties).

Voor commentaar, suggesties, of een interessante vraagstelling: