We modelleren de data met een hoofdeffect voor gif en behandeling en een gif \(\times\) behandeling interactie.
\[\begin{array}{lcl} y_i &=& \beta_0 + \\ &&\beta_{II} x_{iII} + \beta_{III} x_{iIII} + \\ && \beta_{B} x_{iB} + \beta_{C} x_{iC} + \beta_{D} x_{iD} + \\ &&\beta_{II:B}x_{iII}x_{iB} + \beta_{II:C}x_{iII}x_{iC} + \beta_{II:D}x_{iII}x_{iD} + \\ &&\beta_{III:B}x_{iIII}x_{iB} + \beta_{III:C}x_{iIII}x_{iC} + \beta_{III:D}x_{iIII}x_{iD} + \epsilon_i \end{array}\]met \(i = 1, \ldots, n\), \(n=48\), en, \(x_{iII}\), \(x_{iIII}\), \(x_{iB}\), \(x_{iC}\) en \(x_{iD}\) dummy variabelen voor respectievelijk gif II, III, behandeling B, C, en D.
rats1 <- lm(time~poison*treat, rats)
summary(rats1)
##
## Call:
## lm(formula = time ~ poison * treat, data = rats)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2500 -0.4875 0.0500 0.4312 4.2500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.1250 0.7457 5.532 2.94e-06 ***
## poisonII -0.9250 1.0546 -0.877 0.3862
## poisonIII -2.0250 1.0546 -1.920 0.0628 .
## treatB 4.6750 1.0546 4.433 8.37e-05 ***
## treatC 1.5500 1.0546 1.470 0.1503
## treatD 1.9750 1.0546 1.873 0.0692 .
## poisonII:treatB 0.2750 1.4914 0.184 0.8547
## poisonIII:treatB -3.4250 1.4914 -2.297 0.0276 *
## poisonII:treatC -1.0000 1.4914 -0.671 0.5068
## poisonIII:treatC -1.3000 1.4914 -0.872 0.3892
## poisonII:treatD 1.5000 1.4914 1.006 0.3212
## poisonIII:treatD -0.8250 1.4914 -0.553 0.5836
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.491 on 36 degrees of freedom
## Multiple R-squared: 0.7335, Adjusted R-squared: 0.6521
## F-statistic: 9.01 on 11 and 36 DF, p-value: 1.986e-07
plot(rats1)
De errors zijn heteroscedastisch. De residu plot suggereert een relatie tussen gemiddelde en variantie. De QQ-plot suggereert dat de verdeling mogelijks bredere staarten heeft dan de normale verdeling.
rats %>%
ggplot(aes(x=treat,y=log2(time))) +
geom_boxplot(outlier.shape=NA) +
geom_jitter() +
facet_wrap(~poison)
rats2 <- lm(time %>% log2~poison*treat, rats)
plot(rats2)
Log transformatie verwijdert heteroscedasticiteit niet volledig.
rats %>%
ggplot(aes(x=treat,y=1/time)) +
geom_boxplot(outlier.shape=NA) +
geom_jitter() +
facet_wrap(~poison) +
ylab ("rate of dying (1/h)")
rats3 <- lm(1/time~poison*treat, rats)
plot(rats3)
De reciproke transformatie lijkt beter. Transformaties bemoeilijken soms de interpretatie. Hier kan de reciproke transformatie echter worden geïnterpreteerd als de “snelheid van sterven” (rate of dying).