--- title: "Additief lineair model: uitgewerkt voorbeeld Dojovissen" 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 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 (`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. ## Laad de libraries ```{r} library(dplyr) library(ggplot2) #install.packages("GGally") library(GGally) library(car) library(multcomp) ``` # Data-exploratie Lees de dataset poison.dat in via `read.table`. We zagen in de voorgaande practica al dat de overlevingstijd beter op log2-schaal wordt gemodelleerd. ```{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)) poison %>% ggpairs ``` - De overlevingstijd lijkt geassocieerd met het gewicht, soort en de dosis. - We observeren een sterke positieve associatie tussen de log2-overlevingstijd en het gewicht. - Bij lage gewichten lijkt de log2-overlevingstijd wat af te vlakken. - Daarnaast zien we ook dat het gewicht niet gelijk verdeeld is binnen elke dosis. - Er lijkt ook een associatie tussen de gewicht en soort. We bestuderen de effecten verder: ```{r} poison %>% ggplot(aes(x = dosis, y = log2minsurv, col = soort)) + geom_point() + stat_smooth(method = "loess",col="red") + # fit een kromme door de puntenwolk (rode lijn) ylab("Overlevingstijd (log2)") + xlab("Dosis (mg)") + theme_bw() poison %>% ggplot(aes(x = gewicht, y = log2minsurv, col = soort)) + geom_point() + stat_smooth(method = "loess",col="red") + # fit een kromme door de puntenwolk (rode lijn) ylab("Overlevingstijd (log2)") + xlab("Gewicht") + theme_bw() scatterplot1 <- poison %>% ggplot(aes(x = dosis, y = gewicht)) + geom_point() + theme_bw() scatterplot1 boxplot2 <- poison %>% ggplot( aes(x = soort, y = log2minsurv, col = soort) ) + geom_boxplot(outlier.shape = NA) + geom_jitter() + theme_bw() boxplot2 boxplot3 <- poison %>% ggplot( aes(x = soort, y = gewicht, col = soort) ) + geom_boxplot(outlier.shape = NA) + geom_jitter() + theme_bw() boxplot3 ``` - We zien een heel sterke relatie tussen het gewicht en de log2-overlevingstijd. Het gewicht is eveneens ongelijk tussen de soorten. - We observeren dat gewichten niet uniform verdeeld zijn over de verschillende dosissen van het gif in de study. - Het patroon wordt eveneens weerspiegeld in de plot waarbij de log2survival wordt uitgezet t.o.v. dosis. - Deze plot geeft ook weer dat er een effect is van de soort op de overleving. Gezien de onderzoekers op basis van voorgaande studies vermoeden dat het effect van de dosis kan variëren van soort tot soort zouden we in de modellen interacties moeten voorzien tussen dosis en soort, zodat elke soort een verschillende dosisrespons kan tonen. Wegens didactische redenen komt het modelleren van interacties pas later aan bod. Om de data toch correct te modelleren, zullen we het dosiseffect in dit practicum apart modelleren per soort. We starten hierbij met de dojovissen. # Schat het effect van de dosis na correctie voor het gewicht bij dojovissen. ## Aanpak van vorig practicum met enkelvoudige regressie Model met enkel dosis effect: $$ y_i=\beta_0+\beta_d x_d + \epsilon_i, $$ ```{r} lm0 <- lm(log2minsurv~dosis, data = poison %>% filter(soort == 0)) summary(lm0) ``` ## Aanpak waarbij we corrigeren voor het gewicht We voegen nu een extra predictor (onafhankelijke variabele of covariaat) toe aan het model, namelijk 'gewicht'. $$ y_i=\beta_0+\beta_d x_d + \beta_g x_g + \epsilon_i, $$ met $\epsilon_i \text{ i.i.d. } N(0,\sigma^2)$. We schatten de parameters van het model, evalueren de aannames en doen besluitvorming in R: ```{r} lmDojo <- lm(log2minsurv~ dosis + gewicht, data = poison %>% filter(soort == 0)) summary(lmDojo) plot(lmDojo) ``` ### Aannames - We zien geen problemen m.b.t. lineariteit. De residuen zijn mooi verspreid rond 0 voor alle gefitte waarden. - We zien geen problemen mbt de aanname van gelijkheid van variantie. De residuen hebben een gelijke spreiding voor alle gefitte waarden. - De data lijken wat afwijkingen van normaliteit te vertonen in de linkerstaart. - We gaan eerst na welke afwijkingen we kunnen verwachten in een dataset met evenveel residuen (`r sum(poison$soort == 0)`) als in onze dataset wanneer de residuen worden getrokken uit een normale verdeling met dezelfde standaard deviatie als wat we observeerden in onze dataset. ```{r} set.seed(1406) nobs <- sum(poison$soort == 0) data.frame( y = c(lmDojo$res, rnorm(nobs*8, sd = sigma(lmDojo) ) ), label = rep( c("orig", paste0("sim",1:8)), each = 39)) %>% ggplot(aes(sample = y)) + geom_qq() + geom_qq_line() + facet_wrap(~ label) ``` - We observeren dat de afwijkingen die we zien in de QQ-plot voor de residuen van het lineaire model voor dojovissen in de lijn ligt van de afwijkingen die we kunnen verwachten wanneer de residuen normaal verdeeld zouden zijn. De afwijkingen zijn dus niet extreem en bovendien beschikken we reeds over `r nobs` observaties zodat we er zeker van uit mogen gaan dat de parameterschatters benaderend normaal verdeeld zijn. ### Inferentie We observeren dat het effect van de dosis veel significanter is na correctie voor gewicht. 1. Dat is logisch want uit de data-exploratie bleek dat gewicht geassocieerd was met de overleving. Wanneer we een extra variabele opnemen in het model die veel van de variabiliteit in de uitkomstvariabele kan verklaren, zal de variabiliteit van de residuen afnemen, waardoor de standaard errors op de parameterschatters sterk af kunnen nemen. 2. Bovendien observeerden we in de figuren van gewicht en dosis, dat de gewichten niet goed gebalanceerd zijn tussen de verschillende dosissen, waardoor het effect van gewicht de schatting van het dosiseffect wat verstoort. Aangezien we ook een associatie tussen gewicht en overlevingstijd en dosis en overlevingstijd zien, zullen we het effect van de dosis dus beter kunnen schatten als we corrigeren voor het gewicht. Uit (1) en (2) volgt dus dat als we corrigeren voor gewicht dat we het dosiseffect dus beter en met een grotere precisie kunnen schatten! Het dosiseffect krijgt in dit model de betekenis van de gemiddelde verandering in log2-overlevingstijd tussen twee groepen dojovissen met hetzelfde gewicht die worden blootgesteld aan een dosis die 1 mg/l verschilt: \begin{eqnarray} \hat{\mu_1}&=& \beta_0 + \beta_d x_{1d} + \beta_g x_g \text{ (gemiddelde log2-overlevingstijd voor dosis 1 bij een bepaald gewicht)}\\ \hat{\mu_2}&=& \beta_0 + \beta_d x_{2d} + \beta_g x_g \text{ (gemiddelde log2-overlevingstijd voor dosis 2 bij hetzelfde gewicht)}\\ \hat \mu_2- \hat \mu_1&=&\beta_0 + \beta_d x_{2d} + \beta_g x_g - (\beta_0 + \beta_d x_{1d} + \beta_g x_g) \text{ (verschil in gemiddelde log2-overlevingstijd tussen dosis 2 en dosis 1)}\\ \hat \mu_2-\hat \mu_1&=&\beta_d (x_{2d}-x_{1d}) \end{eqnarray} Als het concentratieverschil $(x_{2d}-x_{1d})=1$ mg/l dan bekomen we inderdaad dat $\hat \mu_2-\hat \mu_1=\beta_d$. $\beta_d$ kwantificeert dus inderdaad het verschil in gemiddelde log2-overlevingstijd bij een dosisverschil van 1 mg/l voor dojovissen met hetzelfde gewicht. Dus de log2-overlevingstijd bij dojovissen is gemiddeld `r format(abs(lmDojo$coef["dosis"]),digits=3)` lager bij vissen die blootgesteld worden aan een dosis van het gif die 1mg/l hoger is, na correctie voor gewicht. Als we data modelleren op log2-schaal bekomen we geometrische gemiddelden, $\widehat{minsurv}$ na terugtransformatie: \begin{eqnarray} \log2 \widehat{minsurv}&=&\beta_0+ \beta_d x_{d} + \beta_g x_g\\ \widehat{minsurv}&=&2^{(\beta_0+ \beta_d x_{d} + \beta_g x_g)}\\ \widehat{minsurv}&=&2^{\beta_0}2^{\beta_d x_{d}}2^{\beta_g x_g}\\ \end{eqnarray} Voor twee groepen dojovissen met hetzelfde gewicht die aan een verschillende dosis worden blootgesteld bekomen we dan: \begin{eqnarray} \frac{\widehat{minsurv}_2}{\widehat{minsurv}_1}&=& \frac{2^{\beta_dx_{2d}}}{2^{\beta_dx_{1d}}}\\ &=& 2^{(\beta_dx_{2d} - \beta_dx_{1d})} \\ &=&2^{\left[\beta_d (x_{2d}-x_{1d})\right]}\\ \end{eqnarray} ```{r} 2^(lmDojo$coef) 2^(confint(lmDojo)) ``` ## Conclusie De dosis van gif heeft een extreem significant effect op de overlevingstijd bij dojovissen ($p<<0.001$). Het geometrisch gemiddelde van de overlevingstijd van dojovissen met hetzelfde gewicht wordt gehalveerd (factor $2^{\beta_d}=$ `r format(2^(lmDojo$coef["dosis"]),digits=3)`) als de dosis van het gif met 1 mg/l toeneemt (95% BI [`r paste(format(2^(confint(lmDojo)["dosis",]),digits=3),collapse=",")`]). Het effect van het gewicht op de overlevingstijd is eveneens extreem significant ($p<<0.001$). Het geometrisch gemiddelde van de overlevingstijd van een dojovis die een 1 gram meer weegt dan een andere dojovis en die aan dezelfde dosis wordt blootgesteld, is dubbel zo groot (factor $2^{ \beta_g}$= `r format(2^(lmDojo$coef["gewicht"]),digits=3)`, 95% BI [`r paste(format(2^(confint(lmDojo)["gewicht",]),digits=3),collapse=",")`]).