Data science blog #4: depressie in Amerika

Data opzoeken

Data is steeds vaker verzameld met publiek geld. Die data wordt regelmatig openbaar beschikbaar gesteld, zodat nieuwsgierige onderzoekers daar analyses op kunnen uitvoeren.

Ik zal eerst op zoek gaan naar data over depressie in de USA om te downloaden. Voor Amerikaanse gegevens is NHANES data beschikbaar. Dat is een acronym voor de National Health and Nutrition Examination Survey – CDC. Daar worden met regelmaat enorme steekproeven uit de bevolking aangeschreven die worden uitgenodigd om een gezondheidsonderzoek te ondergaan. Het is beschikbaar in een SAS format, en zo’n dataset kan je met korte sasxport.get-instructies inlezen in R.

library(Hmisc)
NHANES <- sasxport.get("DPQ_I.XPT")

FALSE Processing SAS dataset DPQ_I   ..

Dit is data over de periode van 2015 tot 2016. Het is afkomstig van mensen met een leeftijd van 18 of ouder, afgenomen met een screenings-instrument voor depressie bestaande uit 9 items. Het gaat over de symptomen in de afgelopen 2 weken. De 10e vraag gaat over de impact van deze klachten. De scores van de 9 items loopt van 0 tot 3:

• 0: “not at all”
• 1: “several days”
• 2: “more than half the days”
• 3: “nearly every day”

Ik ben ten eerste benieuwd of het aantal variabelen dat erin zit klopt, en van de hoeveel mensen data van wie data is verzameld. Dat doe je met de dim instructie (voor ‘dimensions’).

dim(NHANES)

## [1] 5735   11

Dat zijn dus 5735 rijen (mensen) en 11 kolommen (variabelen). Met de glimpse-instructie krijg je een net overzicht van de data.

glimpse(NHANES)

## Observations: 5,735
## Variables: 11
## $ seqn   <int+lbl> 83732, 83733, 83734, 83735, 83736, 83737, 83741, 83...
## $ dpq010 <int+lbl> 0, 1, 0, 1, 1, 0, 0, 0, NA, 0, 0, 0, 0, 1, 1, 1, 1,...
## $ dpq020 <int+lbl> 0, 0, 0, 1, 1, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 2, 1,...
## $ dpq030 <int+lbl> 0, 0, 0, 2, 1, 0, 0, 1, NA, 0, 2, 0, 0, 1, 1, 0, 3,...
## $ dpq040 <int+lbl> 1, 0, 1, 2, 1, 0, 1, 0, NA, 0, 0, 0, 1, 1, 1, 1, 3,...
## $ dpq050 <int+lbl> 0, 1, 0, 1, 3, 0, 0, 0, NA, 0, 0, 0, 0, 1, 1, 0, 2,...
## $ dpq060 <int+lbl> 0, 0, 0, 3, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 2, 1,...
## $ dpq070 <int+lbl> 0, 0, 0, 2, 1, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 1, 1,...
## $ dpq080 <int+lbl> 0, 0, 0, 0, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 2, 1,...
## $ dpq090 <int+lbl> 0, 0, 0, 1, 0, 0, 0, 0, NA, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ dpq100 <int+lbl> 0, 0, 1, 0, 0, NA, 0, 0, NA, NA, 0, NA, 0, 0, 0, 0,...

Summary is ook een eenvoudige instructie die werkt voor de meeste variabelen in de dataframe, om de mediaan, gemiddelde, range, en percentielen te krijgen voor de continue variabelen. Zo krijg je de belangrijkste info voor elke variabele in een handomdraai. Om niet teveel informatie te tonen, selecteer ik alleen de eerste 4 variabelen.

Het is handig om gebruik te maken van de pipe-(%>%)-operator, daarmee stuur je het dataresultaat door naar de volgende bewerking, om zo een ‘pipeline’ van verschillende datamanipulaties uit te kunnen voeren. Dat verbetert de leesbaarheid van je code enorm. Daardoor is het eenvoudiger om complexe datamanipulatie uit te voeren, door de ene eenvoudige instructie de andere te laten volgen. De volgende drie acties worden opeenvolgende uitgevoerd:

1. neem het dataframe ‘NHANES’
2. selecteer de kolommen 2 tot en met 5
3. geef de ‘summary’ statistiek weer van deze 4 variabelen

NHANES %>% 
  select (2:5) %>% 
  summary()

##      dpq010           dpq020           dpq030           dpq040      
##  Min.   :0.0000   Min.   :0.0000   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.0000   Median :0.0000   Median :1.0000  
##  Mean   :0.4225   Mean   :0.3466   Mean   :0.6172   Mean   :0.7836  
##  3rd Qu.:1.0000   3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:1.0000  
##  Max.   :9.0000   Max.   :7.0000   Max.   :9.0000   Max.   :7.0000  
##  NA's   :571      NA's   :571      NA's   :571      NA's   :573

Dat is equivalent met deze instructie die van binnen naar buiten gelezen moet worden (zonder de ‘pipe’):

summary(select(NHANES,2:5))

Dus er zijn voor de meeste variabelen 571 tot 573 NA’s (not available) ofwel missing values.

Over de eerste variabele staat in het online codebook het volgende geschreven:

DPQ010 – Have little interest in doing things
Variable Name: DPQ010
SAS Label: Have little interest in doing things
English Text: Over the last 2 weeks, how often have you been bothered by the following problems: little interest or pleasure in doing things?
Would you say…

Dus deze vraag gaat over een key symptoom van depressie, namelijk over interesseverlies of gebrek aan plezier.

De antwoordcategorieën zijn:

0 Not at all
1 Several days
2 More than half the days
3 Nearly every day
7 Refused
9 Don’t know

