--- title: "Enkelvoudige lineaire regressie" output: html_document --- # PC-les 4: Enkelvoudige lineaire regressie ## Enkelvoudige lineaire regressie met een continue predictor ### Overlevingstijd van vissen Bij 96 vissen (dojovissen, goudvissen en zebravissen) werd de resistentie tegen het gif EI-43,064 getest door elke vis individueel in een vat met 2 liter water en een bepaalde `dosis` (in mg) van het vergif te steken. Naast de overlevingstijd in minuten (de uitkomst, `minsurv`) werd ook het `gewicht` van de vis gemeten (in gram). De onderzoeksvragen zijn: "Wat is de associatie tussen dosis en overlevingstijd?" en "Is er gemiddeld gezien een verschil in overlevingstijd tussen dojovissen, goudvissen en zebravissen?" Lees de dataset poison.dat in via `read.table`. ```{r} poison <- read.table("poison.dat", sep = "", header = TRUE) # We gebruiken de volgende variabelen: soort <- poison$soort gewicht <- poison$gewicht dosis <- poison$dosis minsurv <- poison$minsurv sum(soort == 0) sum(soort == 1) sum(soort == 2) ``` Enkelvoudige lineaire regressie is een regressie waarbij de afhankelijke variabele gemodelleerd wordt in functie van slechts 1 onafhankelijke variabele, i.e. van de vorm $E[Y_i] = \beta_0 + \beta_1X_i$ (waarbij de $E$ staat voor de verwachte waarde.). Hierbij stelt $Y$ de afhankelijke variabele (of responsvariabele) en $X$ de onafhankelijke variabele (of predictorvariabele) voor. In deze opgave zullen we een antwoord op de onderzoeksvraag **"Wat is de associatie tussen dosis en overlevingstijd?"** proberen te formuleren. We gaan er voor deze oefening van uit dat, voor een gelijke dosis, het effect van het vergif op elke vis volledig gelijkaardig is (gelijke distributies van de overlevingstijden voor gelijke dosis). Merk op dat deze aanname in de praktijk waarschijnlijk niet zo realistisch is. ### 1. Maak een gepaste grafiek om de associatie tussen de dosis en de overlevingstijd te bekijken. Ga na of het realistisch is om een lineaire associatie tussen dosis en overlevingstijd te onderstellen door eerst eens de best passende kromme doorheen de puntenwolk te tekenen en vervolgens de best passende rechte. ```{r} plot(x=dosis,y=minsurv, xlab="dosis", ylab="overlevingstijd") # fit een kromme door de puntenwolk (volle lijn) lines(lowess(dosis,minsurv), lty=1) # lty=1 tekent een volle lijn # fit een rechte door de puntenwolk aan de hand van de kleinstekwadratenmethode (gestippelde lijn) abline(lsfit(dosis,minsurv), lty=2) # lty=2 tekent een stippellijn ``` De kromme benadert vrij goed een rechte. Op basis van de figuur kunnen we besluiten dat het realistisch is om een lineair verband tussen `minsurv` en `dosis` te veronderstellen. ### 2. Postuleer op basis van voorgaande grafiek een zinvol model voor de gemiddelde overlevingstijd in functie van de dosis. Model: $minsurv_i = \beta_0 + \beta_1 dosis_i + \epsilon_i$ met $\beta_0$ het (werkelijke) **intercept**, $\beta_1$ de werkelijke **helling** of meer specifiek het (werkelijk) **effect van dosis op de gemiddelde overlevingstijd** en $\epsilon_i$ een foutterm ("error term") We kunnen ook zeggen: $E[minsurv_i] = \beta_0 + \beta_1 dosis_i$ Waarbij de $E$ staat voor de verwachte waarde. ### 3. Voer een lineaire-regressieanalyse uit om de parameters in het model te schatten. ```{r} #fit een lineair regressiemodel met 'minsurv' als afhankelijke en 'dosis' als onafhankelijke variabele model1 <- lm(minsurv~dosis) model1 summaryM1 <- summary(model1) summaryM1 ``` Belangrijk! Alvorens we het lineaire regressiemodel kunnen gebruiken om geldige conclusies te trekken, moeten we de voorwaarden nagaan die moeten voldaan zijn om deze analyse te mogen uitvoeren en de resultaten te mogen vertrouwen. **De assumpties van een lineaire-regressieanalyse zijn:** 1. Onafhankelijke gegevens 2. Lineariteit tussen respons en predictor (impliceert dat residuen rond nul verdeeld zijn, zonder merkbaar resterend patroon tussen de residuen en de geschatte respons variabele) 3. Normaal verdeelde residuen 4. Homoscedasticiteit Is er hier aan de voorwaarden voldaan? 1. Onafhankelijkheid Onafhankelijkheid van de observaties wil zeggen dat kennis over de responswaarde van 1 observatie geen informatie oplevert over de responswaarde van een andere variabele. Vanwege onze (onrealistische) aanname dat de vissoort eb het gewicht, bij gelijke dosis, geen informatie geeft over hoe een vis zal reageren op het gif, kunnen we onafhankelijkheid aannemen. Immers, de onderzoekers bestudeerden in hun steekproef 96 vissen door elke vis individueel in een vat met water en een bepaalde dosis vergif te steken. Als je ervan uitgaat dat de randomisatie correct gebeurde en het selecteren van de vissen willekeurig gebeurde, kun je in principe onafhankelijkheid aannemen omdat bij een goede randomisatie observaties onafhankelijk van elkaar worden gemeten. 2. Lineariteit De lineariteitsassumptie vereist een lineair verband tussen predictorvariabele en responsvariabele in het volledige bereik van het model. Dit impliceert dat de residuen willekeurig rond nul verdeeld zijn, onafhankelijk van waar we ons op de rechte bevinden (en dus onafhankelijk van de waarde van zowel predictor als respons). Om dit na te gaan, plotten we de waarden van de residuen in functie van de door het model geschatte waarden (*fitted values*) $\widehat{minsurv}_i = \hat{\beta}_0 + \hat{\beta}_1 dosis_i$, dus de geschatte gemiddelde overlevingstijd voor elke geobserveerde dosis. We plotten ook de best passende smoother doorheen de punten. Als de assumptie van lineariteit opgaat, zullen de residuen over heel het bereik van de gefitte waarden rond 0 liggen en zal de smoother min of meer een rechte zijn en zeker geen duidelijke trends (zoals duidelijke minima of maxima) vertonen die kunnen wijzen op een kwadratisch of hogere-ordeverband dat niet in rekening werd genomen. Hier lijkt goed aan de lineariteitsassumptie voldaan te zijn, aangezien de residuen mooi rond nul liggen. De smoother ligt daarom ook dicht rond de nul-lijn zonder duidelijke trends. **Merk op:** soms kunnen 1 of enkele punten aan het uiteinde van het bereik het verloop van de smoother sterk beinvloeden, vooral als er weinig datapunten zijn. Het criterium om te bepalen of de lineariteitsassumptie geschonden is, is echter dat de residuen over heel het bereik mooi rond nul moeten liggen. De smoother is hierbij (enkel) een hulpmiddel om eventuele trends duidelijker te zien. ```{r} # Lineariteit: plot(fitted(model1),resid(model1)) # Best passende smoother: lines(lowess(fitted(model1),resid(model1))) # Best passende rechte in stippellijn: abline(h=0,lty=2) ``` 3. Normaal verdeelde residuen De residuen moeten normaal verdeeld zijn. De beste manier om dit na te gaan is aan de hand van een QQ-plot. Bij een QQ-plot worden percentielen die men heeft berekend voor de gegeven reeks observaties uitgezet t.o.v. de overeenkomstige percentielen die men verwacht op basis van de normale distributie Als de onderstelling correct is dat de gegevens normaal verdeeld zijn, dan komen beide percentielen telkens vrij goed met elkaar overeen en verwacht men bijgevolg een reeks punten min of meer op een rechte te zien. Systematische afwijkingen van een rechte wijzen op systematische afwijkingen van normaliteit. Lukrake afwijkingen van een rechte kunnen het gevolg zijn van sampling variabiliteit en/of toevallige biologische variatie en zijn daarom niet indicatief voor afwijkingen van normaliteit. ```{r} # Normaal verdeelde residuen: qqnorm(resid(model1)) qqline(resid(model1)) #lange rechtse staart + korte linkse staart. Niet echt OK ``` De residuen hebben een scheef rechtse staart en een korte linkse staart. De normaliteitsassumptie lijkt geschonden. Je kan dit ook zien op het histogram van de residuen: ```{r} hist(resid(model1)) ``` 4. Homoscedasticiteit Homoscedasticiteit betekent "gelijkheid van varianties" (heteroscedasticiteit betekent dat de varianties niet gelijk zijn). Bij lineaire regressie wil de homoscedasticiteitsassumptie zeggen dat de variantie van de residuen onafhankelijk is van waar we ons op de rechte bevinden (dus onafhankelijk van de predictor- en responsvariabele). We plotten het kwadraat van de residuen in functie van de gefitte waarden. Als we hier een smoother door trekken, zou de smoother een horizontaal verloop moeten hebben. **Waarom plotten we het kwadraat van de residuen in functie van de gefitte waarden?** De residuele variantie wordt als volgt berekend: $\frac{\sum{\hat{e}_i^2}}{n-2}$. Hier willen we kijken of de residuele variantie afhankelijk is van waar we ons op de rechte bevinden. We plotten dus de teller $\sum{\hat{e}_i^2}$ in functie van de gefitte waarden. **Waarom staat er n-2 in de noemer?** We verliezen 2 vrijheidsgraden omdat we (1) het intercept schatten en (2) het effect voor de predictor `dosis` schatten. Voor de studenten biochemie en biotechnologie en biomedische wetenschappen: wanneer we aan het hoofdstuk over multivariate lineaire regressie toekomen, zal dit $n-p$ worden, met $p$ het aantal parameters in het model. ```{r} # Homoscedasticiteit plot(fitted(model1),resid(model1)^2) lines(lowess(x=fitted(model1), y=resid(model1)^2), col="red") ``` Het nadeel van deze plot is dat het kwadraat van de residuen scheef verdeeld is, waardoor eventuele trends moeilijk te zien zijn (intuitief: als je waarden tussen -1 en 1 kwadrateert, worden ze nog kleiner, en hoe dichter bij 0, hoe kleiner ze worden. Waarden in absolute waarde groter dan 1: als je ze kwadrateert worden ze groter, en hoe groter ze zijn, hoe groter ze worden na kwadratering. Waarden die initieel standaard normaal verdeeld waren, zullen na kwadrateren dus scheef naar rechts verdeeld zijn (chi-kwadraatverdeling met 1 vrijheidsgraad).). De **vierkantswortel van de absolute waarde** van de gestandaardiseerde residuen ligt bij homoscedasticiteit ook rond nul, maar is veel minder scheef verdeeld. **Gestandaardiseerde residuen** ("internally studentized residuals") zijn een transformatie van de residuen, die een standaard normale distributie zouden vormen gegeven dat de werkelijke fouttermen een gelijke variantie hebben. Indien men dus de absolute waarde van deze residuen zou plotten in functie van de geschatte respons, dan zouden deze geen systematische afwijking van een horizontale lijn mogen volgen. Indien wel, bv. er is een systematische trend waarbij de gestandaardiseerde residuen hoger/lager worden naarmate de *fitted values* hoger/lager worden, dan betekent het dat de variantie van de residuen hoger/lager wordt naarmate de geschatte respons hoger/lager wordt. ```{r} # Gestandaardiseerde residuen standardisedResiduals <- rstandard(model1) # Plot de vierkantswortel van de absolute waarde van de gestandaardiseerde residuen in functie van de gefitte waarden plot(x=fitted(model1), y=sqrt(abs(standardisedResiduals))) # Teken een smoother door de puntenwolk lines(lowess(x=fitted(model1), y=sqrt(abs(standardisedResiduals))), col="red") ``` Uit deze plot blijkt duidelijk dat de residuele variantie toeneemt naarmate de gefitte waarden toenemen. De homoscedasticiteitsassumptie lijkt dus twijfelachtig... **Conclusie:** de residuen volgen geen normale verdeling en de assumptie van homoscedasticiteit is twijfelachtig. We kunnen eventueel beroep doen op een transformatie van de afhankelijke variabele om toch aan deze assumpties te voldoen. Dit zullen we verder behandelen in de aparte oefeningenlessen. Je kan alle assumpties tegelijk checken voor een model. Gebruik hiervoor $plot(model1)$. Merk op dat je de assumptie van homoscedasiticiteit ook kan nagaan op basis van de eerste plot (residuen vs. fitted values): als de spreiding hoger/lager wordt over de fitted values, is hieraan niet voldaan. ```{r} plot(model1) ``` ### 4. Schrijf het model neer voor de gemiddelde uitkomst met de geschatte regressiecoefficienten ingevuld. Geef de nul- en alternatieve hypothese voor de testen in de summary output. Interpreteer de geschatte regressiecoefficienten. **Opgelet: de assumpties van normaliteit en homoscedasticiteit zijn niet voldaan!** De assumpties van onafhankelijkheid en lineariteit lijken niet geschonden. Dit laatste betekent dat de schattingen van de effectgroottes correct zullen zijn. Echter: doordat de normaliteits- en/of homoscedasticiteitsassumptie geschonden zijn, zal onze besluitvorming (m.b.t. significantie) mogelijks incorrect zijn doordat onze teststatistiek niet langer een t-verdeling zal volgen (zie cursus onder "5.5. Nagaan van modelveronderstellingen"). We werken de oefening hier uit ter illustratie, maar hou dit in uw achterhoofd! ```{r} summaryM1 ``` ```{r} sigma(model1) ``` Ons model is als volgt gespecifieerd: $E[minsurv_i] = \beta_0 + \beta_1 dosis_i$ Hierbij is $E[minsurv_i]$ de verwachte (gemiddelde) overlevingstijd bij een dosis $dosis_i$. We hebben ons model zodanig geschat zodat: $minsurv_i = \hat{\beta}_0 + \hat{\beta}_1 dosis_i + e_i$ Hierbij zijn $e_i$ de residuen. We hebben de residuen zodanig gekozen dat de som van hun kwadraten zo klein mogelijk is (kleinstekwadratenschatter). Hieruit volgt dat de som van de residuen nul is. In ons geval bekomen we hetvolgende: $minsurv_i = 7.8187 - 2.2672 * dosis_i + e_i$ De geschatte gemiddelde overlevingstijd $\widehat{minsurv_i}$ gegeven een bepaalde dosis $dosis_i$ is gelijk aan: $\widehat{minsurv_i} = 7.8187 - 2.2672 * dosis_i$ #### 1. Nulhypothese, alternatieve hypothese en interpretatie van het intercept $\hat{\beta}_0$: Wat gebeurt er als $dosis_i$ = 0? Dan krijgen we dit: $\widehat{minsurv_i} = \hat{\beta}_0$ Hieraan zie je dat het intercept dient geïnterpreteerd te worden als de gemiddelde waarde van de responsvariabele gegeven dat de predictorvariabele gelijk is aan nul. We spreken steeds over de **gemiddelde waarde** van de responsvariabele omdat we de gemiddelde waarde van de responsvariabele hebben gemodelleerd in functie van de predictor. Intuïtief: herinner u de definitie van het intercept van een rechte: het is de waarde waar de rechte de y-as snijdt (dus waar x = 0). In de summary output wordt een t-test uitgevoerd om na te gaan of het werkelijke intercept verschilt van nul. De **nulhypothese** voor deze test in ons voorbeeld luidt: de gemiddelde overlevingstijd van vissen bij een dosis van 0 mg is gelijk aan 0 minuten. De **alternatieve hypothese** is: de gemiddelde overlevingstijd van vissen bij een dosis van 0 mg verschilt van 0 minuten. In ons voorbeeld luidt de correcte **interpretatie van het intercept** als volgt: **De gemiddelde overlevingstijd van vissen bij een dosis van 0 mg is gelijk aan 7,8 minuten. De p-waarde is gelijk aan 4,5e-10. Het intercept is dus extreem significant verschillend van 0 minuten op het 5%-significantieniveau.** Merk op dat deze schatting niet realistisch is: vissen die geen vergif krijgen, zouden volgens het model gemiddeld gezien na 7,8 minuten spontaan sterven. Dit komt omdat een dosis van 0 mg ver buiten het bereik van ons model ligt. In de realiteit is het veronderstelde lineair verband tussen predictor- en responsvariabele bijna altijd slechts geldig in een beperkt bereik, het bereik van de verklarende variabelen die we gemeten hebben (ook al omdat bv. veel variabelen in de praktijk niet negatief kunnen worden). We kunnen het model dus niet gebruiken om schattingen te maken buiten het bereik! Je kan dit ook zien op onderstaande figuur. Het intercept is in het rood aangeduid. Hierdoor is een nuttige biologische interpretatie van het intercept niet steeds voorhanden. ```{r} plot(x=dosis,y=minsurv, xlab="dosis", ylab="overlevingstijd", xlim=c(0,2.2)) #fit een rechte door de puntenwolk aan de hand van de kleinstekwadratenmethode (gestippelde lijn) abline(lsfit(dosis,minsurv), lty=2) #lty=2 tekent een stippellijn #duid het intercept aan (rood punt) points(0,summaryM1$coefficients["(Intercept)","Estimate"],col="red", pch=16) ``` #### 2. Nulhypothese, alternatieve hypothese en interpretatie van de predictorvariabele $\hat{\beta}_1$: Vertrek opnieuw van deze vergelijking: $\widehat{minsurv_i} = \hat{\beta}_0 + \hat{\beta}_1 dosis_i$ Wat gebeurt er als we een nieuwe $dosis_j$ definiëren die één eenheid hoger ligt? Dan krijgen we dit: $\widehat{minsurv_j} = \hat{\beta}_0 + \hat{\beta}_1 dosis_j$ Aangezien $j$ één eenheid hoger ligt dan $i$, kunnen we dit ook als volgt schrijven: $\widehat{minsurv_j} = \hat{\beta}_0 + \hat{\beta}_1 * (dosis_i+1)$ Als we nu $\widehat{minsurv_i}$ aftrekken van $\widehat{minsurv_j}$, dan bekomen we exact $\hat{\beta}_1$. Dit betekent dat $\hat{\beta}_1$ gelijk is aan de **gemiddelde toename** van de repsonsvariabele gegeven dat de predictorvariabele met één eenheid toeneemt. Intuïtief: herinner u de definitie van de richtingscoëfficiënt van een rechte: als de waarde op de x-as met één eenheid toeneemt, neemt de overeenkomstige waarde op de y-as met $\hat{\beta}_1$ eenheden toe. In de summary output wordt een t-test uitgevoerd om na te gaan of de werkelijke $\beta_1$ verschilt van nul. De **nulhypothese** voor deze test in ons voorbeeld luidt: er is geen verband tussen de dosis vergif en de gemiddelde overlevingstijd van vissen. De **alternatieve hypothese** is: er is een lineair effect van de dosis op de gemiddelde overlevingstijd van vissen. In ons voorbeeld luidt de correcte **interpretatie van de predictorvariabele $\hat{\beta}_1$** als volgt: **Vissen die 1 miligram vergif meer toegediend krijgen, leven gemiddeld gezien 2,26 minuten minder lang dan vissen die die extra miligram vergif niet toegediend kregen. Dit effect is sterk significant verschillend van 0 minuten op het 5%-significantieniveau (p = 0,0018).** Aangezien we hier met een **experimentele studie** te maken hebben, kunnen we dit ook **causaal** interpreteren binnen het bereik van het model: **Elke extra 1 miligram vergif die vissen toegediend krijgen, leidt tot een gemiddelde afname in overlevingstijd van 2,26 minuten. Dit effect is sterk significant verschillend van nul op het 5%-significantieniveau (p = 0,0018).** Merk op dat deze laatste interpretatie enkel kan bij **experimentele studies** en niet bij observationele studies. Als je bv. een lineair verbandt vindt tussen het BBP van een land en het aantal tv's per inwoner, kan je niet zeggen dat een toename van 1 tv per gezin "leidt tot" een toename in het BBP. Merk ook op dat de assumpties van normaliteit en homoscedasticiteit geschonden waren! Hierdoor volgt onze teststatistiek geen t-verdeling meer en zijn de bekomen p-waarden onnauwkeurig! ###5. Schat de gemiddelde overlevingstijd voor vissen die een dosis van 2 mg kregen en geef een bijhorend 95%-betrouwbaarheidsinterval. ```{r} predict(model1, newdata=data.frame(dosis=2.0), interval="confidence") ``` De geschatte overlevingstijd van vissen die een dosis vergif van 2 mg kregen is gelijk aan 3,28 minuten. **Interpretatie van het 95%-betrouwbaarheidsinterval:** Met een waarschijnlijkheid van 95% zal het interval 2,48 tot 4,08 de werkelijke gemiddelde overlevingstijd omvatten voor vissen die een dosis van 2 mg vergif toegediend kregen. ###6. Interpreteer de determinatiecoëfficiënt. ```{r} summary(model1) ``` **In het finale model is de determinatiecoëfficiënt gelijk aan `0.09854`. Wat betekent dit? Wat zegt dit over de kwaliteit van het model?** In het finale model is de determinatiecoëfficiënt gelijk aan 0,099. Dit betekent dat 9,9% van de variatie in de responsvariabele overlevingstijd wordt verklaard door zijn associatie met de verklarende variabele (= predictorvariabele) dosis. Het model heeft geen sterke voorspellende waarde. Dit zegt echter niets over de kwaliteit van het model! ## Enkelvoudige lineaire regressie met een categorische predictor We focussen ons nu op de tweede onderzoeksvraag: **Is er gemiddeld gezien een verschil in overlevingstijd tussen dojovissen, goudvissen en zebravissen?** We gaan ervan uit dat dojovissen gecodeerd zijn als "0", goudvissen als "1" en zebravissen als "2". Wegens beperkte tijd gaan we de assumpties voor de volgende modellen hier niet nagaan. Dit is analoog aan wat we voor `model1` getoond hebben. Dit is echter wel zeer belangrijk wanneer je zelf een lineaire regressie uitvoert! Merk op dat als er niets gegeven is, we er normaal gezien van uitgaan dat de gegevens in elke groep willekeurig verzameld zijn (dus ook willekeurig over de relevante gewichten in elke groep). Als er niets gegeven is dat op het tegendeel wijst, gaan we ervan uit dat de randomisatie goed gebeurd is en dat er aan de onafhankelijkheidsassumptie voldaan is. ###1. Fit een lineair regressiemodel voor de gemiddelde overlevingstijd in functie van de soort. ```{r} model2 <- lm(minsurv ~ soort) summary(model2) ``` **Wat is het gemiddelde verschil in overlevingstijd tussen goudvissen en dojovissen?** De gemiddelde overlevingstijd voor goudvissen ligt 0,05 minuten lager dan de gemiddelde overlevingstijd voor dojovissen. De p-waarde bedraagt 0,89. Het effect van soort op de overlevingstijd is dus niet significant verschillend van nul op het 5%-significantieniveau. **Wat is het gemiddelde verschil in overlevingstijd tussen zebravissen en goudvissen?** De gemiddelde overlevingstijd voor zebravissen ligt 0,05 minuten lager dan de gemiddelde overlevingstijd voor goudvissen. De p-waarde bedraagt 0,89. Het effect van soort op de overlevingstijd is dus niet significant verschillend van nul op het 5%-significantieniveau. **Wat is het gemiddelde verschil in overlevingstijd tussen zebravissen en dojovissen?** De gemiddelde overlevingstijd voor zebravissen ligt 0,04698*2 = 0,09 minuten lager dan de gemiddelde overlevingstijd voor dojovissen. De p-waarde bedraagt 0,89. Het effect van soort op de overlevingstijd is dus niet significant verschillend van nul op het 5%-significantieniveau. **Is deze codering zinvol? Waarom wel/niet?** Deze codering is weinig zinvol, dit model onderstelt impliciet dat het verschil in responsvariabele dosis 2 keer groter is tussen zebravissen en dojovissen dan tussen goudvissen en dojovissen, omdat het de variabele `soort` aanziet als een numerieke variabele. We zien dit in door het model grafisch voor te stellen: ```{r} plot(minsurv~soort) abline(lm(minsurv~soort),col=2) ``` ###2. Maak nu eerst 2 dummyvariabelen. Fit dit model in R en interpreteer de coefficienten. Oplossing: dummy-variabelen: - variabele $soortgoud$ is 1 als de observatie overeenkomt met een goudvis en nul anders. - variabele $soortzebra$ is 1 als de observatie overeenkomt met een zebravis en nul anders. We kunnen nu afstappen van de onderstelling dat de overlevingstijd lineair is in het soort vis door een model te fitten met als verklarende variabelen deze twee nieuwe variabelen. Merk op dat we met telkens werken met *k*-1 dummy-variabelen waarbij *k* het aantal niveaus is van de oorspronkelijke factor-variabele. Het resterende niveau wordt dan gerepresenteerd door het intercept. Elk niveau van de categorische variabele is eenduidig gedefinieerd door een combinatie van de dummies: **Tabel 1.** Waarden van de dummyvariabelen "soortgoud" en "soortzebra" voor elk level van de factorvariabele "soort". Soort soortgoud soortzebra --------- --------- ---------- dojovis 0 0 goudvis 1 0 zebravis 0 1 ```{r} soortgoud <- ifelse(soort==1,1,0) # soortgoud is 1 als soort==1 en anders 0 soortzebra <- ifelse(soort==2,1,0) # soortzebra is 1 als soort==2 en anders 0 model3 <- lm(minsurv ~ soortgoud + soortzebra) summary(model3) ``` Bovenstaande aanpak kan vrij inefficiënt zijn als je veel levels hebt (veel vissoorten in dit voorbeeld), omdat je zelf voor ieder level een nieuwe variabele moet aanmaken. Een efficiëntere manier is om soort als een (discrete!) **factorvariabele** te definiëren en dit mee te nemen in het model. Je vertelt R op deze manier dat soort een factorvariabele is (en geen numerieke variabele). Wanneer je je regressiemodel fit, zal R automatisch een dummycodering gebruiken om om je model te fitten. Vergelijk de parameterschattingen! ```{r} fSoort <- factor(soort) fmodel <- lm(minsurv ~ fSoort) summary(fmodel) ``` We zien dat het model nu drie parameters heeft geschat; waarin de gemiddelde overlevingstijd voor de drie vissoorte in vervat zitten. ```{r} beta=coef(fmodel) # geschatte regressiecoefficienten plot(minsurv~soort) points(x=0, y=beta[1], pch=17, cex=3/2, col="red") # intercept: dojovis points(x=1, y=beta[1]+beta[2], pch=17, cex=3/2, col="red") # intercept + fsoort1: goudvis points(x=2, y=beta[1]+beta[3], pch=17, cex=3/2, col="red") # intercept + fsoort2: zebravis legend("topright",c("data","geschatte gemiddelde"),pch=c(1,17), col=1:2, bty="n", cex=2/3) ``` #### 1. Interpretatie van het intercept Zoals voorheen moet het intercept geinterpreteerd worden als de gemiddelde respons wanneer de predictoren gelijk zijn aan nul. Uit tabel 1 blijkt dat dit overeenkomt met dojovissen. Het intercept moet dus als volgt geinterpreteerd worden: **De gemiddelde overlevingstijd van dojovissen is gelijk aan 3,4 minuten. Dit effect is extreem significant verschillend van nul (p < 2e-16).** #### 2. Interpretatie van "fSoort1" (of "soortgoud") We weten dat dojovissen voorgesteld worden door het intercept, en dat `soort==1` voor goudvissen. De parameter die `R` `fSoort1` noemt, slaat dus op de goudvissen. `R` zal in de output voor een factor variabele in een lineair model namelijk steeds de naam van de variabele gebruiken gevolgd door het niveau dat de parameter voorstelt. Je kon dit ook bekijken via de tabel hierboven: Als we dojovissen met goudvissen vergelijken, zien we dat de enige variabele die verandert `soortgoud` is. Deze variabele moet dus als volgt geinterpreteerd worden: **De gemiddelde overlevingstijd voor goudvissen ligt 2,7 minuten hoger dan de gemiddelde overlevingstijd voor dojovissen. Dit effect is extreem significant op het 5%-significantieniveau (p = 4,2e-08).** Merk op dat we deze parameter enkel kunnen interpreteren ten opzichte van de referentieklasse, namelijk het intercept. #### 3. Interpretatie van "fSoort2" (of "soortzebra") Analoog als hierboven weten we dat dojovissen voorgesteld worden door het intercept, en dat `soort==2` voor zebravissen. Deze variabele moet dus als volgt geïnterpreteerd worden: **De gemiddelde overlevingstijd voor zebravissen ligt 1,0 minuten lager dan de gemiddelde overlevingstijd voor dojovissen. Dit effect is niet significant maar echter wel suggestief op het 5%-significantieniveau (p = 0,06).** ####Wat als we de gemiddelde overlevingstijd van zebravissen willen vergelijken met goudvissen? Er is geen parameter in het model die het verschil in gemiddelde overlevingstijd tussen zebravissen en goudvissen aangeeft. We kunnen uiteraard wel de grootte van dit effect berekenen via: -1,0452 - 2,7117 = -3,8. De gemiddelde overlevingstijd voor zebravissen ligt dus 3,8 minuten lager dan de gemiddelde overlevingstijd voor goudvissen. Om echter statistische inferentie te kunnen doen op dit effect, moeten we een combinatie van modelparameters testen, wat volgende week aan bod komt in het practicum ANOVA. Een andere oplossing om wel inferentie te kunnen doen op dit effect, is je dummyvariabelen anders te coderen, waarbij we "goudvis" als referentieklasse nemen. Dit kan als volgt: ```{r} fSoortB <- relevel(fSoort, ref="1") # 1 staat voor "goudvis", dit wordt nu de nieuwe referentieklasse fmodel2 <- lm(minsurv ~ fSoortB) summary(fmodel2) ``` `fSoortB0` moet nu geinterpreteerd worden als het verschil in gemiddelde overlevingstijd van dojovissen t.o.v. goudvissen. `fSoortB2` is dan het verschil in gemiddelde overlevingstijd van zebravissen t.o.v. goudvissen. Dit laatste effect is uiteraard exact gelijk aan de -3,8 die we op basis van `fmodel` berekenden. **Merk op:** We willen hier eigenlijk drie nulhypothesen testen (verschil in gemiddelde overlevingstijd tussen goudvissen en dojovissen, tussen zebravissen en dojovissen, en tussen zebravissen en goudvissen). Hoewel elke test hier apart op het 5\%-significantieniveau gecontroleerd is, is de kans om minstens één nulhypothese onterecht te verwerpen, gegeven dat er geen verschil is tussen de drie vissoorten, groter dan 5\%. Om hieraan tegemoet te komen, verwijzen we opnieuw naar het volgende practicum: ANOVA.