1 Gewicht von Hühnerküken

library(Rcmdr)
data("ChickWeight") # load data set

Nachschlagen, was die Spalten in den Daten bedeuten

#?ChickWeight

2 Aufgaben

2.1 Grafiken

Unten ist bereits ein Code, der eine Abbildung erzeugt, die den Zusammenhang zwischen Zeit und Gewicht für jedes Küken abbildet.

  • Kopieren Sie den Code, um eine zweite Abbildung zu erstellen. Ändern Sie den Code dabei so ab, dass der Zusammehang zwischen Zeit und Gewicht für die verschiedenen Diäten sichtbar wird. Verändern Sie den Code auch so, dass die Grafik nur eine Zeile (und 4 Spalten, also eine pro Diät) hat. Färben Sie die Grafik dunkelgrün.

  • Erstellen Sie danach einen Strip-Chart, in dem lediglich die Gewichte der Küken für die verschiedenen Arten von Diät abgebildet werden. Alle vier Strips sollen verschiedene Farben haben. Für die Einfärbung sollen Sie den RColorBrewer verwenden. Wie das geht, können Sie hier lernen. Verwenden Sie für die Grafik die Farbpalette ‘RdBu’.

  • Erstellen Sie einen Plot für Mittelwerte und Konfidenzintervalle (vernachlässigen Sie hier den Faktor “Zeit”).

Weiter unten im Code ist bereits ein Zusat-Plot eingefügt, der mit dem Paket ggplot2 erstellt wurde. Das ist ein Paket, mit dem man (theoretisch) professionelle “Hochglanzabbildungen” erstellen kann. Wenn Sie mehr über das erstellen von Grafiken mit ggplot lernen möchten, schauen Sie sich dieses frei verfügbare Online-Buch an.

2.2 Data-Handling

  • Erstellen Sie ein Subset des Datensatzes, das für die Zeitvariable nur den ersten und den letzten Tag der Messungen enthält.

  • Konvertieren Sie die Zeitvariable in einen Faktor. Vergeben Sie Labels für die Faktorstufen (“first day” und “last day”).

  • Machen Sie dann wieder eine Mittelwertsabbildung, die die Körpergewichte zum ersten und letzten Tag zeigt und dabei auch die vier verschiedenen Diäten berücksichtigt. Ändern Sie den Code für die Grafik so ab, dass Zeit auf der X-Achse ist und “Diet” die verschiedenen Linienfarben sind.

  • Erstellen Sie nun zwei Subsets: Eins mit “first day” und eins mit “last day”.

2.3 Statistische Auswertung

  • Prüfen Sie mittels einfaktorieller ANOVA die Hypothese, dass die Küken in den verschiedenen Diätgruppen an Tag 1 ähnliches Gewicht haben. Levene-Test nicht vergessen.

  • Prüfen Sie mittels einfaktorieller ANOVA die Hypothese, dass die Küken in den verschiedenen Diätgruppen am letzten Tag unterschiedliche Gewichte haben. Levene-Test nicht vergessen.

3 Analysen

3.1 Grafiken

Zusammenhang von Zeit und Gewicht für jedes Küken. Der Code für die Abbildung ist in der R-Beschreibung für den Datensatz zu finden.

require(graphics)
coplot(weight ~ Time | Chick, data = ChickWeight,
       type = "b", show.given = FALSE)

Abändern:

require(graphics)
coplot(weight ~ Time | Diet, data = ChickWeight,
       type = "b", show.given = FALSE, columns = 4, col = "darkgreen")

Strip-Chart:

library("RColorBrewer")
stripchart(weight ~ Diet, vertical=TRUE, method="jitter", jitter = 0.25, xlab="Diet", 
  ylab="weight", col=brewer.pal(n = 4, name = "RdBu"), data=ChickWeight)

Mittelwerte nach Diäten:

with(ChickWeight, plotMeans(weight, Diet, error.bars="conf.int", level=0.95,
   connect=TRUE))

### Zusatzgrafik mit ggplot2

library(ggplot2) # das Paket ggplot2 müssen Sie sich ggf. installieren


myTheme <- theme(plot.title = element_text(face="bold", size = 22),
        axis.title.x = element_text(size = 20),
        axis.title.y = element_text(size = 20),
        axis.text.x = element_text(size = 16, angle = 0), 
        axis.text.y = element_text(size = 16, angle = 0),
        legend.text = element_text(size = 18),
        legend.title = element_text(face = "bold", size = 18),
        strip.text.x = element_text(size = 18),
        #panel.grid.major = element_blank(), 
        #panel.grid.minor = element_blank(), 
        panel.background = element_blank(), 
        panel.grid.major = element_line(color = "#bdbdbd", size = 0.5,linetype = 2),
        axis.line.x = element_line(colour = "black"), 
        axis.line.y = element_line(colour = "black"),
        axis.text = element_text(colour ="black"), 
        axis.ticks = element_line(colour ="black"))

