PC-les Algemene lineaire regressie: laatste oefening: regressie met drie onafhankelijke variabelen

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.

Data-exploratie

Lees de dataset poison.dat in via read.table. Verander de directory naar de folder waarin je de poison.dat file hebt opgeslagen. We zagen in de voorgaande practica al dat de overlevingstijd beter op logschaal wordt gemodelleerd

poison <- read.table("poison.dat", sep="", header = TRUE)
##we vormen de vissoort om in een factor. 
poison$soort <- as.factor(poison$soort)
#log transformatie
poison$logminsurv <- log(poison$minsurv)

plot(poison)

  • De overlevingstijd lijkt geassocieerd met het gewicht, soort en de dosis.

  • We observeren een sterke positieve associatie tussen de log-overlevingstijd en het gewicht.

  • Bij lage gewichten lijkt de log-overlevingstijd wat af te vlakken.

  • Daarnaast zien we ook dat het gewicht niet gelijk verdeeld is binnen elke dosis.

  • Er is ook een associatie tussen de gewicht en soort.

We bestuderen de effecten verder:

with(poison,scatter.smooth(logminsurv~gewicht,pch=as.double(soort),col=as.double(soort)))
legend("topleft",col=1:3,pch=1:3,legend=c("dojo","goud","zebra"))

boxplot(gewicht~dosis,poison,outline=FALSE,xlab="dosis",ylab="gewicht")
#outline = FALSE, geen outliers omdat we punten via stripchart toevoegen
stripchart(gewicht~dosis,poison,vertical=TRUE,pch=1,method = "jitter", ,add=TRUE)

plot(logminsurv~dosis,data=poison, col=as.double(soort),pch=as.double(soort))

plot(logminsurv~soort,poison,outline=FALSE)
stripchart(logminsurv~soort,poison,vertical=TRUE,col=as.double(poison$soort),pch=as.double(poison$soort),method = "jitter", ,add=TRUE)

  • We zien een heel sterke relatie tussen het gewicht en de log-overleving. Het gewicht is eveneens ongelijk tussen de soorten.

  • We observeren dat gewichten niet uniform verdeeld zijn over de dosisgroepen.

  • Het patroon wordt eveneens weerspiegeld in de plot waarbij de logsurvival 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 verschillende dosisrespons kan tonen.

We voeren nu een analyse uit die vergelijkbaar is met de drie aparte modellen die we hebben geschat tijdens het vorige practicum.

Schat alle effecten a.d.h.v. een model met interacties.

\[y_i=\beta_0+\beta_d x_{id} + \beta_g x_{ig} +\beta_{sa} x_{isa} +\beta_{sz} x_{isz} + \beta_{d-sa} x_{id}x_{isa} + \beta_{d-sb} x_{id}x_{isb} + \beta_{g-sa} x_{ig}x_{isa} + \beta_{g-sb} x_{ig}x_{isb} + \epsilon_i,\] met \(x_{isa}\) en \(x_{isz}\) dummy-variabelen die aangeven of een vis een goudvis is of een zebravis, respectievelijk. De referentieklasse is dus voor de soort dojovissen (als \(x_{isa}=0\) en \(x_{isz}=0\)). Verder is \(\epsilon_i \text{ i.i.d. } N(0,\sigma^2)\).

In de oefeningenles toonden we op bord aan dat dit model eigenlijk drie regressievlakken oplevert: 1 voor dojovissen, 1 voor goudvissen en 1 voor zebravissen. \[\begin{array}{ll} \text{dojovis } (x_{isa}=0\text{ en }x_{isz}=0):&E[y\vert \text{dojovis}]=\beta_0+\beta_d x_{id} + \beta_g x_{ig}\\ \text{goudvis } (x_{isa}=1\text{ en }x_{isz}=0):&E[y\vert \text{goudvis}]=\beta_0+\beta_{sa}+(\beta_d+ \beta_{d-sa}) x_{id} + (\beta_g+\beta_{g-sa}) x_{ig}\\ \text{zebravis } (x_{isa}=0\text{ en }x_{isz}=1):&E[y\vert \text{zebravis}]=\beta_0+\beta_{sz}+(\beta_d+ \beta_{d-sz}) x_{id} + (\beta_g+\beta_{g-sa}) x_{ig}\\ \end{array}\]

