--- title: "Enkelvoudige lineaire regressie" output: html_document --- # PC-les 4: Enkelvoudige lineaire regressie ## Enkelvoudige lineaire regressie met een continue predictor: log-transformatie Herinner u de opgave van de gezamelijke oefeningenles: Bij 96 vissen (dojovissen, goudvissen en zebravissen) werd de resistentie tegen het vergif 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). In deze aparte oefeningenles gaan we verder in op de onderzoeksvraag: "Wat is de associatie tussen dosis en overlevingstijd?" In de gezamelijke oefeningenles zagen we immers dat de residuen geen normale verdeling volgen en ook assumptie van homoscedasticiteit geschonden is. In deze les zullen we beroep doen op een (log-)transformatie van de afhankelijke variabele om toch aan deze assumpties te voldoen. Lees de dataset poison.dat opnieuw 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 ``` #### 1. Waarom een log-transformatie? Herinner u: bij het model van vorige les hadden de residuen een korte linkse staart en een lagen rechtse staart. ```{r} model1 = lm(minsurv~dosis) qqnorm(resid(model1)) qqline(resid(model1)) ``` 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. Dit zie je ook wanneer je een histogram maakt van de uitkomstvariabele voor en na logtransformatie. ```{r} hist(minsurv, main="Histogram van de overlevingstijd") ``` ```{r} #Log-transformatie van de afhankelijke variabele minsurv log.minsurv <- log(minsurv) hist(log.minsurv, main="Histogram van de log-overlevingstijd") ``` 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(fitted(model1),sqrt(abs(resid(model1)))) lines(lowess(x=fitted(model1),y=sqrt(abs(resid(model1)))),col="red") ``` Dit is een extra indicatie dat een log-transformatie zinvol kan zijn: een logtransformatie zal een additieve errorstructuur omzetten naar een multiplicatieve structuur. #### 2. Voer een log-transformatie uit van de afhankelijke variabele en voer een lineaire-regressieanalyse uit voor de getransformeerde variabele. ```{r} #fit een lineair regressiemodel met 'log-minsurv' als afhankelijke en 'dosis' als onafhankelijke variabele logModel <- lm(........) logModel ``` 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. ```{r} log(100, base=10) ``` ### 3. Ga zelf na of de voorwaarden voor lineaire regressie nu wel voldaan zijn. ```{r} par(mfrow=c(2,2)) plot(logModel) ``` Ga hier allevier de voorwaardes na! ... #### 4. Geef de nul- en alternatieve hypothese voor de t-test in de modeloutput die geassocieerd is met het dosis-effect enerzijds en de nul- en alternatieve hypothese voor de t-test in de modeloutput die geassocieerd is met het intercept anderzijds. **Nulhypothese voor het dosiseffect:** ... **Alternatieve hypothese voor het dosiseffect:** ... **Nulhypothese voor het intercept:** ... **Alternatieve hypothese voor het intercept:** ... #### 5. Interpretatie van het dosis-effect voor het model met een log-getransformeerde uitkomstvariabele. **1. Interpretatie op de schaal van getransformeerde waarden:** $log(minsurv_i) = \hat{\beta}_0 + \hat{\beta}_1 dosis_i + \epsilon_i$ $log(minsurv_i) = 2.105 - 0.5112*dosis_i + \epsilon_i$ met $\epsilon_i \sim N(0,0.5104^2)$. **Interpretatie:** als de dosis met 1 miligram stijgt, dan daalt het gemiddelde van het natuurlijk logaritme van de overlevingstijd met 0.5112. **2. Interpretatie op de originele schaal:** Vanwege de transformatie modelleren we het gemiddelde van het logaritme van de responsvariabele, en interpreteer je de parameters in deze context, zoals hierboven. Als men de exponent neemt van de geschatte parameters, gelden ze in termen van het **geometrisch gemiddelde** van de responsvariabele op de originele schaal. Het geometrisch gemiddelde van $Y$ is gelijk aan: $e^{rekenkundig\ gemiddelde(log(y_i))} = exp(\frac{\sum_{i=1}^n log(y_i)}{n}).$ **Opmerking:** "exp(...)" betekent hetzelfde als "e^..."! Zoals reeds aangehaald duidt de **`log`-transformatie** in de statistiek typisch op het **natuurlijk logaritme** (i.e. met **basis gelijk aan $e$** en niet met basis $10$). Op deze manier interpreteren we dan ook de parameters. $\widehat{minsurv}_i = exp(2.105 - 0.5112*dosis_i) = \frac{e^{2.105}}{e^{0.5112*{dosis}_i}}$ (we maakten gebruik van de rekenregel: $e^{a-b} = e^a/e^b$) **Interpretatie:** als de dosis van het vergif stijgt met één miligram, dan daalt het **geometrisch gemiddelde** van de overlevingstijd van dojovissen, goudvissen en zebravissen met een **factor** van $e^{-0.5112} = 1/e^{0.5112} = 0.60$. Doe het volgende zelfstandig om bovenstaande interpretatie te verduidelijken: - Bereken het verwachte geometrische gemiddelde van de overlevingstijd voor een dosis van 1 miligram. ```{r} ``` - Doe hetzelfde voor een dosis van 2 miligram. ```{r} ``` - Vergelijk deze waarden. Wat merk je op? .... #### 6. Probeer nu zelf eens een correcte en volledige interpretatie van het effect van de dosis op de verwachte uitkomst te geven, samen met zijn bijhorende 95%-betrouwbaarheidsinterval, zowel op de log-schaal als op de originele schaal. **Interpretatie op de log-schaal** ```{r} summary(logModel) confint(logModel) ``` Interpretatie van de parameter `dosis` op de log-schaal: ... Interpretatie van het 95%-betrouwbaarheidsinterval voor `dosis` op de log-schaal: ... **Interpretatie op de originele schaal** ```{r} exp(summary(logModel)$coefficients[,"Estimate"]) exp(confint(logModel)) ``` Interpretatie van de parameter `dosis` op de originele schaal: ... Interpretatie van het 95%-betrouwbaarheidsinterval voor `dosis` op de originele schaal: ... **Gegeven een model waarvan de uitkomstvariabele log-getransformeerd is, als ik u voor een bepaalde parameter een 95%-betrouwbaarheidsinterval op de originele schaal zou geven en u zou zeggen dat de p-waarde voor die parameter kleiner is dan 0,05, welke waarde ligt er dan met zekerheid NIET in het gegeven 95%-betrouwbaarheidsinterval?** ... Merk ook op dat de betrouwbaarheidsintervallen op een geometrisch gemiddelde (dus na terugtransformatie) niet meer symmetrisch rond de geschatte parameterwaarde liggen! #### 7. Geef een correcte en volledige interpretatie van het intercept met zijn bijhorende 95%-betrouwbaarheidsinterval, zowel op de log-schaal als op de originele schaal ```{r} summary(logModel) confint(logModel) ``` **Interpretatie op de log-schaal** Interpretatie van het intercept op de log-schaal: ... Interpretatie van het 95%-betrouwbaarheidsinterval voor het intercept op de log-schaal: ... ```{r} exp(confint(logModel)) exp(summary(logModel)$coefficients[,"Estimate"]) ``` **Interpretatie op de originele schaal** Interpretatie van het intercept op de originele schaal: ... Interpretatie van het 95%-betrouwbaarheidsinterval voor het intercept op de originele schaal: ... **Is deze interpretatie van het intercept biologisch relevant?** ... #### 8. Schat het geometrisch gemiddelde van de overlevingstijd bij een dosis van 2 mg. Geef een bijhorend 95%-betrouwbaarheidsinterval. ```{r} hulpdata<-data.frame(dosis = c(2)) #predictie op log-schaal: predict mean(log(y)) predict(logModel,hulpdata, interval="confidence") ``` **Hoe interpreteren we dit 95% betrouwbaarheidsinterval?** ... ```{r} #geometrisch gemiddelde: exp(mean(log(y))) exp(predict(logModel,hulpdata, interval="confidence")) ``` Merk op dat $e^{0.9122891}=2.490016$ en $e^{1.252174}=3.497941$. **Hoe interpreteren we het 95% betrouwbaarheidsinterval na terugtransformatie?** ... #### 9. Wat kunnen we besluiten uit de waarde voor de meervoudige correlatiecoëfficiënt? Deze waarde staat in de output van het lineaire regressiemodel (als u de "summary" opvraagt) bij "multiple R-Squared". ```{r} summary(logModel) ``` ... Intuïtief: de regressielijn capteert 10,88% van de totale variatie in het natuurlijk logaritme van de overlevingstijd in minuten. #### 10. Uitwerken van effectgroottes, standaardfouten, T-waarden en p-waarden In dit laatste deel zullen we zelf de effectgroottes, standaardfouten, T-waarden en p-waarden van de parameters van het lineair regressiemodel berekenen. **Voor het effect van dosis** Op basis van hoofdstuk 5.3 van de cursus weten we dat: $$\hat{\beta}_1 = \frac{\sum_{i=1}^{n}(y_i-\bar{y})(x_i-\bar{x})}{\sum_{i=1}^{n}(x_i-\bar{x})^2}$$ Voor ons voorbeeld: $$\hat{\beta}_1 = \frac{\sum_{i=1}^{n}(log.minsurv_i-\overline{log.minsurv_i})(dosis_i-\overline{dosis})}{\sum_{i=1}^{n}(dosis_i-\overline{dosis})^2}$$ De grootte van het dosis-effect is dus gelijk aan: ```{r} beta1 = sum((log.minsurv-mean(log.minsurv))*(dosis-mean(dosis)))/sum((dosis-mean(dosis))^2) beta1 summary(logModel)$coefficients["dosis","Estimate"] ``` De geschatte variantie op het dosis-effect kunnen we als volgt berekenen: $$\hat{\sigma}^2_{\hat{\beta_1}} = \frac{\hat{\sigma}^2}{\sum_{i=1}^{n}(x_i-\bar{x})^2}$$ De geschatte standaardfout op $\hat{\beta_1}$ is dus: $$\hat{\sigma}_{\hat{\beta_1}} = \frac{\hat{\sigma}}{\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2}}$$ ```{r} sigma_beta1 = sigma(logModel)/sqrt(sum((dosis-mean(dosis))^2)) summary(logModel)$coefficients["dosis","Std. Error"] ``` Onze geschatte t-waarde is simpelweg: $$t = \frac{\hat{\beta}_1}{\hat{\sigma}_{\hat{\beta_1}}}$$ ```{r} tval_beta1 = beta1/sigma_beta1 tval_beta1 summary(logModel)$coefficients["dosis", "t value"] ``` De p-waarde is de kans om een even extreem of nog extremer resultaat te observeren als ons geveven resultaat `r tval_beta1`. We weten dat als aan al onze assumpties voldaan is, de teststatistiek onder de nulhypothese ($\beta_1 = 0$) een t-verdeling volgt met "# observaties - 2" = 96 - 2 = 94 vrijheidsgraden. De kans dat een t-verdeelde variabele met 94 vrijheidsgraden kleiner is dan `r tval_beta1` is gelijk aan: ```{r} pt(tval_beta1, df=94, lower.tail=TRUE) ``` Aangezien we tweezijdig testen, moeten we hier ook de kans bij optellen dat onze t-verdeelde variabele groter is dan `r -tval_beta1`. Deze kans is gelijk aan: ```{r} pt(-tval_beta1, df=94, lower.tail=FALSE) ``` De p-waarde is dus gelijk aan de kans dat een t-verdeelde toevalsveranderlijke extremer (d.w.z. groter dan 3.38 of kleiner dan -3.38) is: ```{r} pt(tval_beta1, df=94, lower.tail=TRUE)+pt(-tval_beta1, df=94, lower.tail=FALSE) ``` Aangezien een t-verdeling symmetrisch is, kan de p-waarde ook gewoon als volgt berekend worden: ```{r} pval_beta1 = 2*pt(-abs(tval_beta1), df=94, lower.tail=TRUE) #Met -abs(tval_beta1) nemen we de absolute waarde en maken de t-waarde dan negatief, zodat je onafhankelijk van het teken van de t-waarde altijd in de linkerstaart test. pval_beta1 summary(logModel)$coefficients["dosis", "Pr(>|t|)"] ``` **Probeer dit nu zelf eens voor het intercept** Op basis van hoofdstuk 5.3 van de cursus weten we dat: $\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}$ Voor ons voorbeeld: $\hat{\beta}_0 = \overline{log.minsurv} - \hat{\beta}_1 \overline{dosis}$ ```{r} beta0 = ... beta0 summary(logModel)$coefficients["(Intercept)","Estimate"] ``` De geschatte standaardfout op het intercept kunnen we als volgt berekenen: $$\hat{\sigma}_{\hat{\beta_0}} = \hat{\sigma}*\sqrt{\frac{\sum_{i=1}^{n}(x_i)^2}{n*\sum_{i=1}^{n}(x_i-\bar{x})^2}}$$ ```{r} sigma_beta0 = sigma(logModel)*sqrt(sum(dosis^2)/(96*sum((dosis-mean(dosis))^2))) sigma_beta0 summary(logModel)$coefficients["(Intercept)","Std. Error"] ``` Onze geschatte t-waarde is opnieuw: $$t = \frac{\hat{\beta}_0}{\hat{\sigma}_{\hat{\beta_0}}}$$ ```{r} tval_beta0 = ... tval_beta0 summary(logModel)$coefficients["(Intercept)", "t value"] ``` Onze p-waarde berekenen we opnieuw als de kans dat een t-verdeelde variabele met 94 vrijheidsgraden nog extremer is dan `r tval_beta0`: ```{r} pval_beta0=... pval_beta0 summary(logModel)$coefficients["(Intercept)", "Pr(>|t|)"] ```