 |
 |
 |
|
{"In mathematics you don't understand
things. You just get used to them." John von Neumann.} |
Moest je net als ik, tijdens je opleiding lineaire regressie op de ingewikkelde manier leren, via taaie calculus? Dit zou kunnen verklaren waarom veel mensen de neiging hebben om een ‘bereken-nou-maar-gewoon’-houding aan te nemen bij het gebruik van R of een ander statistisch pakket. Misschien ben je echter ook wel eens in aanraking gekomen met de lineaire-algebra-benadering van lineaire regressie, die naar mijn bescheiden mening intuïtiever is. Er zijn formele artikelen hierover op het web te vinden; dit artikel is mijn bescheiden bijdrage. Wiskundegoeroes zullen terecht opmerken dat verschillende benaderingen die tot hetzelfde resultaat leiden, conceptueel toch van elkaar kunnen verschillen. Maar hoezeer ik ook van esoterische discussies houd, ik ben nu eenmaal een pragmatisch type.
Stel dat we een regressormatrix X hebben. Zie X als een matrix waarin de waarnemingen per steekproef in de kolommen zijn opgeslagen, terwijl we één rij per steekproef gebruiken. In ANOVA zou dit je ontwerpmatrix kunnen zijn. We zouden ook een reeks responswaarden hebben, georganiseerd als vector y. Het idee is dat zodra we de waarnemingen van een toekomstige steekproef kennen, we de responswaarde van deze steekproef voorspellen.
Stel dat we reden hebben om aan te nemen dat de respons op de een of andere manier lineair is, of lineair kan worden gemaakt, in de waarnemingen. In dat geval willen we een lineair model berekenen: y = X β + ε. Het model kan dan worden weergegeven door (β, ε), en daarmee kunnen we de uitkomst voor een nieuwe steekproef voorspellen door het inproduct te nemen van de waarnemingen met β en daar ε bij op te tellen.
We streven ernaar om X β = y, een exact model met ε = 0, op te lossen. Dit betekent dat we, simplistisch gezien, β = X-1 y zouden willen berekenen. Bedenk nu dat X en y doorgaans bestaan uit experimentele gegevens, behept met ruis en fouten. Dan ligt y waarschijnlijk buiten het bereik R van X, het vlak van alle exacte voorspellingen van de vorm X β, en bestaat er geen exact verband tussen X en y. Bovendien, als we minder steekproeven hebben dan waarnemingen per steekproef, is X niet vierkant en bestaat X-1 niet eens. We moeten iets slims bedenken.
Een veelgebruikte oplossing is de volgende: een matrix vermenigvuldigd met zijn over de diagonaal omgekeerde getransponeerde, is vierkant, symmetrisch en inverteerbaar. Dus vermenigvuldigen we beide zijden van de vergelijking met de getransponeerde van X, en lossen we XT X β = XT y op en vinden dat β = (XT X)-1 XT y en ε = y - X β. Je zult dit resultaat herkennen, aangezien alle afleidingen in de literatuur, hoe ingewikkeld ook, hierop uitkomen. Maar dit is nog steeds niet intuïtief – waarom is dit een kleinste-kwadraten benadering?
Merk op dat we een mogelijk niet-bestaande inverse X-1 hebben vervangen door een bestaande pseudo-inverse (XT X)-1 XT die je wellicht herkent als de Moore-Penrose inverse X+. Dan geldt: β = X+ y. Praktische lineaire regressie komt neer op het zoeken naar de beste benaderde oplossing voor een stelsel lineaire vergelijkingen dat geen exacte oplossing heeft.
Om in te zien dat de Moore-Penrose inverse een kleinste-kwadraten benadering oplevert, moeten we ons perspectief enigszins aanpassen. Om de redenering beter voorstelbaar te maken, beschouwen we het speciale geval waarin y een vector met één kolom is met N elementen en X een vierkante N×N-matrix is. Deze vereenvoudigde redenering kan worden gegeneraliseerd naar meerdere dimensies.
Stel je y niet voor als een vector, een reeks getallen, maar als een enkel punt in een N-dimensionale ruimte. Zie je y net boven R zweven? Het dichtstbijzijnde punt y' op R is de orthogonale projectie van y op R. Als je loodrecht licht op R zou schijnen, dan is y' de schaduw van y op R.
De restterm ε = y - y' is de kortste vector die elk punt op R met y verbindt, dus hebben we duidelijk een kleinste-kwadraten benadering, aangezien de som van de gekwadrateerde residuen, dat wil zeggen het kwadraat van de lengte van de vector ε, de laagst mogelijke waarde heeft.
We lossen X β = y' op in plaats van X β = y. We kunnen de onbekende y' buiten beschouwing laten, aangezien ε, loodrecht op R, in de nulruimte van de getransponeerde van X ligt. Dit betekent dat XT ε = 0. Substitutie levert XT (y - X β) = 0 op, wat opnieuw leidt tot β = X+ y en ε = y - X β.
We kunnen de onbekende y' eenvoudig berekenen met behulp van bovenstaande uitdrukkingen, namelijk als y' = X X+ y. Dit komt je misschien bekend voor, aangezien in de statistiek de matrix X X+ de projectiematrix is (nu begrijp je waarom) of de invloedsmatrix.
In een paar intuïtieve stappen zijn we tot dezelfde oplossing gekomen als die welke uit moeizame calculus is afgeleid. Het enige wat we hebben gedaan, is een punt buiten een vlak vervangen door de dichtstbijzijnde projectie op dat vlak, en basis lineaire algebra toepassen.
