--- title: "Algemeen lineair model met interactietermen: half uitgewerkte oefening" output: html_document: code_download: yes highlight: tango number_sections: yes theme: cosmo toc: yes toc_float: yes pdf_document: toc: yes word_document: toc: yes --- # Achtergrond, libraries laden en data-exploratie We herhalen aan het begin van deze oefening nogmaals de situatie en bijhorende data-exploratie. Net zoals in de vorige oefening zullen jullie hier een additief model opstellen, waarbij het effect van de dosis na correctie voor het gewicht bij **goudvissen** (soort 1) geschat zal worden. Resistentie tegen het gif EI-43,064 wordt getest bij 96 vissen (dojovissen (0), goudvissen (1) en zebravissen (2)). Elke vis wordt apart in een aquarium gestopt die een bepaalde `dosis` (in mg) van het gif bevat. Naast de overlevingstijd in minuten (de uitkomst, `minsurv`) werd ook het `gewicht` van de vis gemeten (in gram). De onderzoekers weten uit vorige experimenten dat de overlevingstijd vaak sterk afhangt van het gewicht en dat de resistentie dikwijls soortafhankelijk is. De onderzoekers wensen inzicht te krijgen in het effect van de dosis op de overlevingstijd en of resistentie tegen het gif verschillend is bij de verschillende soorten. ***Vorige oefening hebben we enkel de hoofdeffecten in het model opgenomen, zonder interactietermen. Dit was louter voor didactische redenen. Normaal gezien moet je altijd eerst de interactie meenemen in het model, en bekijken of deze significant is. Indien je deze interactie niet weerhoudt, dan pas kan je een model opstellen met enkel de hoofdeffecten. We verwachten dus, wanneer je een analyse doet, dat je de methode van deze file gebruikt, in plaats van het lineair additief model.*** ```{r} #gebruik de install.packages() functie voor packages die nog niet geïnstalleerd zijn. library(dplyr) library(ggplot2) library(GGally) library(car) library(multcomp) ``` ```{r} poison <- read.table("https://raw.githubusercontent.com/statOmics/statistiekBasisCursusData/master/practicum8/poison.dat", sep="", header = TRUE) # We vormen de vissoort om in een factor en log2 transformeren de overlevingstijd poison <- poison %>% mutate(soort = as.factor(soort), log2minsurv = log2(minsurv)) ``` In de vorige oefening onderzochten we de invloed van dosis en gewicht op de overlevingstijd voor dojovissen We zullen nu een gelijkaardige analyse uitvoeren, maar we zullen deze toepassen op de goudvissen (soort==1), maar zullen nu interactietermen toevoegen. # Vraag 1: Model met interactietermen Stel het model op met interactietermen. ```{r} lmGoudInt <- lm(..., data = poison %>% filter(soort == 1)) ``` # Vraag 2: Ga de voorwaarden van het model na (Meerkeuzevraag) ```{r} plot(lmGoudInt) ``` In het vorige model was aan alle voorwaarden voldaan. Welke van de voorwaarden is niet meer voldaan voor dit model? 1. Onafhankelijke gegevens 2. Lineariteit 3. Normaal verdeelde residuen 4. Geen van bovenstaande # Vraag 3: Interpretatie (Meerkeuzevraag) Welke interpretatie van de schattingen is juist ```{r} summary(lmGoudInt) ``` 1. De verwachte log2-overlevingstijd van een vis van 2 g zal gemiddeld 1.66 lager liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt. 2. De verwachte log2-overlevingstijd van een vis van 2 g zal gemiddeld 0.37 hoger liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt. 3. De verwachte log2-overlevingstijd van een vis van 2 g zal gemiddeld 0.51 lager liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt. 4. De verwachte log2-overlevingstijd van een vis van 2 g zal gemiddeld 0.66 lager liggen, vergeleken met een vis met hetzelfde gewicht, indien de dosis 1 eenheid hoger ligt. # Vraag 4: Test of er een effect is van dosis. Test of er een effect is van dosis (hoofdeffect en interactie-effect samen). Stel daarvoor eerst een lineair model op zonder dosis en vergelijk dit model met het lmGoudInt model via de anova-functie ```{r} lmGoudH0 <- lm(..., data = ... %>%filter(...)) tabel_anova<-anova(...) tabel_anova ``` # Vraag 5: Interactie testen Test of er een interactie-effect aanwezig is. Geef ook de tabel Type III kwadratensommen. Tip: gebruik de `Anova` functie van de `car` package. ```{r} tabel_kwadratensommen<-... tabel_kwadratensommen ``` # Vraag 6: Marginale dosis-effect (Leesopdracht) De interactieterm is niet significant. Om verder te gaan met de analyse hebben we nu twee opties. Volgens de conventionele methode zouden we de interactie term uit het model kunnen laten. In dit geval zouden we hetzelfde additief model opstellen als in de vorige oefening. Merk op het additief model dat we opstelden in de vorige oefening dus wel een correct model was om deze data te analyseren, maar dat we altijd eerst moeten nagaan of er al dan niet een interactie-effect aanwezig is. Een andere optie is om de interactie term toch in het model te laten. In dat geval kunnen we een analyse uitvoeren met een **marginaal model**. Het marginaal model schat dan het marginale dosis effect; Dit is het gemiddelde effect over alle vissen heen en is gelijk aan het effect voor een vis met gemiddeld gewicht, m.a.w. $\beta_d +\beta_{d:g}*\bar{x}_g$. Beide opties zijn correct en zullen voor deze analyse bijna exact hetzelfde resultaat bekomen. In sommige (andere) toepassingen verkiest het marginale model de voorkeur. Bijvoorbeeld, bij "high-throughput" experimenten (experimenten met honderden of duizenden predictorvariabelen, zoals transcriptomics en proteomics experimenten) weten we dat we zelden genoeg statistische power hebben om interactie-effecten te detecteren, ook al zijn deze in aanwezig in de realiteit. In zulke gevallen zal het vaak correcter zijn om de interactie-term toch in het model te laten en verder te gaan met een marginaal model, dan de interactietermen gewoon te verwijderen. Om het marginaal model voor deze oefening op te stellen, moeten we het dosis effect berekenen voor een vis met een gemiddeld gewicht: ```{r} wBar <- ... #Bereken het gemiddelde gewicht van de goudvissen summary(lmGoudInt) 0.3691+wBar*(-0.5146) ``` # Vraag 7: Schat nu het marginaal effect met behulp van de glht-functie en geef de summary met de p-waarde. ```{r} marginaleDosisEffect <- glht(..., linfct = c("...")) summary1<-... summary1 ``` # Vraag 8: Geef de geschatte coëfficienten met betrouwbaarheidsinterval op originele schaal (niet log-getransformeerd) en interpreteer. (Meerkeuzevraag) ```{r} 2^confint(marginaleDosisEffect)$confint ``` 1. De gemiddelde overlevingstijd zal 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$). 2. Het geometrisch gemiddelde van de overlevingstijd zal een factor 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$). 3. Het marginaal gemiddelde van de overlevingstijd zal 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$). 4. Het marginaal geometrisch gemiddelde van de overlevingstijd zal een factor 0.578 lager liggen bij vissen die onderworpen worden aan een dosis die 1 mg/l hoger ligt (95% CI [0.4385462 0.7629894`]), na correctie voor gewicht ($p<<0.001$).