Het hoofdeffect voor soort zorgt dus dat elke soort een verschillend intercept heeft (\(\beta_0\), \(\beta_0 + \beta_{sa}\), \(\beta_0+\beta_{sz}\) voor dojo-, goud- en zebravissen, respectievelijk). De interactie tussen soort en dosis zorgt ervoor dat het dosiseffect (helling) verschillend kan zijn voor elke soort (\(\beta_d\), \(\beta_d+\beta_{d-sa}\), \(\beta_d+\beta_{d-sz}\) voor dojo-, goud- en zebravissen, respectievelijk). De interactie tussen soort en dosis zorgt ervoor dat het gewichtseffect (helling) verschillend kan zijn voor elke soort (\(\beta_g\), \(\beta_g+\beta_{g-sa}\), \(\beta_g+\beta_{g-sz}\) voor dojo-, goud- en zebravissen, respectievelijk). Het model heeft verder als voordeel dat het alle data gebruikt om de variantie te schatten van de residuen.

lmSam <- lm(logminsurv~ soort + dosis + gewicht + soort:dosis + soort:gewicht ,data=poison)
plot(lmSam)

We zien in de residuenplot opnieuw geen problemen mbt lineariteit en gelijkheid van variantie. De QQ-plot laat wel lichte afwijkingen zien van normaliteit. Maar deze zijn symmetrisch. Bovendien zijn er veel observaties (96) waardoor we kunnen verwachten dat de gemiddelden approximatief normaal verdeeld zullen zijn. We verwachten dus geen problemen mbt de inferentie.

Interpretatie van de parameters?

library(car)
## Warning: package 'car' was built under R version 3.4.4
## Loading required package: carData
Anova(lmSam, type = "III")
summary(lmSam)
## 
## Call:
## lm(formula = logminsurv ~ soort + dosis + gewicht + soort:dosis + 
##     soort:gewicht, data = poison)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62547 -0.21517 -0.04538  0.21384  0.63564 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      0.5749     0.4471   1.286 0.201855    
## soort1          -0.4090     0.6634  -0.616 0.539219    
## soort2           0.8956     0.7520   1.191 0.236924    
## dosis           -0.6648     0.1307  -5.086 2.08e-06 ***
## gewicht          0.7474     0.1933   3.866 0.000213 ***
## soort1:dosis     0.1171     0.1944   0.602 0.548640    
## soort2:dosis     0.1761     0.2196   0.802 0.424863    
## soort1:gewicht   0.3311     0.2660   1.245 0.216506    
## soort2:gewicht  -0.6964     0.3861  -1.804 0.074754 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.285 on 87 degrees of freedom
## Multiple R-squared:  0.7429, Adjusted R-squared:  0.7192 
## F-statistic: 31.42 on 8 and 87 DF,  p-value: < 2.2e-16

De model parameters kunnen dus als volgt worden geinterpreteerd:

  • Intercept: log getransformeerde geometrische gemiddelde van de overlevingstijd bij dojovissen met een gewicht van 0 gram die zich in een aquarium zonder gif bevinden. We zien dat het opnieuw een interpretatie betreft die biologisch geen steek houdt gezien het op extrapolatie berust!
  • Hoofdeffect van dosis: De log-overlevingstijd bij dojovissen is gemiddeld 0.665 lager als de dosis van het gif 1mg/l hoger is, na correctie voor gewicht. (Het hoofdeffect bij een interactie tussen een continue predictor en een factor heeft dus betrekking op het lineaire effect van de continue predictor in de referentie groep)
  • Hoofdeffect van gewicht: De log-overlevingstijd bij dojovissen is gemiddeld 0.747 hoger bij vissen die 1 g meer weegt, na correctie voor dosis.
  • soort1: Het gemiddeld verschil in log overlevingstijd tussen goudvissen en dojovissen die beiden een gewicht hebben van 0 gram en zich in een aquarium zonder gif bevinden. Het geeft m.a.w. aan hoe het intercept wijzigt tussen goud en dojovissen.
  • Dosis x soort1 interactie: De log-overlevingstijd bij goudvissen neemt gemiddeld met 0.117 minder snel af per mg/l toename van de dosis dan de afname bij dojovissen, na correctie voor gewicht. Het effect van de dosis is dus minder groot bij goudvissen. We dat dit verschil niet significant is. (Een interactie tussen een continue predictor en een factor heeft dus betrekking op hoe de helling van de continue predictor wijzigt tussen een groep en de referentie groep).
  • Dosis x soort2 interactie: analoge interpretatie als voor soort1
  • De interacties voor gewicht en soort hebben een analoge interpretatie.

! Stuk hieronder bevat modelselectie en dit is nog niet gezien, dus dit heb ik verwijderd uit de oefening.