Algoritme voor een longitudinale, representatieve steekproef

Nemen we twee types steekproeven:

• Longitudinaal: steekproef in jaar t, eenheden opgevolgd door de tijd
• Cross-sectioneel: steekproef in jaar t, steekproef in jaar t+1, t+2, ...

Men zou denken dat deze benaderingen onverzoenbaar zijn. Nochtans is er een algoritme dat toelaat steeds representatief te blijven, én over longitudinale gegevens te beschikken.

Dit voorbeeld kan toegepast worden op de selectie van werknemers, waar de populatie afhankelijk is van de werkloosheid.

Jaar t

Men neemt een steekproef, bijvoorbeeld 25% van de populatie.

Jaar t+1

Er is uitstroom uit de populatie. Dit is geen probleem: aangezien de steekproef aselect is verwachten we dat de sample een gelijke uitstroom heeft. We behouden eenvoudigweg de cases uit jaar t.

Er is ook instroom. Hieruit moeten we opnieuw 25% gaan selecteren.

Jaar t+2

Uitstroom en instroom zullen zich telkens herhalen zoals in jaar t+1.

Er is ook herinstroom in de populatie. Dit betekent voor onze steekproef dat we alle ooit geselecteerde eenheden moeten opvolgen, en opnieuw opnemen als ze terug in de populatie komen. Zij blijven namelijk altijd 25% van het gelijkaardig segment (de herinstroom) vertegenwoordigen. De steekproef blijft dus representatief.

Heteroskedasticiteit

deze post is herzien op 11-2-2013

Theorie

Homoskedasticiteit is een centrale assumptie in OLS, en betekent de fout niet varieert volgens de waarden van het model, dat is wenselijk (als in 'homofilie'). Bij heteroskedasticiteit wordt hier dus aan gezondigd. Deze assumptie kan met andere technieken (GMM, ML) weggenomen worden.

Testen

Er zijn twee populaire testen voor heteroskedasticiteit: Brausch-Pagan en White. Het idee is gelijkaardig: bij Brausch-Pagan gaan we de gekwadrateerde fout $ ( y - \hat{y})^2 $ (de restvariantie) regresseren op een quadratische functie van de voorspelde waarde $ \hat{y} $.

Bij White regresseren we op alle andere (unieke) verklarende variabelen in het model. De R² van dit model geeft aan of er heteroskedasticiteit is. Idealiter is R² hier zeer laag. De teststatistiek voor de nulhypothese is:

nR² -> Chi²(m), met m = aantal variabelen in het model

Referenties

• White test: Hayashi 2000, p. 131

Factoranalyse

Een zeer korte handleiding voor factoranalyse:

Waarom factoranalyse

Factoranalyse is een datareductiemethode, waarbij gezocht wordt naar latente variabelen die een set van variabelen samenvatten. Het vertrekt daarom van de correlatiematrix van deze geobserveerde variabelen.

Voorbeeld: leerlingen leggen op het eind van het jaar twaalf testen af. Een exploratieve factoranalyse zou aan het licht kunnen brengen dat dit eigenlijk slechts drie vaardigheden toetst.

1. Taalvakken: Nederlands, Frans, Engels, Duits
2. Analytische vakken: wiskunde, natuurkunde, scheikunde, logica
3. Kennisvakken: geografie, biologie, geschiedenis, esthetica

Een confirmatorische factoranalyse werk op de omgekeerde manier: 

Let op!

• Controleer eerst de correlatiematrix: enkel variabelen die een correlatie boven .30 met andere variabelen hebben dienen te worden weerhouden.
• Gebruik geen factoranalyse als de correlatie tussen variabelen niet aan een latente factor te wijten is maar bijvoorbeeld aan de distributie van de populatie op irrelevante kenmerken.

Scree plot

Assumpties

Dezelfde als voor correlaties: normaal verdeelde, continue variabelen. Heb je ook of enkel ordinale of binaire data, dan moet je vluchten naar polychorische principale componentenanalyse (het trage en onoverzichtelijke polychoricpca in Stata of de complexe explorative factoranalyse in Mplus).

Eigenwaardes

