--- title: "Oefensessie 3: Enkelvoudige lineaire regressie - full" output: html_document: code_download: yes theme: cosmo toc: yes toc_float: yes highlight: tango number_sections: yes pdf_document: toc: yes --- 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?" # Laad de libraries ```{r} library(dplyr) library(ggplot2) library(readr) ``` # Importeer de data Lees de dataset "poison.dat" met de `read_delim` functie. ```{r} poison <- read_delim("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum4/poison.dat", col_types = list(col_factor(), col_double(), col_double(),col_double()), col_names = TRUE, delim = " ") head(poison) ``` Zoals aangegeven in de introducerende tekst werd het experiment uitgevoerd op drie verschillende vissoorten: ```{r} poison %>% pull(soort) %>% table ``` Het doel van deze analyse is het bestuderen van de associatie tussen dosis en overlevingstijd in deze 96 vissen. Echter, de verschillende observaties in deze dataset zijn niet onafhankelijk van elkaar. We verwachten dat vissen van eenzelfde soort meer gelijkaardig zullen reageren op een bepaalde dosis gif dan vissen van een andere soort. Wanneer we deze dataset analyseren zonder hiermee rekening te houden schenden we de assumptie van onafhankelijke observaties (zie verder voor de assumpties van het lineaire model). Dit kunnen we op twee manieren "oplossen". 1. Een lineair model opstellen dat expliciet in rekening brengt tot welke soort een bepaalde vis behoort. Deze strategie zal aan bod komen in het hoofdstuk meervoudige lineaire regressie, *enkel voor de studenten biochemie en biomedische wetenschappen*. 2. Een subset van de data nemen zodat er nog maar 1 soort overblijft. Aangezien meervoudige lineaire regressie nodig is voor optie 1, zullen we ons in deze oefening beperken tot het analyseren van de data van soort "0" (dojovissen). ```{r} poison <- poison %>% filter(soort=="0") # filter de data zodat enkel soort "0" overblijft ``` In deze opgave zullen we een antwoord op de onderzoeksvraag **"Is er een lineaire associatie tussen dosis en overlevingstijd?"** proberen te formuleren. # Data exploratie We gaan 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} poison %>% ggplot(aes(x=dosis,y=minsurv)) + geom_point() + stat_smooth(method = "loess",col="red") + # fit een kromme door de puntenwolk (rode lijn) stat_smooth(method='lm',col="black") + # fit een rechte door de puntenwolk aan de hand van de kleinstekwadratenmethode ylab("Overlevingstijd (min)") + xlab("Dosis (mg)") + theme_bw() ``` 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. # Enkelvoudige lineaire regressie analyse ## Formulering van een lineair model voor de gemiddelde overlevingstijd in functie van de dosis. Enkelvoudige lineaire regressie is een regressie waarbij de afhankelijke variabele gemodelleerd wordt in functie van 1 onafhankelijke variabele, i.e. van de vorm $E[Y_i] = \beta_0 + \beta_1X_i$ (waarbij de $E$ staat voor de verwachtingswaarde.). Hierbij stelt $Y$ de afhankelijke variabele (of responsvariabele) en $X$ de onafhankelijke variabele (of predictorvariabele) voor. In dit voorbeeld krijgen we dus het volgende model: $$minsurv_i = \beta_0 + \beta_1 dosis_i + \epsilon_i$$ - $\beta_0$ is het (werkelijke) **intercept**, - $\beta_1$ is de werkelijke **helling** of meer specifiek het (werkelijk) **effect van dosis op de gemiddelde overlevingstijd** - $\epsilon_i$ is een foutterm ("error term"), waarbij $\epsilon_i$ i.i.d. normaal verdeeld zijn met gemiddelde 0 en (constante) variantie $\sigma^2$. We kunnen ook schrijven: $$E[minsurv_i] = \beta_0 + \beta_1 dosis_i$$ waarbij $E$ opnieuw staat voor de verwachtingswaarde. ## 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,data = poison) model1 summary(model1) ``` ## Assumpties 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. Bekijk hiervoor zeker nog eens de kennisclip "6.5 Nagaan van modelveronderstellingen". **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? ### 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. Door de analyse slechts op 1 vissoort uit te voeren, hebben we het probleem dat responswaarden binnen eenzelfde vissoort meer gelijkaardig zijn dan responswaarden tussen vissoorten (zie onderstaande figuur) reeds "opgelost". Er is nog echter nog een tweede predictor die een rol kan spelen; het is te verwachten dat vissen met een hoger gewicht beter bestand zullen zijn tegen het gif. Wanneer er in het experiment een systematische bias zou zijn, bijvoorbeeld wanneer alle lichtere vissen en hogere of lagere dosis van het gif kregen, dan kunnen we het lineaire verband tussen dosis en overlevingstijd niet juist schatten, aangezien er dan confounding is met gewicht. ```{r} poison %>% ggplot(aes(x = dosis, y = gewicht)) + geom_point() + ylab("Gewicht (g)") + xlab("Dosis (mg)") + theme_bw() ``` We kunnen er echter vanuit gaan dat de vissen volledig willekeurig werden toegewezen aan een bepaalde dosis. Visueel zien we geen systematische trend tussen het gewicht van de vissen en de dosis gif die ze kregen. Wel merken we dat gewichten niet volledig uniform verdeeld zijn over de dosissen. Als we de vissen at random toekennen aan de verschillende dosissen kan dergelijke variabiliteit wel door toeval voorkomen. Wanneer men dat wil verhinderen had met de vissen via gebalanceerde randomisatie kunnen toewijzen waarbij men de gewichten zou balanceren per dosis. Het design van het experiment is dus een belangrijke factor bij het nagaan van de onafhankelijkheidsvoorwaarde. ### 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. Deze plot wordt by default weergegeven door de plot functie voor een lm object ```{r} # Lineariteit: plot(model1, which = 1) ``` 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. ### 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. Voor meer uitleg rond QQ-plots verwijzen we naar de kennisclip "4.4 De Normale benadering van gegevens". Ook deze plot wordt by default gegeneerd door de plot functie van een lm object. ```{r} # Normaal verdeelde residuen: plot(model1, which = 2) ``` De residuen hebben een scheef rechtse staart en een korte linkse staart. De normaliteitsassumptie lijkt geschonden. ### 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 hiervoor de **vierkantswortel van de absolute waarde** van de gestandaardiseerde residuen in functie van de geschatte respons. **Gestandaardiseerde residuen** 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 plot(model1, which = 3) ``` Deze plot suggereert 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 later in deze oefening. 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, ,fig.show="hold", out.width="50%"} plot(model1) ``` ## Model voor de gemiddelde uitkomst met de geschatte regressiecoefficienten ingevuld. **Opgelet: de assumpties van normaliteit en homoscedasticiteit zijn hier geschonden!** 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 zodat u ook de interpretatie ziet voor een lineair model waarbij geen transformatie is uitgevoerd. Maar hou steeds in uw achterhoofd dat p-waarden en betrouwbaarheidsintervallen niet correct zijn! ```{r} summary(model1) ``` Ons model is als volgt gespecifiëerd: $$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 het volgende: $$minsurv_i = 6.9572 - 2.3029 \times dosis_i + e_i$$ De geschatte gemiddelde overlevingstijd $\widehat{minsurv_i}$ gegeven een bepaalde dosis $dosis_i$ is dus gelijk aan: $$\widehat{minsurv_i} = 6.9572 - 2.3029 \times dosis_i$$ ## 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. **De interpretatie van het intercept $\hat{\beta}_0$:** *De gemiddelde overlevingstijd van vissen bij een dosis van 0 mg is gelijk aan 6.95 minuten.* Merk op dat deze interpretatie **niet wetenschappelijk relevant** is: vissen die geen gif krijgen, zouden volgens het model gemiddeld gezien na 6.9 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 de interpretatie van het intercept dus dikwijls wetenschappelijk irrelevant. We zullen dus ook meestal geen formele statistische test uitvoeren voor het intercept. ```{r} poison %>% ggplot(aes(x=dosis,y=minsurv)) + geom_point() + xlim(0,2.1) + stat_smooth(method='lm',col="black",fullrange=TRUE,lty=2) + # fit een rechte door de puntenwolk aan de hand van de kleinstekwadratenmethode (gestippelde lijn) geom_point(aes(x=0, y=model1$coefficients["(Intercept)"]), colour="red",size=3) + ylab("Overlevingstijd (min)") + xlab("Dosis (mg)") + theme_bw() ``` We zullen dus ook meestal geen formele statistische test uitvoeren voor het intercept. Ter illustratie doen we dit hier toch, maar merk op dat dit wetenschappelijk niet relevant is. In de summary output wordt een t-test uitgevoerd om na te gaan of de werkelijke $\beta_0$ verschilt van nul. - De **nulhypothese** luidt: *het intercept is gelijk aan nul* (of in symbolen: $\beta_0=0$ ) - De **alternatieve hypothese** is: *het intercept is verschillend van nul* (of in symbolen: $\beta_0\neq0$) We zien dat de p-waarde zeer klein is (p = 2.32e-08) en we kunnen dus besluiten dat het intercept extreem significant verschilt van 0. ## Nulhypothese, alternatieve hypothese en interpretatie voor het effect van dosis ($\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 ($dosis_j=dosis_i+1$)? Dan krijgen we dit: \begin{aligned} \widehat{minsurv_j} &= \hat{\beta}_0 + \hat{\beta}_1 dosis_j\\ &= \hat{\beta}_0 + \hat{\beta}_1 (dosis_i+1)\\ &= \hat{\beta}_0 + \hat{\beta}_1 dosis_i+\hat{\beta}_1)\\ &= \widehat{minsurv_i} + \hat{\beta}_1 \end{aligned} Dit betekent dat $\hat{\beta}_1$ gelijk is aan de **gemiddelde toename** van de responsvariabele 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 gif en de gemiddelde overlevingstijd van vissen* (of in symbolen: $\beta_1=0$). - De **alternatieve hypothese** is: *er is een lineair effect van de dosis op de gemiddelde overlevingstijd van vissen* (of in symbolen: $\beta_1\neq 0$). We zien dat de p-waarde zeer klein is (p = 0.000767) en we kunnen dus besluiten dat het effect extreem significant verschilt van 0. **De interpretatie van de hellingsparameter $\hat{\beta}_1$** *Vissen die 1 milligram gif meer toegediend krijgen, leven gemiddeld gezien 2,30 minuten minder lang dan vissen die die extra milligram gif niet toegediend kregen. Dit effect is extreem significant verschillend van 0 minuten op het 5%-significantieniveau (p = 0.000767).* Omdat deze studie een **experimentele studie** is, kunnen we de associatie tussen het gif en de overlevingstijd **causaal** interpreteren: *De overlevingstijd van dojo-vissen neemt gemiddeld gezien met 2,30 minuten af per gram gif die wordt toegediend. Dit effect is extreem significant verschillend van nul op het 5%-significantieniveau (p = 0.000767).* Merk op dat deze laatste interpretatie enkel kan bij **experimentele studies** en niet bij observationele studies. Als je bv. een lineair verband 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! ## Geef en interpreteer de betrouwbaarheidsintervallen voor het intercept en de hellingsparameter. De betrouwbaarheidsintervallen vinden we met: ```{r} confint(model1) ``` **De interpretatie van het 95%-betrouwbaarheidsinterval voor intercept:** *Het interval [`r confint(model1)[1,] %>% round(.,2) %>% sort`]) bevat met 95% waarschijnlijkheid het werkelijke intercept.* **De interpretatie van het 95%-betrouwbaarheidsinterval voor het effect van `dosis`:** *Het interval [`r -confint(model1)[2,] %>% round(.,2) %>% sort`]) bevat met 95% waarschijnlijkheid het werkelijke gemiddelde verschil in overlevingstijd in minuten tussen dojovissen die 1 mg gif minder krijgen dan vissen die meer gif krijgen.* ## 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 gif van 2 mg kregen is gelijk aan 2.35 minuten. **Interpretatie van het 95%-betrouwbaarheidsinterval:** *Het interval 1.60 min tot 3.10 min omvat met 95% waarschijnlijkheid de werkelijke gemiddelde overlevingstijd voor vissen die een dosis van 2 mg gif toegediend kregen.* ## Interpreteer de determinatiecoëfficiënt. ```{r} summary(model1) ``` In het finale model is de determinatiecoëfficiënt gelijk aan 0.2665. Dit betekent dat 26.65% van de variatie in de responsvariabele overlevingstijd wordt verklaard door zijn associatie met de verklarende variabele (= predictorvariabele) dosis. Het model zal hierdoor niet bruikbaar voor predictiedoeleinden. De determinatiecoëfficiënt zegt echter niets over de kwaliteit van het model! **Merk op dat de voorwaarden niet zijn voldaan voor dit model, we kunnen de output van de statistische testen en betrouwbaarheidsintervallen dus niet vertrouwen. We zullen daarom nagegaan of de aannames voldaan zijn voor het model na transformatie van de uitkomst.** # Analyse na log-transfromatie ## Waarom een log-transformatie? - Transformeren van de uitkomst kan zinvol zijn als er niet is voldaan aan normaliteit of gelijkheid van variantie. - Als niet is voldaan aan lineariteit kan men nagaan of het zinvol is om de predictor te transformeren. - Herinner u: dat de verdeling van de residuen een korte linkse staart en een lange rechtse staart had en dat er een stijgende trend leek van de variabiliteit bij toenemend gemiddelde en dat de uitkomst, de overlevingstijd steeds positief is. Voor een positieve uitkomst kunnen dergelijke afwijkingen kunnen vaak geremedieerd worden d.m.v. log-transformatie van de response variabele. ```{r} plot(model1, which = 2) ``` Een log-transformatie maakt grote waarden heel veel kleiner terwijl kleine waarden slechts een beetje kleiner zullen worden. Een log-transformatie heeft dus het potentieel om verdelingen die scheef naar rechts zijn meer symmetrisch te maken. Herinner u ook de plot van de vierkantswortel van de absolute waarde van de residuen i.f.v. de gefitte waarden. Hieruit besloten we dat de variantie toeneemt naarmate de geschatte uitkomstvariabele groter wordt. ```{r} plot(model1, which = 3) ``` Dit is een extra indicatie dat een log-transformatie zinvol kan zijn: een logtransformatie zal een additieve errorstructuur omzetten naar een multiplicatieve structuur. ## Voer een log-transformatie uit van de overlevingstijd ```{r} poison <- poison %>% mutate(logminsurv = log10(minsurv)) ``` Merk op dat de functie `log` in `R` het natuurlijk logaritme (d.w.z. het logaritme met grondtal $e$) berekent. Om een logaritme met een ander grondtal te berekenen, kan je het argument `base = ...` gebruiken. De volgende code berekent het logaritme van 100 met grondtal 10. Voor het logaritme met basis 2 en 10 bestaat ook de `log2` en `log10` functie. ```{r} log(100, base=10) log10(100) ``` ```{r} #fit een lineair regressiemodel met het 10-delig logartime van minsurv als afhankelijke en 'dosis' als onafhankelijke variabele logModel <- lm(logminsurv~dosis,data=poison) logModel ``` ## Assumpties ```{r} plot(logModel, which=1:3) ``` ### Onafhankelijkheid Idem als voordien: als je ervan uitgaat dat de randomisatie correct gebeurde en het toekennen van de vissen aan de concentratie willekeurig gebeurde, kun je in principe onafhankelijkheid aannemen omdat bij een goede randomisatie observaties onafhankelijk van elkaar worden gemeten. ### Lineariteit Lineariteit gaan we na op basis van de "Residuals vs Fitted" plot. Deze toont de waarden van de residuen in functie van de door het model gefitte waarden. Hier lijkt goed aan de lineariteitsassumptie voldaan te zijn, aangezien de residuen mooi rond nul verdeeld zijn over heel het bereik van de gefitte waarden. De smoother ligt daarom ook dicht rond de nul-lijn zonder duidelijke trends. ### Normaal verdeelde residuen Hiervoor kijken we naar de QQ-plot. De residuen lijken niet sterk af te wijken van wat men zou verwachten als ze normaal verdeeld zouden zijn (de percentielen van de residuen komen goed overeen met de percentielen die men verwacht op basis van de normale verdeling). De iets kortere linkse staart is niet van die aard dat we ons zorgen zouden moeten maken over schendingen van normaliteit. De normaliteitsassumptie lijkt dus voldaan. ### Homoscedasticiteit Hiervoor kijken we naar de plot die de gefitte waarden uitzet in functie van de vierkantswortel van de absolute waarde van de gestandaardiseerde residuen. Als de data homoscedastisch is, dan zal het gemiddelde van de absolute waarde van de gestandaardiseerde residuen altijd rond dezelfde waarde liggen, onafhankelijk van de gefitte waarde. Een smoother door de puntenwolk zal dan vrijwel horizontaal zijn en geen duidelijke trends vertonen. Er is geen patroon te zien wanneer de vierkantswortel van de gestandaardiseerde residuen wordt uitgezet in functie van de gefitte waarden. Onze smoother loopt dan ook vrijwel horizontaal, de zeer beperkte stijging is niet van die aard dat we ons zorgen moeten maken over de homoscedasticiteitsassumptie. De assumptie van homoscedasticiteit lijkt dus voldaan. In de plot van de residuen in functie van de gefitte waarden zien we dat de residuen mooi gelijk gespreid liggen rond 0 en dat de spreiding gelijk blijft in functie van de gefitte uitkomst. **Besluit**: alle assumpties lijken voldaan te zijn. We gaan dus een correcte inferentie kunnen doen op basis van dit model. ## Schrijf het model neer voor de gemiddelde log-getransformeerde overlevingstijd met de geschatte regressiecoefficienten ingevuld, zowel op log-schaal als op originele schaal ```{r} summary(logModel) ``` ### Model op log-schaal Ons model is als volgt gespecifiëerd: $$E[log(minsurv_i)] = \beta_0 + \beta_1 dosis_i$$ Hierbij is $E[log(minsurv_i)]$ de verwachtingswaarde of het (rekenkundig) gemiddelde is van de log10-overlevingstijd bij een dosis $dosis_i$. De geschatte gemiddelde log10-getransformeerde overlevingstijd, gegeven een bepaalde dosis $dosis_i$ is gelijk aan: $$\widehat{\log_{10}(minsurv_i)} = `r round(logModel$coef[1],3)` - `r round(logModel$coef[2],3) %>% abs` \times dosis_i$$ ### Model op originele schaal Om het model te geven op originele schaal, transformeren we de vergelijking $E[log(minsurv_i)] = \beta_0 + \beta_1 dosis_i$ met de inverse transformatie, in dit geval dus door links en rechts de macht van 10 te nemen \begin{aligned} 10^{E[log10(minsurv_i)]} &= 10^{\beta_0 + \beta_1 dosis_i}\\ &= 10^{\beta_0} 10^{\beta_1 dosis_i}\\ &= 10^{\beta_0} 10^{\beta_1 dosis_i}\\ &= 10^{\beta_0} \left(10^{\beta_1}\right)^{dosis_i} \end{aligned} Uit de theorie weten we dat het linkerlid staat voor het geometrisch gemiddelde van de overlevingstijd, wat we zullen noteren met GM(minsurv_i). Stellen we nu $B_0=10^{\beta_0}$ en $B_1=10^{\beta_1}$, dan krijgen we $$GM(minsurv_i)=B_0\, B_1^{dosis_i} $$ We berekenen nu $\hat{B_0}$ en $\hat{B_1}$ ```{r} 10^logModel$coefficients #B0 en B1 ``` Het geschatte geometrisch gemiddelde van de overlevingstijd, gegeven een bepaalde dosis $dosis_i$ is dus gelijk aan: $$\widehat{GM(minsurv_i)} = 8.26 * 0.53^{dosis_i}$$ ## Nulhypothese, alternatieve hypothese en interpretatie van het intercept $\hat{\beta}_0$, inclusief 95%-betrouwbaarheidsinterval, zowel op de log-schaal als op de originele schaal ### Nulhypothese en alternatieve hypothese - De **nulhypothese** is: *het gemiddeld 10-delig logaritme van de overlevingstijd in minuten voor dojovissen die geen gif toegediend kregen is gelijk aan 0* (of in symbolen: $\beta_0=0$). - De **alternatieve hypothese** is: *het gemiddeld 10-delig logaritme van de overlevingstijd in minuten voor dojovissen die geen gif toegediend kregen is niet gelijk aan 0* (of in symbolen: $\beta_1=0$). ### Interpretatie van het intercept op log-schaal ```{r} summary(logModel) ``` ```{r} logModel$coefficients confint(logModel) ``` **Interpretatie van het intercept op log-schaal** *Het gemiddeld 10-delig logaritme van de overlevingstijd in minuten voor dojovissen die geen gif toegediend kregen is gelijk aan `r logModel$coef[1] %>% round(.,2)` (95%-betrouwbaarheidsinterval: [`r confint(logModel)[1,] %>%round(.,2)`]).* **Interpretatie van het 95%-betrouwbaarheidsinterval voor het intercept op de log-schaal:** *Het interval [`r confint(logModel)[1,] %>%round(.,2)`] bevat met 95% waarschijnlijkheid het werkelijke gemiddelde 10-delige logaritme van de overlevingstijd in minuten voor vissen die geen gif toegediend kregen.* ### Interpretatie van het intercept op de originele schaal ```{r} 10^logModel$coefficients #schattingen voor B0=10^beta0 en B1=10^beta1 10^confint(logModel) #Betrouwbaarheidsintervallen voor B0 en B1 ``` Ons geschatte model was $$\widehat{GM(minsurv_i)} = \hat{B_0}\hat{B_1}^{dosis_i}=8.26 \times 0.53^{dosis_i}$$ waar $\hat{B_0}=10^{\hat{\beta_0}}$ en $\hat{B_1}=10^{\hat{\beta_1}}$ Als we de $dosis_i$ gelijk stellen aan nul dan krijgen we $$\widehat{GM(minsurv_i)} = \hat{B_0}= 8.26$$ **De interpretatie van het intercept op de originele schaal wordt dus:** *Het geometrisch gemiddelde van de overlevingstijd voor dojovissen die geen gif toegediend kregen is gelijk aan `r 10^logModel$coef[1] %>% round(.,1)` minuten (95%-betrouwbaarheidsinterval:[`r 10^(confint(logModel)[1,]) %>%round(.,1)`]). * **Interpretatie van het 95%-betrouwbaarheidsinterval voor het intercept op de originele schaal:** *Het interval [`r 10^(confint(logModel)[1,]) %>%round(.,1)`] omvat met 95% waarschijnlijkheid het werkelijke geometrische gemiddelde van de overlevingstijd in minuten voor vissen die geen gif kregen toegediend.* Merk nogmaals op dat de interpretatie van het intercept wetenschappelijk niet relevant is. Het berust op een sterke extrapolatie en het schatten van de overlevingstijd van vissen die geen gif kregen op basis van de data in de steekproef is dus niet relevant. Daarom hebben we ook de graad van significantie achterwege gelaten bij de interpretatie van het intercept. ```{r} range(poison$dosis) #print de minimale en maximale waarde van de dosis-variabele ``` ## Nulhypothese, alternatieve hypothese en interpretatie van het effect van de dosis op de overlevingstijd, inclusief 95%-betrouwbaarheidsinterval, zowel op de log-schaal als op de originele schaal. ### Nulhypothese en alternatieve hypothese - **Nulhypothese voor het dosiseffect:** *er is geen verband tussen de dosis gif die de dojovissen toegediend kregen en het gemiddeld 10-delig logaritme van hun overlevingstijd* (of in symbolen: $\beta_1=0$). - **Alternatieve hypothese voor het dosiseffect:** *er is een lineair effect van de dosis gif op het gemiddeld 10-delig logaritme van hun overlevingstijd in minuten van dojovissen* (of in symbolen: $\beta_1\neq 0$). ### Interpretatie op de log-schaal ```{r} summary(logModel) ``` ```{r} logModel$coefficients #schattingen voor beta0 en beta1 confint(logModel) #betrouwbaarheidsintervallen voor beta0 en beta1 ``` **De interpretatie van de het effect van dosis op de log10-overlevingstijd is dus:** *Voor dojovissen die aan een hogere concentratie van het gif worden blootgesteld is het 10-delig logaritme van de overlevingstijd in minuten gemiddeld `r logModel$coef[2] %>% abs %>% round(.,2)` per mg gif lager dan dat van vissen die aan een lagere dosis worden blootgesteld (95%-betrouwbaarheidsinterval: [`r confint(logModel)[2,] %>% round(.,2)`]). Dit verschil is extreem significant (p < 0.001).* **De interpretatie van het 95%-betrouwbaarheidsinterval voor het effect van dosis op de log10-overlevingstijd:** *Het interval [`r -confint(logModel)[2,] %>% round(.,2) %>% sort`]) bevat met 95% waarschijnlijkheid het werkelijke gemiddelde verschil in 10-delig logaritme van de overlevingstijd in minuten tussen dojovissen die 1 mg gif minder krijgen dan vissen die meer gif krijgen.* ### Interpretatie op de originele schaal ```{r} 10^logModel$coefficients #schattingen voor B0=10^beta0 en B1=10^beta1 10^confint(logModel) #Betrouwbaarheidsintervallen voor B0 en B1 ``` Voor de interpretatie op originele schaal, vertrekken we terug van $$\widehat{GM(minsurv_i)} = \hat{B_0} * \hat{B_1}^{dosis_i}$$ Wat gebeurt er als we een nieuwe $dosis_j$ definiëren die één eenheid hoger ligt ($dosis_j=dosis_i+1$)? Dan krijgen we dit: \begin{aligned} \widehat{GM(minsurv_j)} &= \hat{B_0} * \hat{B_1}^{dosis_j}\\ &= \hat{B_0} * \hat{B_1}^{dosis_i+1}\\ &= \hat{B_0} * \hat{B_1}^{dosis_i}*\hat{B_1}\\ &= \widehat{GM(minsurv_i)}*\hat{B_1}\\ \end{aligned} Als de dosis een éénheid hoger ligt, zal het geometrisch gemiddelde dus vermenigvuldigd worden met een factor $\hat{B_1}= 10^{\hat{beta}}=10^{-0.27}=0.53$. Merk op dat vermeningvuldigen met een factor $0.53$ hetzelfde is als delen door een factor ${1\over0.53}=1.9$. ```{r} 1/(10^(logModel$coefficients)["dosis"]) # 1/B1 1/(10^(confint(logModel))["dosis",]) # Betrouwbaarheidsinterval voor 1/B1 ``` **De interpretatie van het effect van dosis op de overlevingstijd wordt dus:** *Het geometrisch gemiddelde van de overlevingstijd van dojovissen die worden blootgesteld aan een dosis van het gif die één milligram hoger is, is een factor `r 10^logModel$coef[2] %>% round(.,2)` van het geometrisch gemiddelde van de overlevingstijd van vissen die aan de lagere dosis worden blootgesteld (95% BI [`r 10^(confint(logModel)[2,]) %>%round(.,2) %>% sort`]). Dit effect is extreem significant (p < 0.001).* of anders gezegd *Het geometrisch gemiddelde van de overlevingstijd van dojovissen die aan een hogere dosis van het gif worden blootgesteld neemt gemiddeld met een factor `r 10^abs(logModel$coef[2]) %>% round(.,1)` af per mg extra gif die wordt toegediend (95% BI [`r 10^(abs(confint(logModel)[2,])) %>%round(.,1) %>% sort`]). Dit verschil is extreem significant (p < 0.001).* **Interpretatie van het 95%-betrouwbaarheidsinterval voor het effect van dosis op de overlevingstijd:** *Het interval [`r 10^(confint(logModel)[2,]) %>%round(.,2) %>% sort`] bevat met 95% waarschijnlijkheid de werkelijke factor waarmee het geometrisch gemiddelde van de overlevingstijd van dojovissen vermenigvuldigd wordt per extra mg gif die extra wordt toegediend.* of anders gezegd *Het interval [`r 10^(abs(confint(logModel)[2,])) %>%round(.,1) %>% sort`] bevat met 95% waarschijnlijkheid de werkelijke factor waarmee het geometrisch gemiddelde van de overlevingstijd van dojovissen afneemt per extra mg gif die extra wordt toegediend.* **Merk op dat de waarde 1 niet in het 95% BI ligt**. De nulhypothese van onze t-test veronderstelt dat de parameter op de log-schaal gelijk is aan nul. Aangezien p < 0,05, wordt de nulhypothese verworpen op het 5%-significantieniveau. Een 95%-betrouwbaarheidsinterval voor deze parameter op de log-schaal zal dus zeker niet de waarde 0 bevatten. Als we de parameter terugtransformeren, zal de waarde $10^0$ = 1 dus zeker niet in het gegeven interval liggen. Merk ook op dat de betrouwbaarheidsintervallen op een geometrisch gemiddelde (dus na terugtransformatie) niet meer symmetrisch rond de geschatte parameterwaarde liggen! Voor extra uitleg rond het verschil in interpretatie bij log-transformaties verwijzen we naar kennisclip "6.6 Afwijkingen van Modelveronderstellingen (Vervolg 1)" ## Schat het geometrisch gemiddelde van de overlevingstijd bij een dosis van 2 mg, met bijhorend 95%-betrouwbaarheidsinterval. ```{r} hulpdata <- data.frame(dosis = c(2)) # predictie op log-schaal: predict mean(log(y)) predict(logModel,hulpdata, interval="confidence") # interval="c" duidt aan dat we een confidence interval willen verkrijgen. ``` **Interpretatie van de voorspelling op logaritmische schaal** *Het gemiddelde 10-delig logaritme van de overlevingstijd in minuten van dojovissen die worden blootgesteld aan een dosis van 2 mg gif bedraagt `r predict(logModel,hulpdata, interval="confidence")[1] %>% round(.,2)` (95% BI [`r predict(logModel,hulpdata, interval="confidence")[-1] %>% round(.,2)`]). * ```{r} # geometrisch gemiddelde: 10^(mean(log(y))) 10^(predict(logModel,hulpdata, interval="confidence")) ``` **Interpretatie van de voorspelling op originele schaal** *Het geometrische gemiddelde van de overlevingstijd van dojovissen die worden blootgesteld aan een dosis van 2 mg gif bedraagt `r 10^predict(logModel,hulpdata, interval="confidence")[1] %>% round(.,1)` minuten (95% BI [`r 10^predict(logModel,hulpdata, interval="confidence")[-1] %>% round(.,1)`]).* ## Determinatiecoefficient Wat kunnen we besluiten uit de waarde voor de determinatiecoefficient? Deze waarde staat in de output van het lineaire regressiemodel (als u de "summary" opvraagt) bij "multiple R-Squared". ```{r} summary(logModel) ``` Met de predictorvariabele "dosis" kunnen we 31.28% van de totale variatie in de responsvariabele "10-delig logaritme van de overlevingstijd in minuten" verklaren.Het model heeft dus geen sterke voorspellende waarde, maar dit zegt echter niets over de kwaliteit van het model! Een grote determinatiecoefficient is een indicatie dat het lineaire model mogelijks tot goede predicties kan leiden. Een model met een kleine determinatiecoefficient blijft wel nuttig om associatie te bestuderen, zolang het de associatie correct modelleert (model assumpties). Intuïtief: de regressielijn capteert 31,28% van de totale variatie in het 10-delig logaritme van de overlevingstijd in minuten.