Laat ik daar eens een bar-plot van maken.

NHANES$dpq010 %>% 
  ftable() %>% 
  barplot(main="dpq010: little interest or pleasure", 
    xlab="Antwoordcategorie")

Dat ziet er al aardig uit, maar het is niet zo handig dat je geen informatie krijgt te zien over de labels op de x-as. Deze data kan verder met de Data Documentation, Codebook, and Frequencies worden opgeschoond en gelabeld.

Data opschonen

Ik denk dat het wat beter werkt om de assen om te draaien (met coord_flip()), en het op te maken in het handige data-visualisatie package ggplot2.

NHANES %>%
    mutate (var1 = dpq010 %>% as.factor() %>% 
            recode("0" = "Not at all", "1" = "Several days", 
                   "2" = "More than half the days", 
                   "3" = "Nearly every day", 
                   "7" = "Refused", "9" = "Don't know")) %>% 
   ggplot(aes(var1)) + geom_bar() + 
   coord_flip() + xlab("") +
   ggtitle(attr(NHANES$dpq010, "label")) +
   ylab("Aantal") +
   theme_economist()

Dat wordt inderdaad overzichtelijker. Nu kan ik deze vrij lange instructie 8 maal herhalen voor elke variabele, maar een basisregel in programmeren is om nooit iets te herhalen. Maak daar dan een functie voor en roep die aan. R is een ‘functional programming’ taal. Het is heel eenvoudig om de voorgaande instructie om te zetten in een functie, en deze toe te kennen aan een naam voor die functie met <-. Deze omgekeerde pijl is het toekenningssymbool.

maak_figuur <- function(var_naam) {
   NHANES %>%
    mutate (var1 = NHANES[[var_naam]] %>% as.factor() %>% 
            recode("0" = "Not at all", "1" = "Several days", 
                   "2" = "More than half the days", 
                   "3" = "Nearly every day", 
                   "7" = "Refused", "9" = "Don't know")) %>% 
   ggplot(aes(var1)) + geom_bar() + 
   coord_flip() + xlab ("") +
   ggtitle(attr(NHANES[[var_naam]], "label")) +
   ylab("Aantal") +
   theme_economist()
}

Nu kan voortaan met één korte instructie een figuur worden gemaakt van elk van de 9 variabelen.

maak_figuur("dpq060")

maak_figuur("dpq050")

maak_figuur("dpq090")

In deze aselecte steekproef komt suïcidaliteit dus maar weinig voor.

Data ordenen

Veel data kan op verschillende manieren in een database worden geordend. Om je ‘data-cleaning’ eenvoudig en effectief te laten verlopen is een opgeschoonde dataset nodig. Dit houdt meestal in dat je data in een lang (long) format staat. Zo’n lange gegevensset is gemakkelijk te manipuleren, te modelleren, en te visualiseren. Het heeft een specifieke structuur: elke variabele is een kolom en elke waarneming is een rij.

De volgende twee datasets zijn in essentie hetzelfde, maar de ene is wide (breed) en de ander long (of tall). Deze fictieve tabellen bevatten informatie over depressie die wel of niet aanwezig is in opeenvolgende jaren. Via een enkele spread of een gather instructie zijn deze twee ‘formats’ in elkaar om te zetten.

In de long data-frame zijn geen rijen voor de missing values. In de wide data-frame daarentegen zie je dat sommige deelnemers missing data (NA’s) hebben voor bepaalde jaren.

Hadley Wickham en collega’s hebben veel packages gemaakt met gereedschappen die nodig zijn om deze ‘tidy data’ aan te pakken. Het is vrij eenvoudig om met de functies van de dplyr package deze data te manipuleren en er data manipilatie ermee uit te voeren. Dat pakket zit in de tidyverse. Nadat je eerst dit pakket hebt geinstalleerd, zal je het met de library instructie kunnen activeren voor jouw R-sessie.

install.packages("tidyverse")
library(tidyverse)

Netwerk analyse

Nu kunnen we op deze data bijvoorbeeld een netwerkanalyse doen met het qgraph package. Prof. dr. Denny Borsboom en zijn onderzoeksgroep heeft hier veel over geschreven. Ik zal proberen een netwerk te maken, volgens de online instructies beschreven door Sacha Epskamp in Network Analysis in R Cookbook.

De Spearman’s correlation tabel in de vorm van een heatmap kan gemaakt worden met het package corrplot en toont bijvoorbeeld dat item 9 (‘suïcidaliteit’) hier slechts zwak geassocieerd is met de andere symptomen. Echter hier wordt niet gecorrigeerd (‘partial correlations’) voor andere associaties in het netwerk van symptomen.

library(corrplot)
corMatrix1 <- netwerk %>% 
  cor(use = "pairwise.complete.obs",method = c("spearman")) %>%
  corrplot(method=c("color"), order="hclust",addrect = 5, tl.col = "black")

Zo’n adjustering gebeurt wel in deze netwerkanalyse:

library(qgraph)
corMat <- cor_auto(netwerk)
qgraph(corMat, graph = "pcor", layout = "spring", threshold = "bonferroni",
                     sampleSize = nrow(netwerk), alpha = 0.05)

Item 6 “Feeling bad about yourself” is dus bijvoorbeeld het sterkst gelinkt aan item 2 “Feeling down, depressed, or hopeless” en aan item 9 “Thought you would be better off dead”.

Een mooi artikel over dit onderwerp van dr. Angélique Cramer, waar ik ook aan meegeschreven heb, is getiteld: Major depression as a complex dynamic system.

Voor commentaar, suggesties, of een interessante vraagstelling:
Contact
Afdeling Psychiatrie, LUMC