Als een vuistregel gebruikt men Kaisers criterium (Kaiser's criterion) dat zegt dat eigenwaardes groter dan 1 moeten zijn. Alternatief kan je kijken naar relatieve sterke dalingen van de eigenwaarde, naar de cumulatieve voorspelde variantie, naar de factorsterkte (Mplus: > .9) of naar bepaalde significantiematen (Mplus: RMSE). 

Rotatie

Om de geobserveerde variabelen te scoren op een bepaalde matrix zijn er verschillende mogelijkheden. Orthogonale rotaties of oblique rotaties. Voor zover ik het begrijp zijn de dimensies in een orthogonale rotatie niet gecorreleerd (denk aan een kruis met rechte hoeken), terwijl dit in oblique rotaties wel het geval is (een kruis met scherpe en stompe hoeken). De vergelijking met kruisen gaat natuurlijk enkel op in twee- of driedimensionele omgevingen. De richting van de rotatie wordt bepaald door de te maximaliseren of minimaliseren afstand:

• Varimax: maximaliseert per factor de variantie van de ladingen op deze factor
• Quartimax (~oblimax): maximaliseert per variabele de variantie van de ladingen op de verschillende factoren
• Geomin (Mplus), oblimin, equamax, ...

Er zijn slechts vuistregels omtrent de interpretatie van de factorladingen. Vanaf .400 en zeker vanaf .600 spreken we van hoge ladingen, al hangt het telkens af van de rotatiemethode.

Links
Yahoo Answers: What are difference between varimax, quartimax and equamax rotation in factor analysis? - http://answers.yahoo.com/question/index?qid=20080616044600AAnsCc4

Cluster analysis: yeah yeah yeah

Context
For a recent working paper, I did some cluster analyses. Plural, my friend, because there's not one way. Before reading further, please understand that cluster analysis is an explorative, non-inferential method.

Problem
The Stata manual says it all: "Some researchers claim that there are as many clustering methods as there are researchers. This is untrue, there are many more methods than researchers." The wording may not be exact, but I agree with the statement. In this post, I will address some difficulties.

Issues

• Linkage method: there are five linkage methods, which define how cases get grouped: which distances to look after. There is no optimal method: single linkage causes linking patterns (one cases after another joining the same cluster), average linkage and ward's linkage are sensitive to outliers, ward's linkage distance measurement is not possible to interprete and centroid linkage will even refuse to return dendrograms.
• Similarity distance measurement: there are many distance measurements: simple, Euclidian, city block, Mahalanobis, ... Which one to choose? There's no truth.
• Cases: excluding cases may skew central values of clusters, resulting in wrongly added cases
• Order: in two-step cluster analysis, some pre-clustering is done before tackling the full data set (because it may be too large). Shuffle your deck and you'll get different results.
• Variables: ideally you cluster orthogonal factors, still, and in the general case, it doesn't mean that each variable/factor is equally important for clustering. They will all have the same a priori impact on clustering though.
• Succes garanteed: you'll always get a result. Does that mean you found something? No. 
• Stopping rules: there are many rules to determine the optimal number of clusters (see Mulligan, 1985). Calinski is default in Stata. There are some issues though, such as multiple optima or none at all. And does it make sense to choose a 23 cluster solution?

Conclusion

Cluster analysis is a wonderful way to reduce and explore data. However, I would recommend to experiment with different ways and keep my hands of it if they bear no similarity. Heck, a cluster analysis of outcomes would be useful!

Judgement: not to be trusted

Mixed Models

Andere termen: multilevelanalyse, longitudinale analyse, growth models en varianten.

Voor een eenvoudig multilevelmodel is de handleiding hier welkom.

Om statistieken te verkrijgen over het model (BIC, AIC) is er het commando postestimation
http://www.stata.com/help.cgi?xtmixed+postestimation

Een interessante toepassing is het opvatten van een mixed model als een SEM, zie:
http://blog.stata.com/2011/09/28/multilevel-random-effects-in-xtmixed-and-sem-the-long-and-wide-of-it/
(andere Stata blog posts i.v.m. SEM vind je hier).

Pearson correlatie: interpretatie en statistische intuïtie

Update 30 december 2014

"Hoe moeten we dit nu interpreteren?" - Het is een terechte vraag bij elke parameter die je voorlegt en dus zeker bij de onvermijdelijke Pearson correlatiecoëfficiënt. Deze geeft weer hoe sterk variabelen samenhangen. Het teken toont in welke richting. Zo'n simpele correlatie zal je typisch gebruiken bij een exploratieve analyse. Het is ook een belangrijke bouwsteen van allerhande regressieanalyses, factoranalyses, en SEM-modellen. Een tussenstap tussen eenvoudige correlaties en deze technieken zijn de (semi)partiële correlaties, waar je abstractie maakt van bepaalde andere variabelen die een effect mediëren.

Niemand kan echter precies zeggen hoe je dit cijfer, dat tussen -1 en 1 ligt, moet interpreteren. Het is een uitdrukking van de sterkte van een lineair verband meet. Plot de variabelen op een assenstelsel, en je zou een puntenwolk moeten zien waar zich een rechte lijn in aftekent. Kwadrateer je de correlatie, dan krijg je te zien hoeveel procent van de totale variantie door die rechte lijn wordt gevat. Dat is de technische interpretatie, maar meestal kwadrateren we niet. We hanteren in de plaats een vuistregel die sterk varieert van discipline tot discipline - geneeskundige statistiek, sociometrie, psychometrie of econometrie. In de handboeken zeilt men ook liever om de interpretatie heen. Mijn interpretatie (in sociologie & arbeidseconomie) is dit:

• r < .10: geen verband
• r < .30: zwak verband
• r < .50: matig verband
• r < .80: sterk verband
• r > .80: zeer sterk verband

Een experiment

Als alle punten van de puntenwolk exact op een rechte lijn liggen, dan is de correlatie perfect (r = 1). Hoge correlaties krijg je echter ook bij verbanden die helemaal niet zo sterk zijn, dus laat je niet misleiden. Hier een klein experiment dat je in Excel kan doen:

Maak een dataset voor 37 eenheden die elk 10 keer geobserveerd worden (dus 370 observaties in totaal, voor de snelle rekenaar). Eén maat is random veranderlijk over alle observaties, de andere is het gemiddelde per eenheid. We kunnen enige correlatie verwachten, want de twee variabelen zijn niet onafhankelijk, maar voor elke eenheid zijn er ongeveer evenveel cases onder als boven het gemiddelde, wat de correlatie laag zou moeten houden.

Welnu, de Pearson correlatie tussen de veranderlijke en haar gemiddeldes is - met de random getallen die ik verkreeg - 0.9974. Onze statistische intuïtie is duidelijk feilbaar.

Germán Rodríguez, Office of Population Research, Princeton University

http://data.princeton.edu/stata/

Een geweldige site, met informatie over General Linear Models, Mixed Models, Stata en R.

Multilevel analyse

Er is al veel te veel over multilevel analyse geschreven, vandaar een selectie van termen die voor mij nuttig zijn/waren:

Fixed vs random effects

• Fixed effects: interindividuele effecten
• Random effects: individuele effecten

Random effecten laten toe dat een bepaalde parameter voor elke eenheid verschilt, waardoor de overige parameters "onafhankelijk" van deze parameter geschat worden.

Voor intercepts is dit heel nuttig: een relatie kan voor elke L2 case gelijk zijn, maar om tal van andere redenen kan het intercept variëren. Een random effect intercept only model toont a.d.h.v. de variantie rond dit intercept (variance components), hoeveel variantie op level 2 bestaat, en dus door variabelen op dat niveau moet worden wegverklaard.

Andere parameters (slopes) kunnen ook op random gezet worden om na te gaan of er nood is aan het inbrengen van verklarende variabelen op level 2 die interageren met slope.

Bij wijze van voorbeeld deze grafiek met data voor acht fictieve landen op een fictieve Y-variabele tussen het jaar 2 en het jaar 12. De groeivoet van Y is zo geconstrueerd dat deze steeds 3% is.

Het belang zit 'm in het verschil tussen de oranje en de rode fitted line. De rode lijn volgt exact de helling van elk land. In dit geval is het gemiddelde van de trends gelijk aan de gemiddelde trend. De oranje trendlijn is gebaseerd op een steekproef uit alle punten (in het groen gemarkeerd): elk land behoudt dezelfde trend, maar hier en daar vallen er cases weg. Door het intercept te randomiseren, zou je ook deze lijn correct inschatten.

Intraclass correlatie

Betekenis: het deel van de variantie dat door groepsverschillen wordt verklaard.

Rho = sigma2(u)/som(sigma2u,sigma2res)

De sigma's worden geschat met een intercept only model (xtreg, xtmixed, anova, ... ). Vergeet niet de sigma's te kwadrateren. Voor Bernouilli-schattingen wordt een andere schatter gebruikt 

wiki: http://en.wikipedia.org/wiki/Intraclass_correlation

Three level vs cross classification

Wanneer een case tot twee niet geneste groepen behoort (vb. land en sector of vakbond en voetbalclub), spreekt men van cross classification op het tweede niveau. Termen als non hierarchical of non nested worden hier ook voor gebruikt.

Er is een kleine omweg voor de case sector en land. Als je verondersteld dat sector een landspecifiek effect heeft, dan kan je de sectoren hernoemen zodat bouw-land A verschilt van bouw-land B. Wat je dan niet onderzoekt is de algemene impact van de bouwsector, maar wel de sectorale impact binnen een land versus de internationale verschillen op level 3.

Interessant, hé.

Links

• http://faculty.chass.ncsu.edu/garson/PA765/multilevel.htm - multilevel analyse met SPSS, veel aandacht voor covariantiestructuur

Stata weights

Opnieuw iets bijgeleerd over Stata waar de gemiddelde SPSS-gebruiker niet van wakker lag: gewichten. Ik had het er al eerder over (hier). Nu blijkt dat je tóch een robuuste schatting van de standaardfouten kunt bekomen.

Even overlopen:

• pweight : doorgaans de beste keuze voor samples die herwogen zijn. Een grote waarde betekent ondervertegenwoordiging en dus een onnauwkeurige schatting. De case zal een groter belang krijgen in het berekenen van de gezochte parameters, maar de fout zal vergroten. (vb. project wage gap, decompositie)
• aweight : wanneer waarden voor een case gemiddeldes zijn: een groot gewicht betekent een meer nauwkeurige schatting. De case zal enerzijds een groter belang kijgen, en daarnaast zal de fout kleiner worden. (vb. project loonvorming)
• fweight : te gebruiken wanneer de case in feite meermaals voorkomt, maar de frequencies gecollapsed zijn (vb. datalevering Eurostat voor project Walqing)
• iweight : don't bother

Interessant, is het niet?

Links

• http://www.stata.com/help.cgi?weight
• http://www.cpc.unc.edu/research/tools/data_analysis/statatutorial/sample_surveys/weight_syntax

Threefold decomposition

Ik heb niet erg veel tijd, dus bij wijze van nota:

De decompositiemethode, met als z'n goede en kwade kanten, heeft een aantal nevenvormen. Gewoonlijk, en liefst, gebruik ik de twofold- of omegamethode. Hier bepaal je de referentiegroep voor de coëfficiënten, eventueel via gewichten.

Een andere methode is threefold, waarbij dezelfde groep de geïsoleerde bèta en gemiddelde waarde heeft. Dit heeft het voordeel dat het gemakkelijker te plotten is (zie presentatie). Het nadeel is dat er een restterm is die men CE of CxE noemt, en die interacties aangeeft. Helaas paashaas, die term is nauwelijks op een zinvolle manier te interpreteren.

Decompositie met Stata

Indertijd heb ik een macro geschreven voor SPSS, waarmee de Cotton, Reimers, Neumark & de twee Oaxaca-Blinder decomposities kunnen berekend en gebootstrapt worden. Achteraf bekeken zot werk.

Het probleem met die macro is dat je zijn regels moet volgen, en dat vraagt veel aanpassingen aan de dataset en input voor de macro. In Stata bestaan twee gelijkaardige commando's die deze klus in één lijn klaren. Het enige nadeel is dat ze elk unieke eigenschappen hebben die we helaas niet kunnen samennemen.

De commando's zijn oaxaca en nldecompose. Laat je niet afschrikken door die laatste: nl staat voor niet lineair, maar eigenlijk is het gewoon een uitbreiding van (bijna) alle regressies die Stata heeft, waaronder ook de normale regress.

Wat zijn de voor- en nadelen:

• oaxaca is erg to the point, maar lineair. Het voordeel is dat je het detail van een decompositie kunt opvragen (verklaring per (set van) verklarende variabelen), het nadeel is dat je de Cotton en Reimers specificatie niet kunt aannemen. Er is enkel: threefold, Neumark, Oaxaca 1 & 2.
• nldecompose is zoals gezegd multifunctioneel, maar heeft het nadeel dat je het detail niet kan opvragen.

Voor iemand zoals ik, die dus graag het detail van de Cotton specificatie te zien zou krijgen, is er geen oplossing. Jammer.

Gewogen data

Context
Je hebt je data gewogen, omwille van een overselectie van bepaalde kleine groepen of omwille van een onderselectie door een slecht design. Kan gebeuren.

Vraag
Wat voor analyses kan ik nog uitvoeren? Zal de weging mijn standaardfouten beïnvloeden? Zijn de coëfficiënten juist als ik niet weeg?

Antwoord
Ik heb mij slechts over twee methodes gebogen: regressies en chi²-toetsen. Mijn conclusie is dat wanneer een analyse zich op variantie baseert, wegen niet aanbevolen is. In het andere geval (chi²) wél. Even nuanceren:

Bij regressie hoef je niet te wegen: het gaat om het effect van één variabele op een andere. Weging zou eventueel outliers nadrukkelijker het effect laten bepalen, en dat kan gewenst zijn: er was immers een onderschatting van een bepaalde groep. Neem je echter de variabelen waarop de weging gebaseerd is in je model op als interactieterm, dan heb je net hetzelfde resultaat, een correcter model én juiste standaardfouten. Met andere woorden: wegen is over het algemeen af te raden.

Nadelen: misschien begrijp je niet veel meer van je coëfficiënten door alle interacties en het intercept only model geeft geen populatiegemiddelde (in het volledige model is het intercept een gecontrolleerd gemiddelde).

Chi² is een ander geval: het gaat hier om een eenvoudigere statistiek. Je vergelijkt verwachte en geobserveerde celpercentages. Die celpercentages zullen preciezer zijn bij een gestratificeerde steekproef. Maar dan moet je wegen, anders kloppen de percentages gewoonweg niet. Voor het aantal vrijheidsgraden gebruik je de ongewogen n (of de gewogen n als het gewicht een gemiddelde heeft van 1), of het aantal categorieën als de test dit vereist.

Wat t-testen betreft denk ik dat het onmogelijk is een correctie test uit te voeren: een t-test is niet anders dan een regressie waar niet voor andere variabelen gecontroleerd wordt. Dat is precies de voorwaarde om ongewogen data te mogen gebruiken. Je zou dus moeten wegen, maar dan verandert de variantie, die precies het betrouwbaarheidsinterval bepaalt. Stata heeft hiervoor speciale weegtechnieken.

Links
http://www.sociology.ohio-state.edu/ptv/faq/weights.htm
http://www.dcs.napier.ac.uk/peas/errors.htm

Gretl

http://gretl.sourceforge.net/

GRETL is een open source programma voor regressies, tijdsreeksanalyse en andere econometrische statistiek. Ik heb het gebruikt voor kwantielregressie. Best een knap programma, maar in de vorige versie zaten toch nog wat haperingen.

Squared R voor ML-analyses

Om de een of andere reden is het voor statistische software moeilijk of onmogelijk om de determinatiecoëfficiënt R² te berekenen bij analyses die niet gebaseerd zijn op een OLS-schatting.

Nochtans is dit eigenlijk niet zo moeilijk. R² is immers het kwadraat van de multiple correlatiecoëfficient R, en dit is de correlatie tussen de verwachte waarde en de geobserveerde waarde.

Zolang software je toelaat om bij een analyse een nieuwe vector met de 'predicted value' aan te maken, kun je dus zelf de Squared R berekenen.

R² = corr(pred,obs)²

Ter info: R² geeft aan welk deel van de variantie verklaard wordt door het model. Een reden waarom deze parameter enkel bij OLS-schattingen gegeven wordt is allicht dat OLS-schattingen deze verklaring maximaliseren (soms met vertekeningen van de effecten of de fout daarrond als kostprijs). Zoals ik al heb geschreven is de GENLIN proc in SPSS handiger dan REGRESSION, en geven ML-schattingen onder bepaalde assumpties exact dezelfde schattingen als OLS-schattingen, zodat bovenstaande werkwijze niet onverantwoord is. Vergeet verder niet dat R² een slechte toets is voor de kwaliteit van het model. Immers: R² zal NOOIT dalen bij wanneer extra verklarende variabelen toegevoegd worden. Nochtans peuzel je dan vrijheidsgraden op en kun je alle mogelijke problemen (vb. multicollineariteit, endogeniteit) gaan verwachten.

Vrijheidsgraden

Niemand weet precies wat vrijheidsgraden zijn, maar in de statistiek kan je er niet omheen.

De intuïtieve interpretatie is dat vrijheidsgraden het maximum zijn van de onafhankelijke rijen en kolommen van je matrix waarop je een statistiek hebt berekend, gegeven de schattingen die reeds gebeurd zijn. Bijvoorbeeld: om een variantie te berekenen heb je eerst een gemiddelde nodig. Welnu, je gemiddelde kan variëren volgens alle eenheden in de steekproef, maar eens je een gemiddelde hebt, kan je steeds de waarde van de laatste eenheid berekenen. Je hebt dus één vrijheidsgraad opgesoupeerd.

Een Chi²-test en een t-test hebben één enkele vrijheidsgraadparameter, een F-test heeft er twee (teller en noemer), omdat dit een ratio-test is.

Meer vrijheidsgraden betekent dat de test minder extreme waarden moet hebben om toch op een significante afwijking van de (nul)hypothese te wijzen.

Kookboek statistiek

De mensen van UCLA hebben een walk through voor de REGRESSION proc in SPSS (en Stata en SAS) gemaakt. Extreem handig is het hoofdstuk Regression Diagnostics. Ze leggen uit hoe je de assumpties toetst:

• Outliers en hun leverage (invloed op de coëfficiënten)
• Homoskedasticiteit
• Lineariteit: Q-Q plot, skewness, kurtosis
• Onafhankelijkheid van de waarnemingen (vs autocorrelatie): Durbin-Watson test

Correlatie- en regressiecoëfficiënten

Verband
Pearson Correlatie = bèta-coëfficiënt (bivariaat model mét intercept)

Pearson Correlatie
r = S(x,y) / S(x)*S(y)

Vraag
Waarom is een betacoëfficiënt niet gelijk bij een model met of zonder intercept?

Stationaire functies

Eng. Stationarity

Een stationaire functie is onafhankelijk van de tijd. Niet-stationaire functies zijn m.a.w. trends doorheen de tijd.

Tussen twee niet-stationaire variabelen kan een stationaire relatie bestaan. De koersen van aandelen en futures zijn bijvoorbeeld beide niet-stationaire functies, maar als één op de andere geregresseerd wordt, krijgen we een stationaire functie. Dit heet coïntegratie (Eng. cointegration). Men test dit met een Dickey-Fuller test voor de residuen, waarbij de nulhypothese is dat deze niet-stationair zijn.

Gelijkaardig aan het voorgaande, is het mogelijk om een niet-stationaire functie stationair te maken door detrending. Een eenvoudige methode is het opnemen van tijd bij de verklarende variabelen. Dit is courant bij variabelen die conjunctuurgevoelig zijn, zoals economische groei.

Heteroskedasticiteit is een vorm van niet-stationariteit van de variantie. Op wikipedia wordt het mooie voorbeeld van een cymbaalslag gegeven. Zowel het gemiddelde volume als de variantie van de frequentie (in akoestische zin) neemt af doorheen de tijd door het uitdoven en verminderen van de trillingen.

Referenties
Maritza López Novella (2001) 'Salaires conventionnels et effectifs en Belgique: une analyse empirique et macroéconomique des écarts', Working Paper 2-01, Bureau fédéral du Plan.

Context
Project Loonvorming (verklaring van de loondrift 1996-2006): we gebruiken dezelfde afhankelijke variabele (zij het niet gelogaritmeerd), maar omwille van de kleinere tijdspanne i.v.m. de studie van het Planbureau zie ik géén voortschrijdend gemiddelde. Bovendien is het sinds de invoering van de loonnorm heel aannemelijk dat er geen trendmatige vergroting van de loondrift is. Ten slotte vonden zij in hun modellen niet eens evidentie voor autoregressie. We kunnen dus concluderen dat niet-stationariteit geen issue is.