data_graph <- ChickWeight

data_graph$Diet <- factor(data_graph$Diet, levels = c("1", "2", "3", "4"), 
                              labels = c("Diet 1", "Diet 2", "Diet 3", "Diet 4"))

pd <- position_dodge(width = 0.3)

# new labes for the facets 

g <- ggplot(data_graph, aes(x=Time, y=weight, group = Chick)) +
  guides(fill=FALSE)+
  facet_grid( ~ Diet)+
  ggtitle("Trajectory of chick weight over time given \ndifferent diets") +
  scale_y_continuous(limits = c(-5, 400), breaks=seq(0, 400, 50), expand = c(0,0)) +
  geom_line(position = pd, color = "black", size = 1, alpha=0.2) +
  geom_jitter(aes(color = Diet), alpha = 0.3, width = 0.1) +
  stat_summary(aes(y = weight,group=1, fill = Diet), fun.data = mean_cl_boot, geom = "ribbon", width = 0, 
               size = 1, alpha = 0.5) +
  stat_summary(aes(y = weight,group=1, col = Diet), fun.y=mean, 
                geom="line",group=1, size = 1.5, linetype = "solid", alpha = 1)+
  #stat_summary(aes(y = weight,group=1, fill = Diet), fun.y=mean, geom="point", 
  #             color = "black", shape = 22, size = 3, group=1, alpha = 0.8)+
  #stat_summary(aes(y = weight,group=1), fun.y=median, geom="point", 
  #             color = "black", shape = 3, size = 4, group=1, alpha = 1, 
  #             position = position_dodge(width = 0.5))+
  labs(x = "Time (in Days)", y = "Chick weight (g)") +
  scale_color_manual(name = "Diet",values=brewer.pal(n = 4, name = "RdBu"))+
  scale_fill_manual(name = "Diet",values=brewer.pal(n = 4, name = "RdBu"))+
  theme(legend.position = "none")+
  myTheme
## Warning: Ignoring unknown parameters: width
## Warning: `fun.y` is deprecated. Use `fun` instead.
  #theme_bw()
g

#ggsave("results_means.svg",width=9,height=7) # mit diesen Codes können Grafiken gespeichert werden
#ggsave("results_means.pdf",width=9,height=9) # speichern als pdf

3.2 Data-Handling

Daten subsetten und in Faktor konvertieren:

dat_1st_last <- subset(ChickWeight, Time == 0 | Time == 21) # Subset erstellen

dat_1st_last$Time <- factor(dat_1st_last$Time, levels = c(0,21), labels = c("first day", "last day")) # faktorisieren

Abbildung:

with(dat_1st_last, plotMeans(weight, Time, Diet, error.bars="conf.int", 
  level=0.95, connect=TRUE, legend.pos="farright"))

Subsets für first und last day

first_day <- subset(dat_1st_last, Time == "first day")
last_day <- subset(dat_1st_last, Time == "last day")

3.3 Statistische Auswertungen

Mit einfakt. ANOVA testen, ob sich die Gewichte zu Tag 1 zwischen den Diäten unterscheiden. Zuvor aber Leven-Test, ob Varianzen zu Tag 1 gleich sind.

library(mvtnorm, pos=19)
library(survival, pos=19)
library(MASS, pos=19)
library(TH.data, pos=19)
## 
## Attaching package: 'TH.data'
## The following object is masked from 'package:MASS':
## 
##     geyser
library(multcomp, pos=19)
library(abind, pos=24)


Tapply(weight ~ Diet, var, na.action=na.omit, data=first_day) # variances by group
##         1         2         3         4 
## 0.9894737 2.2333333 1.0666667 1.1111111
leveneTest(weight ~ Diet, data=first_day, center="mean") # levene test
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  3  1.9615  0.133
##       46

Test nicht signifikant, also H0, dass Varianzen zu Tag 1 gleich sind, beibehalten. Jetzt kan “normale” ANOVA ohne Welch-Korrektur (für ungleiche Varianzen) durchgeführt werden.

