Regressie
Correlatie gebruikt men om het verband tussen 2 [of meer] variabelen te beschrijven en te toetsen. We kunnen die samenhang ook gebruiken om de waarde van de ene variabele, zeg y, te voorspellen uit de andere, x. We gaan er daarbij van uit dat de relatie tussen x en y lineair is. Dan spreken we van lineaire regressie, of [lineaire] regressie-analyse.
Om y te kunnen voorspellen uit x moeten we y als een functie van x kunnen uitdrukken, bijvoorbeeld zo: y = a + bx
Van het paar x,y noemen we x de onafhankelijke variabele, of ook wel de voorspeller, en y de afhankelijke of responsie-variabele; b is de regressiecoëfficient, ofwel de richtingscoëfficient van de regressielijn, en a het intercept. De regressiecoëfficient geeft aan hoeveel verschil in y we kunnen verwachten voor iedere eenheid toename [positieve coëfficient] of afname [negatieve coëfficient] in x.
# lees dataframe in
> bron <- "http://www.mzandee.net/~zandee/statistiek/data/gegevens.txt"
> cohort <- read.table(bron, header=T)
> attach(cohort)
> names(cohort)
[1] "lichaam" "arm" "pols" "geslacht" "hand" "ogen"
> lm(arm~lichaam) # arm is de responsie-variabele
Call:
lm(formula = arm ~ lichaam)
Coefficients:
(Intercept) lichaam
0.3853 0.4128
We kunnen de waarden van de regressiecoëfficient en het intercept controleren door ze ook nog even volgens de rekenwijze in Buijs 'op de hand' uit te rekenen, door gebruik te maken van de kwadraatsommen van x en y en de produktsom van xy:
> SSX <- sum(lichaam^2)-(sum(lichaam)^2)/length(lichaam)
> SSX
[1] 6044.564
# maar het kan natuurlijk ook zo:
> SSX <- sum((lichaam - mean(lichaam))^2)
> SSY <- sum(arm^2)-(sum(arm)^2)/length(arm)
> SSY
[1] 1898.273
# of zo:
> SSY <- sum((arm - mean(arm))^2)
> SSXY <- sum(lichaam*arm)-sum(lichaam)*sum(arm)/length(arm)
> SSXY
[1] 2495.227
> b <- SSXY/SSX # de regressiecoefficient
> b
[1] 0.4128051
> sum(arm)/length(arm)-(SSXY/SSX)*sum(lichaam)/length(lichaam)
[1] 0.3853418 # het intercept a = Ey/n - b.Ex/n
We kunnen nu ook toetsen of de regressiecoëfficient b significant afwijkt van nul, dwz dat de in de steekproef gevonden samenhang, en daarmee de voorspelbaarheid, van x en y niet aan het toeval is te wijten maar op een [biologisch] effect berust. Dat toetsen doen we met behulp van de F-ratio, dwz de verhouding tussen 2 varianties. De totale variatie die we aantreffen in de twee variabelen is deels te wijten aan hun samenhang, de regressie-component, en deels aan toevalsfactoren [afleesfouten, onvolledigheden in het model, etc]. De variatiecomponent die te wijten is aan de regressie van y op x, SSR = b*SSXY , kunnen we uitrekenen, en vervolgens kunnen we via SSE = SSY - SSR ook de variatie veroorzaakt door toevalsfactoren uitrekenen. Door de regressiecomponent daarna te vergelijken met de toevalscomponent kunnen we de rol van het toeval in het tot standkomen van de gevonden regressiecoëfficient b beoordelen:
> b <- SSXY/SSX
> SSR <- b*SSXY
> SSR
[1] 1030.043
> SSE <- SSY - SSR
> SSE
[1] 868.23
|
Bron |
kwadraatsom |
vrijheidsgraden |
Gemiddelde ks |
F ratio |
|
Regressie |
1030 |
1 |
1030 |
75.9 |
|
Toevalsfouten |
868.2 |
64 |
13.57 |
|
|
Totaal |
1898.2 |
65 |
|
|
We kunnen nu uitrekenen wat de kans is om, onder de nul-hypothese dat de regressie-coëfficient nul is, een F-waarde groter dan of gelijk aan 75.9 te vinden:
> 1-pf((1030/13.5656),1,64)
[1] 1.792233e-12
Die kans is nagenoeg nul, en daarmee is de nulhypothese dus verworpen!
Bovenstaande excercitie met de variantieanalyse-tabel als resultaat kan als volgt worden samengevat in R:
> model<- lm(arm ~ lichaam)
> summary.aov(model)
Df Sum Sq Mean Sq F value Pr(>F)
lichaam 1 1030.04 1030.04 75.928 1.792e-12
Residuals 64 868.23 13.57
Nog meer informatie over het model krijgen we als volgt:
> summary(model)
Call:
lm(formula = arm ~ lichaam)
Residuals:
Min 1Q Median 3Q Max
-5.7543 -2.9519 -0.5159 2.3710 8.3738
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.38534 8.41868 0.046 0.964
lichaam 0.41281 0.04737 8.714 1.79e-12
Residual standard error: 3.683 on 64 degrees of freedom
Multiple R-Squared: 0.5426, Adjusted R-squared: 0.5355
F-statistic: 75.93 on 1 and 64 DF, p-value: 1.792e-12
Met behulp van deze tabel kunnen we, behalve over de regressiecoefficient, ook iets zeggen over het intercept a. De nulhypothese zegt dat a niet afwijkt van nul: de H0 kan niet worden verworpen want de overschrijdingskans p=0.964 is groter dan het significantieniveau alpha=0.05.
Met behulp van een strooidiagram kunnen we het verband tussen lichaamslengte en armlengte goed zichtbaar maken:
> plot(lichaam,arm,pch=16) # strooidiagram; x=lichaam, y=arm
> abline(lm(arm~lichaam)) # met regressielijn
Maar de variabele model, het resultaat van de functie lm, maakt het mogelijk om een 4-tal andere plots te maken die inzicht geven in details van de analyse (Crawley, p.144), zoals bijvoorbeeld de residuen tegen de gefitte waarden, en een normaal-kwantielplot van de residuen.
In de plot van de residuen moeten de punten random in de plot verstrooid zijn, als sterren aan een nachtelijke hemel: er mag geen structuur in de plot zichtbaar zijn, want dat zou betekenen dat het lineaire model niet goed op de data past.
Met behulp van de kwantiel plot kunnen we zien of de residuen een normaal verdeling volgen [zoals het model veronderstelt]. In deze plot moeten de punten min of meer op een rechte lijn liggen.
>plot(model)
De waarden van de gefitte data (= voorspelde waarden) en de residuen kunnen we eventueel ook zelf uitrekenen en daarna plotten, en wel als volgt:
> fit<-predict(lm(arm~lichaam)) # voorspelde waarden voor armlengtes
> residu <- arm-fit
> plot(fit, residu)
> qqnorm(residu) # normaal kwantielplot voor residuen
Bron:
Buijs, A. - Statistiek om mee te werken. Stenfert Kroese, Groningen (2003)
Crawley, M.J. - Statistics. An introduction using R. Wiley, Hoboken, NJ, USA (2005)
Slotboom, A. - Statistiek in woorden. Wolters-Noordhoff, Groningen (1987)
