Op deze pagina vind je een demonstratie van een statistische techniek aan de hand van een voorbeeld.
Meer informatie over hoe je deze pagina kan gebruiken vind je in deze handleiding.
De analyse gebeurt met behulp van R en RStudio. Een inleiding tot deze software vind je hier.
Het voorbeeld op deze pagina is geïnspireerd door een studie van Vanheule, Desmet, Groenvynck, Rosseel & Fontaine (2008). Om praktische en didactische redenen werden enkele aanpassingen aangebracht.
Veel verschijnselen zijn niet direct meetbaar. Zeker in de gedragswetenschappen is dat het geval. Denk bijvoorbeeld aan depressie1: er is geen meetinstrument of “lat” voorhanden waarmee je iemands depressie rechtstreeks kan meten. Variabelen die niet rechtstreeks meetbaar zijn noemen we latente variabelen.
In factoranalyse gaat men ervan uit dat zulke latente variabelen toch meetbaar zijn door verschillende indicatoren te combineren. Dat zijn variabelen die wel rechtstreeks meetbaar zijn (ook manifeste variabelen genoemd). Elke indicator meet een aspect of een facet van de onderliggende latente variabele.
Een heel courant voorbeeld van zulke indicatoren zijn afzonderlijke vragen (“items”) uit een vragenlijst.
Je mag er natuurlijk niet zomaar vanuit gaan dat een reeks indicatoren samen een goede meting zijn voor een bepaalde latente variabele. Dat moet je eerst nagaan door de samenhang tussen indicatoren te onderzoeken. Het basisidee van factoranalyse is dat indicatoren die eenzelfde latente variabele meten, gecorreleerd zullen zijn (of anders gezegd, dat ze gedeelde variantie zullen hebben). Bijvoorbeeld personen met een hoge mate van depressiviteit zullen vermoedelijk vrij consistent hoog scoren voor een vraag die peilt naar hun gevoel van wanhoop, hoog voor een vraag naar hun zwartgallige gedachten, laag voor hun aantal activiteiten buitenshuis enzovoort. Mensen die niet of weinig depressief zijn zullen vermoedelijk vrij consistent omgekeerd scoren.
Deze pagina demonstreert hoe je een confirmatorische factoranalyse (CFA) kan uitvoeren. Daarin is het de bedoeling om na te gaan of een vooropgestelde samenhang tussen indicatoren inderdaad in een dataset terug te vinden is. Je hebt in dit geval dus vooraf al een idee van
Het doel van CFA is om na te gaan of een vooropgestelde factorstructuur (= een aantal factoren, elk met hun eigen set van indicatoren) daadwerkelijk terug te vinden is in de data.2
Die confirmatorische analyse staat tegenover een exploratieve aanpak (EFA). Daarin start je enkel vanuit een reeks indicatoren, zoals een vragenlijst. Hierbij vraag je je af welke indicatoren samen horen en dus mogelijk een factor vormen. In zo’n geval heb je vooraf geen hypothese over hoeveel en welke factoren uit de analyse tevoorschijn zullen komen.
De studie van Vanheule et al. (2008) voert onder meer3 een confirmatorische factoranalyse uit op de Beck Depression Inventory-II (BDI). Die vragenlijst meet de mate van depressie. De onderzoekers willen nagaan of er bepaalde subdimensies of deelaspecten van depressie kunnen worden gevonden in de BDI. Die subdimensies zijn in deze studie de latente variabelen.
Zoals steeds in CFA hebben de onderzoekers op voorhand hypothesen over hoeveel en welke factoren gemeten worden met de vragenlijst. Concreet wordt nagegaan of er een somatisch-affectieve (som.aff
) en een cognitieve (cog
) factor terug te vinden zijn in de data. Deze latente variabelen zouden worden gemeten door respectievelijk 9 en 12 items.
Op het diagram hieronder zie je welke items zouden samenhoren en hoeveel/welke factoren men denkt terug te vinden in de data.
\(\eta_1\) en \(\eta_2\) stellen de twee latente variabelen voor. \(Y_i\) zijn alle 21 indicatoren. Bij elke indicator hoort ook een \(\epsilon_i\). Die slaat op de unieke variantie in elke indicator die niet gedeeld wordt met de andere indicatoren.
De dataset bdiKlinisch
bevat gegevens van 21 variabelen geobserveerd bij 404 patiënten uit een Belgisch centrum voor geestelijke gezondheidszorg.
Merk op dat je bij CFA kan vertrekken van de ruwe data of van de covariantiematrix, maar niet van de correlatiematrix.
De dataset kan je inladen met read.csv()
. De data kan je best meteen in een object bdiKlinisch
onderbrengen zodat je die later makkelijk opnieuw kan oproepen.
bdiKlinisch <- read.csv("https://statlas.ugent.be/datasets/bdi.items.klinisch.csv")
Met str()
krijg je een opsomming van alle variabelen in de dataset.
str(bdiKlinisch)
'data.frame': 404 obs. of 21 variables:
$ BDI1 : int 0 0 1 1 2 1 1 1 1 0 ...
$ BDI2 : int 0 1 0 4 2 1 1 1 2 0 ...
$ BDI3 : int 1 0 0 2 0 2 2 2 2 1 ...
$ BDI4 : int 0 0 0 2 3 0 2 2 1 1 ...
$ BDI5 : int 0 0 0 1 2 1 2 0 1 1 ...
$ BDI6 : int 3 0 0 0 0 3 3 3 3 0 ...
$ BDI7 : int 0 0 1 2 2 2 2 2 2 0 ...
$ BDI8 : int 0 1 0 3 0 0 2 1 1 1 ...
$ BDI9 : int 0 0 0 1 0 1 1 1 1 0 ...
$ BDI10: int 1 0 0 2 3 3 3 2 3 3 ...
$ BDI11: int 0 1 0 1 1 0 1 3 1 0 ...
$ BDI12: int 0 0 0 1 3 1 1 2 2 0 ...
$ BDI13: int 1 1 0 3 3 1 3 2 1 0 ...
$ BDI14: int 0 0 0 3 1 1 2 1 2 0 ...
$ BDI15: int 0 0 0 3 3 1 2 3 3 1 ...
$ BDI16: int 0 1 1 2 2 1 3 3 2 1 ...
$ BDI17: int 0 0 0 2 1 0 2 2 1 0 ...
$ BDI18: int 0 0 0 2 2 1 2 3 1 1 ...
$ BDI19: int 0 0 0 3 2 2 2 1 2 1 ...
$ BDI20: int 0 0 0 3 2 1 2 2 3 1 ...
$ BDI21: int 0 1 0 3 3 2 2 0 3 0 ...
In de output van str()
zie je inderdaad dat er 21 variabelen zijn met telkens 404 observaties.
lavaan
Voor de analyse van een CFA-model gebruik je het package lavaan
, wat staat voor latent variable analysis.
install.packages('lavaan') # eenmalig installeren
library(lavaan) # het package laden bij de start van elke sessie
Elke variabele BDI1
, BDI2
, enzovoort is mogelijk een indicator voor een onderliggende factor. Die plaats je bij elke latente variabele in overeenstemming met je hypothesen:
mijnModel <- 'cog =~ BDI1 + BDI2 + BDI3 + BDI5 + BDI6 + BDI7 + BDI8 + BDI9 + BDI14
som.aff =~ BDI4 + BDI10 + BDI11 + BDI12 + BDI13 + BDI15 + BDI16 +
BDI17 + BDI18 + BDI19 + BDI20 + BDI21
'
Deze code weerspiegelt dus de hypothesen die de onderzoekers bij voorbaat hadden over welke indicatoren welke latente variabele meten. Het gelijkheidsteken =
in combinatie met een tilde ~
kan je lezen als “wordt gemeten door”. Dus kan je de eerste regel hierboven begrijpen als “de latente variabele cog
wordt gemeten door de indicatoren BDI1
, BDI2
, BDI3
enzovoort”4.
Let op de enkele aanhalingstekens '
vooraan en achteraan! De modelsyntax vormt immers een “string” (een stuk tekst) in R.
De functie cfa()
voert alle berekeningen uit die je nodig hebt. Deze functie heeft de specificatie van het model nodig (object mijnModel
) en de data.
fit <- cfa(mijnModel, data=bdiKlinisch, se='robust', test='Satorra-Bentler')
Een assumptie bij het schatten is dat de indicatoren \(Y_i\) normaal verdeeld zijn. Als dat niet het geval is, dan zullen de resultaten vertekend zijn. Voor de puntschattingen van de factorladingen is de vertekening klein. Voor de standaardfouten (“standard errors”) is er wel een betekenisvol verschil. Met het argument se='robust'
bouw je een correctie in: de standaardfouten worden gecorrigeerd voor non-normaliteit van de data. Een gelijkaardige redenering gaat op voor de teststatistiek \(\chi^2\). Die wordt gecorrigeerd voor non-normaliteit met het argument test='Satorra-Bentler'
.
Ook goed om te weten: cfa()
vult je modelspecificatie automatisch verder aan. Zo wordt er standaard voor gekozen dat de factoren kunnen correleren. Ook worden de residuele varianties van de indicatoren geschat. In de meeste gevallen zijn dat goede keuzes. Mocht je het toch anders willen5, dan kan je dit expliciet aangeven in jouw modelspecificatie. Bijvoorbeeld, volgende code legt vast dat de correlatie tussen factoren som.aff
en cog
toch 0 moet zijn.
mijnModel <- 'cog =~ BDI1 + BDI2 + BDI3 + BDI5 + BDI6 + BDI7 + BDI8 + BDI9 + BDI14
som.aff =~ BDI4 + BDI10 + BDI11 + BDI12 + BDI13 + BDI15 + BDI16 +
BDI17 + BDI18 + BDI19 + BDI20 + BDI21
som.aff ~~ 0*cog
'
In de rest van deze demonstratie zullen we wel toelaten dat som.aff
en cog
correleren.
Om de resultaten te zien gebruik je de functie summary()
. Het object fit
voeg je toe als eerste argument. Met de andere argumenten kan je verder bepalen wat je precies in de output wil zien. (De output zelf wordt verderop besproken.)
summary(fit, fit.measures = TRUE, standardized = TRUE)
Enkele keuzes die je kan maken bij de summary()
-functie zijn:
fit.measures = TRUE
Indien fit.measures = FALSE
- wat standaard het geval is - wordt enkel een \(\chi^2\)-toets uitgevoerd als maat voor de fit van het model. Deze toets is niet altijd voldoende om de model fit te beoordelen, zeker bij hele grote steekproeven.
Het is aangeraden om bijkomende fitmaten te laten berekenen en te rapporteren in je publicatie, zoals CFI/TLI, RMSEA en SRMR. Al deze fitmaten zijn op een verschillende manier gevoelig voor misspecificatie van het model. Voeg deze fitmaten toe door eenvoudigweg het argument fit.measures = TRUE
te gebruiken. De lezer van je onderzoek kan zo ook zelf een oordeel vellen over de fit.
Meer uitgebreide uitleg over deze fitmaten vind je in:
Kline, R. B. (2015). Principles and practice of structural equation modeling (Fourth Edition). New York: Guilford Press.
standardized = TRUE
Het is niet zomaar mogelijk om de factorladingen (zie verderop) onderling te vergelijken en om uitspraken te doen over welke groter of kleiner is. Dat komt omdat de getalwaarde afhangt van de meeteenheid van de variabelen.
Met het argument standardized = TRUE
krijg je een extra kolom Std.all
in de output met schattingen voor gestandaardiseerde parameters (factorladingen, covarianties tussen factoren en residuele varianties van de indicatoren).
Wat kan je leren uit de output van summary(fit, fit.measures = TRUE, standardized = TRUE)
?
summary()
lavaan 0.6-11 ended normally after 35 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 43
Used Total
Number of observations 367 404
Model Test User Model:
Standard Robust
Test Statistic 531.058 472.044
Degrees of freedom 188 188
P-value (Chi-square) 0.000 0.000
Scaling correction factor 1.125
Satorra-Bentler correction
Model Test Baseline Model:
Test statistic 3211.683 3207.866
Degrees of freedom 210 210
P-value 0.000 0.000
Scaling correction factor 1.001
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.886 0.905
Tucker-Lewis Index (TLI) 0.872 0.894
Robust Comparative Fit Index (CFI) 0.894
Robust Tucker-Lewis Index (TLI) 0.881
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -9293.801 -9293.801
Loglikelihood unrestricted model (H1) -9028.272 -9028.272
Akaike (AIC) 18673.601 18673.601
Bayesian (BIC) 18841.532 18841.532
Sample-size adjusted Bayesian (BIC) 18705.109 18705.109
Root Mean Square Error of Approximation:
RMSEA 0.071 0.064
90 Percent confidence interval - lower 0.063 0.057
90 Percent confidence interval - upper 0.078 0.071
P-value RMSEA <= 0.05 0.000 0.000
Robust RMSEA 0.068
90 Percent confidence interval - lower 0.060
90 Percent confidence interval - upper 0.076
Standardized Root Mean Square Residual:
SRMR 0.056 0.056
Parameter Estimates:
Standard errors Robust.sem
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
cog =~
BDI1 1.000 0.529 0.681
BDI2 1.283 0.103 12.429 0.000 0.679 0.704
BDI3 1.200 0.109 11.006 0.000 0.635 0.642
BDI5 1.020 0.098 10.363 0.000 0.540 0.643
BDI6 1.134 0.127 8.925 0.000 0.600 0.473
BDI7 1.208 0.100 12.109 0.000 0.639 0.715
BDI8 1.140 0.113 10.047 0.000 0.603 0.595
BDI9 0.884 0.078 11.348 0.000 0.468 0.606
BDI14 1.376 0.113 12.191 0.000 0.728 0.741
som.aff =~
BDI4 1.000 0.663 0.740
BDI10 0.818 0.094 8.751 0.000 0.542 0.441
BDI11 0.593 0.068 8.662 0.000 0.393 0.458
BDI12 1.141 0.066 17.182 0.000 0.756 0.744
BDI13 1.040 0.079 13.100 0.000 0.689 0.617
BDI15 0.939 0.061 15.395 0.000 0.622 0.716
BDI16 0.708 0.069 10.279 0.000 0.469 0.523
BDI17 0.789 0.073 10.840 0.000 0.523 0.556
BDI18 0.838 0.078 10.790 0.000 0.555 0.516
BDI19 0.820 0.064 12.868 0.000 0.543 0.647
BDI20 1.079 0.076 14.204 0.000 0.715 0.710
BDI21 0.819 0.081 10.107 0.000 0.543 0.478
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
cog ~~
som.aff 0.292 0.032 9.239 0.000 0.833 0.833
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.BDI1 0.323 0.033 9.932 0.000 0.323 0.536
.BDI2 0.470 0.042 11.181 0.000 0.470 0.505
.BDI3 0.573 0.045 12.857 0.000 0.573 0.587
.BDI5 0.412 0.035 11.648 0.000 0.412 0.586
.BDI6 1.245 0.075 16.616 0.000 1.245 0.776
.BDI7 0.390 0.036 10.982 0.000 0.390 0.489
.BDI8 0.664 0.050 13.298 0.000 0.664 0.646
.BDI9 0.378 0.039 9.802 0.000 0.378 0.633
.BDI14 0.434 0.037 11.667 0.000 0.434 0.450
.BDI4 0.362 0.028 12.959 0.000 0.362 0.452
.BDI10 1.218 0.234 5.197 0.000 1.218 0.805
.BDI11 0.582 0.046 12.637 0.000 0.582 0.790
.BDI12 0.462 0.041 11.264 0.000 0.462 0.447
.BDI13 0.773 0.062 12.566 0.000 0.773 0.620
.BDI15 0.369 0.028 13.188 0.000 0.369 0.488
.BDI16 0.585 0.044 13.226 0.000 0.585 0.727
.BDI17 0.612 0.049 12.497 0.000 0.612 0.691
.BDI18 0.849 0.063 13.506 0.000 0.849 0.734
.BDI19 0.410 0.032 12.808 0.000 0.410 0.581
.BDI20 0.504 0.047 10.735 0.000 0.504 0.496
.BDI21 0.994 0.074 13.527 0.000 0.994 0.771
cog 0.280 0.041 6.889 0.000 1.000 1.000
som.aff 0.439 0.045 9.840 0.000 1.000 1.000
Of het model in zijn geheel goed fit met de data kan je achterhalen met behulp van de verschillende fitmaten.
Je krijgt twee kolommen met fitmaten. De linkerkolom (“Standard”) toont de fitmaten zonder correctie voor normaliteit; de rechterkolom (“Robust”) toont de fitmaten na correctie voor normaliteit. Het zijn deze gecorrigeerde fitmaten in de rechterkolom die je moet interpreteren en rapporteren.
Let op! De \(\chi^2\)-toets onder Model Test User Model
interpreteren is hier anders dan wat je misschien intuïtief zou denken. Hier stelt de nulhypothese dat er een goede fit is. Dus als je een p-waarde kleiner dan 0.05 ziet, dan is de conclusie dat je de nulhypothese van goede fit moet verwerpen. Dat is in ons model het geval.
fit.measures = TRUE
heb je naast de \(\chi^2\)-toets nog andere fitmaten opgevraagd. Er bestaat geen consensus over de precieze waarde vanaf dewelke je zeker kan spreken van een goede fit. Vuistregels voor deze fitmaten die in de literatuur voorkomen zijn:
In de rechterkolom met robuuste fitmaten zie je voor de CFI/TLI en voor de RMSEA nog eens twee “versies”. De voorkeur gaat hier telkens naar de onderste van de twee, dus degene vermeld bij “Robust Comparative Fit Index”, “Robust Tucker-Lewis Index” en “Robust RMSEA”.
Hier zijn de waarden voor deze fitmaten ook niet bijster goed.6
In principe betekent dit dat je niet verder kan gaan met de analyse van de afzonderlijke parameters. Je kan bijvoorbeeld de factorladingen niet beoordelen. Je zou eerst op zoeken moeten gaan naar een beter passend model. Om de demonstratie verder te kunnen zetten zullen we er nu van uitgaan dat er wél een goede fit is.
Als de fit van je model niet goed is, kan je eventueel de functie modificationIndices()
gebruiken om suggesties te krijgen voor verbeteringen. In de kolom mi
lees je welke bijkomende parameters een substantiële impact zouden hebben op de fit.7
modificationIndices(fit,
sort=TRUE, # Sorteer van grootste naar kleinste effect op de model fit (kolom 'mi')
maximum.number=5L # Toon enkel de eerste vijf
)
modificationIndices()
lhs op rhs mi epc sepc.lv sepc.all sepc.nox
260 BDI15 ~~ BDI20 74.697 0.224 0.224 0.520 0.520
50 cog =~ BDI13 34.931 1.301 0.688 0.616 0.616
236 BDI11 ~~ BDI17 25.377 0.164 0.164 0.274 0.274
133 BDI5 ~~ BDI13 20.633 0.144 0.144 0.254 0.254
106 BDI3 ~~ BDI5 16.684 0.115 0.115 0.236 0.236
De factorladingen geven informatie over de samenhang tussen de indicatoren. Je kan ze interpreteren als regressiecoëfficiënten in een model waarbij de indicator de afhankelijke variabele is en de factor de onafhankelijke variabele.
\[y_1 = \lambda_{11}\eta_1 + \epsilon_1\]
Onder Latent Variables
vind je in de kolom Std.all
de gestandaardiseerde factorladingen. Deze kan je met elkaar vergelijken om te bepalen welke groter en kleiner zijn.
Een goede vuistregel stelt dat een gestandaardiseerde factorlading boven 0.70 moet liggen om de indicator als een degelijke meting te beschouwen. Boven 0.80 spreken we van een echt goede meting. In ons model zijn er maar weinig indicatoren die aan deze vuistregels voldoen.
Onder Covariances
vind je in de kolom Std.all
de correlatie tussen de factoren. Hier blijken de twee factoren in ons model sterk te correleren: 0.833. Het is dus erg twijfelachtig dat de twee constructen (de cognitieve en de somatisch-affectieve component van depressie) wel goed uit elkaar te houden zijn.
Onder Variances
vind je in de kolom Std.all
de schattingen van de residuele variantie van elke indicator. Dat is de variantie van de indicator die niet gedeeld wordt met de andere indicatoren. Het vormt het spiegelbeeld van de factorladingen: als de factorlading hoog is, dan moet de residuele variantie laag zijn.
Bijgevolg heeft het rapporteren van de residuele varianties hier geen meerwaarde. Het is wel aangeraden om na te gaan of er in de kolom Estimate
geen negatieve varianties zijn. Dat zou immers een symptoom kunnen zijn van slechte fit, ook al zijn de fitmaten gunstig.
Met het package lavaanPlot
kan je een CFA-model visualiseren.
install.packages("lavaanPlot") # eenmalig het package installeren
library(lavaanPlot) # package laden voor gebruik
lavaanPlot(model = fit,
node_options = list(shape = "box", fontname = "Helvetica"),
edge_options = list(color = "grey"),
)
Een alternatief is het package semPlot
.
install.packages('semPlot') # eenmalig het package installeren
library(semPlot) # package laden voor gebruik
semPaths(fit, layout='circle')
Vanheule S., Desmet M., Groenvynck H., Rosseel Y. & Fontaine J. (2008). The factor structure of the Beck Depression Inventory-II. An evaluation. Assessment 15 (2). doi: 10.1177/1073191107311261
Andere voorbeelden zijn intelligentie, stress, welbevinden,…↩︎
Een CFA kan niet helpen bij het bepalen van de juiste benaming van factoren. De analyse kan enkel bepalen of indicatoren gezamenlijk een “iets” meten. Het is aan jou, de vakspecialist, om te bepalen of die indicatoren metingen zijn van depressie of verveling of apathie of nog iets anders.↩︎
In Vanheule et al. (2008) toetst men veel verschillende modellen met 2, 3 of 4 factoren, met ook andere combinaties van items. Het ultieme doel is om uit te zoeken welk model het best past bij de data.↩︎
Of dat werkelijk het geval is zal natuurlijk moeten blijken uit de analyse. Op dit moment is het nog maar een hypothese!↩︎
Om theoretische, inhoudelijke redenen.↩︎
Mochten deze drie fitmaten allemaal wijzen op een goede fit, dan zou dat de conclusie zijn die we trekken - ook al wijst de \(\chi^2\)-toets eerder in de richting van slechte fit. Zie ook eerder.↩︎
Bij de verbanden die al in je model zitten, staat er logischerwijs een waarde 0 in deze kolom.↩︎