AnovaModel.1 <- aov(weight ~ Diet, data=first_day) # erstellt ANOVA Modell
summary(AnovaModel.1) # Gibt das Ergebnis zurück
##             Df Sum Sq Mean Sq F value Pr(>F)
## Diet         3   4.32   1.440   1.132  0.346
## Residuals   46  58.50   1.272
with(first_day, numSummary(weight, groups=Diet, statistics=c("mean", "sd"))) # Ab hier alles Plot für Mittel-Diffs
##   mean        sd data:n
## 1 41.4 0.9947229     20
## 2 40.7 1.4944341     10
## 3 40.8 1.0327956     10
## 4 41.0 1.0540926     10
local({
  .Pairs <- glht(AnovaModel.1, linfct = mcp(Diet = "Tukey"))
  print(summary(.Pairs)) # pairwise tests
  print(confint(.Pairs, level=0.95)) # confidence intervals
  print(cld(.Pairs, level=0.05)) # compact letter display
  old.oma <- par(oma=c(0, 5, 0, 0))
  plot(confint(.Pairs))
  par(old.oma)
})
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = weight ~ Diet, data = first_day)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0  -0.7000     0.4368  -1.603    0.385
## 3 - 1 == 0  -0.6000     0.4368  -1.374    0.519
## 4 - 1 == 0  -0.4000     0.4368  -0.916    0.795
## 3 - 2 == 0   0.1000     0.5043   0.198    0.997
## 4 - 2 == 0   0.3000     0.5043   0.595    0.932
## 4 - 3 == 0   0.2000     0.5043   0.397    0.978
## (Adjusted p values reported -- single-step method)
## 
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = weight ~ Diet, data = first_day)
## 
## Quantile = 2.6601
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate lwr     upr    
## 2 - 1 == 0 -0.7000  -1.8618  0.4618
## 3 - 1 == 0 -0.6000  -1.7618  0.5618
## 4 - 1 == 0 -0.4000  -1.5618  0.7618
## 3 - 2 == 0  0.1000  -1.2416  1.4416
## 4 - 2 == 0  0.3000  -1.0416  1.6416
## 4 - 3 == 0  0.2000  -1.1416  1.5416
## 
##   1   2   3   4 
## "a" "a" "a" "a"

Die ANOVA mit dem between-subject-Faktor “Diet” ist nicht signifikant, \(F(3,46)~= 1.132\), \(p~= 346\). Die Hypothese, dass sich die Gewichte der Küken zu Tag 1 zwischen den vier verschiedenen Diäten nicht unterscheiden, wird beibehalten.

Jetzt Testen, ob die Gewichte am letzten Tag unterschiedlich sind (aber erste Levene-Test):

Tapply(weight ~ Diet, var, na.action=na.omit, data=last_day) # variances by group
##        1        2        3        4 
## 3445.933 6105.567 5129.789 1879.028
leveneTest(weight ~ Diet, data=last_day, center="mean")
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  3  1.2412 0.3072
##       41

Varianzhomogenität kann angenommen werden.

Jetzt ANOVA:

AnovaModel.2 <- aov(weight ~ Diet, data=last_day)
summary(AnovaModel.2)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## Diet         3  57164   19055   4.655 0.00686 **
## Residuals   41 167839    4094                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
with(last_day, numSummary(weight, groups=Diet, statistics=c("mean", "sd")))
##       mean       sd data:n
## 1 177.7500 58.70207     16
## 2 214.7000 78.13813     10
## 3 270.3000 71.62254     10
## 4 238.5556 43.34775      9
local({
  .Pairs <- glht(AnovaModel.2, linfct = mcp(Diet = "Tukey"))
  print(summary(.Pairs)) # pairwise tests
  print(confint(.Pairs, level=0.95)) # confidence intervals
  print(cld(.Pairs, level=0.05)) # compact letter display
  old.oma <- par(oma=c(0, 5, 0, 0))
  plot(confint(.Pairs))
  par(old.oma)
})
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = weight ~ Diet, data = last_day)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)   
## 2 - 1 == 0    36.95      25.79   1.433   0.4854   
## 3 - 1 == 0    92.55      25.79   3.588   0.0048 **
## 4 - 1 == 0    60.81      26.66   2.281   0.1186   
## 3 - 2 == 0    55.60      28.61   1.943   0.2252   
## 4 - 2 == 0    23.86      29.40   0.811   0.8480   
## 4 - 3 == 0   -31.74      29.40  -1.080   0.7025   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## 
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = weight ~ Diet, data = last_day)
## 
## Quantile = 2.6762
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate  lwr       upr      
## 2 - 1 == 0   36.9500  -32.0748  105.9748
## 3 - 1 == 0   92.5500   23.5252  161.5748
## 4 - 1 == 0   60.8056  -10.5401  132.1512
## 3 - 2 == 0   55.6000  -20.9761  132.1761
## 4 - 2 == 0   23.8556  -54.8190  102.5301
## 4 - 3 == 0  -31.7444 -110.4190   46.9301
## 
##    1    2    3    4 
##  "a" "ab"  "b" "ab"

Effektgröße ausrechnen:

library(lsr)
etaSquared(AnovaModel.2)
##         eta.sq eta.sq.part
## Diet 0.2540591   0.2540591

Die ANOVA ist signifikant, \(F(3,41)~= 4.655\), \(p~= .007\), \(\eta^2~= 0.25\). Die Hypothese, dass sich die Gewichte am letzten Tag zwischen den vier verschiedenen Diäten unterscheiden, wird angenommen. Die post-hoc Tests zeigen, dass der Unterschied durch den Unterschied zw. Diät-1 und Diät-3 entsteht.