Ist bei Raupen die Ernährung bei der Resistenz gegen Raubtiere wichtiger als die Größe?

Ich versuche festzustellen, ob Raupen, die sich natürlich ernähren (Affenblume), resistenter gegen Raubtiere (Ameisen) sind als Raupen, die sich künstlich ernähren (eine Mischung aus Weizenkeimen und Vitaminen). Ich habe eine Versuchsstudie mit einer kleinen Stichprobengröße (20 Raupen; 10 pro Diät) durchgeführt. Ich habe jede Raupe vor dem Experiment gewogen. Ich bot einer Gruppe von Ameisen für einen Zeitraum von fünf Minuten ein Paar Raupen an (eine pro Diät) und zählte, wie oft jede Raupe abgelehnt wurde. Ich habe diesen Vorgang zehnmal wiederholt.

So sehen meine Daten aus (A = künstliche Ernährung, N = natürliche Ernährung):

Trial A_Weight   N_Weight   A_Rejections N_Rejections
1     0.0496     0.1857     0     1 
2     0.0324     0.1112     0     2
3     0.0291     0.3011     0     2
4     0.0247     0.2066     0     3
5     0.0394     0.1448     3     1
6     0.0641     0.0838     1     3
7     0.0360     0.1963     0     2
8     0.0243     0.145      0     3
9     0.0682     0.1519     0     3
10    0.0225     0.1571     1     0

Ich versuche, eine ANOVA in R auszuführen. So sieht mein Code aus (0 = künstliche Ernährung, 1 = natürliche Ernährung; alle Vektoren werden zuerst mit Daten für die zehn Raupen mit künstlicher Ernährung organisiert, gefolgt von Daten für die zehn natürlichen Diäten Raupen):

diet <- factor (rep (c (0, 1), each = 10) 
rejections <- c(0,0,0,0,3,1,0,0,0,1,1,2,2,3,1,3,2,3,3,0) 
weight <- c(0.0496,0.0324,0.0291,0.0247,0.0394,0.0641,0.036,0.0243,0.0682,0.0225,0.1857,0.1112,0.3011,0.2066,0.1448,0.0838,0.1963,0.145,0.1519,0.1571)   
all.data <- data.frame(Diet=diet, Rejections = rejections, Weight = weight)  
fit.all <- lm(Rejections ~ Diet * Weight, all.data)  
anova(fit.all)  

Und so sehen meine Ergebnisse aus:

Analysis of Variance Table  

Response: Rejections  
            Df  Sum Sq Mean Sq F value   Pr(>F)    
Diet         1 11.2500 11.2500  9.8044 0.006444 ** 
Weight       1  0.0661  0.0661  0.0576 0.813432    
Diet:Weight  1  0.0748  0.0748  0.0652 0.801678    
Residuals   16 18.3591  1.1474                     
--- 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Meine Fragen sind:

• Ist ANOVA hier angemessen? Mir ist klar, dass die geringe Stichprobengröße bei jedem statistischen Test ein Problem darstellen würde. Dies ist nur eine Teststudie, für die ich Statistiken für eine Klassenpräsentation erstellen möchte. Ich habe vor, diese Studie mit einer größeren Stichprobe zu wiederholen.
• Habe ich meine Daten korrekt in R eingegeben?
• Sagt mir das, dass Ernährung wichtig ist, Gewicht aber nicht?

tl; dr @whuber hat recht, dass Ernährung und Gewicht in Ihrer Analyse verwechselt werden: So sieht das Bild aus.

Die Fettpunkte + Bereiche zeigen den Mittelwert und das Bootstrap-Konfidenzintervall nur für die Ernährung. Die graue Linie + Konfidenzintervall zeigt die Gesamtbeziehung zum Gewicht. Die einzelnen Linien + CI zeigen die Beziehungen zum Gewicht für jede Gruppe. Es gibt mehr Ablehnung für Diät = N, aber diese Personen haben auch höhere Gewichte.

Gehen Sie auf die blutigen mechanischen Details ein: Sie sind mit Ihrer Analyse auf dem richtigen Weg, aber (1) wenn Sie die Wirkung der Ernährung testen, müssen Sie die Wirkung des Gewichts berücksichtigen und umgekehrt ; Standardmäßig führt R eine sequentielle ANOVA durch, bei der die Wirkung der Ernährung allein getestet wird. (2) Für solche Daten sollten Sie wahrscheinlich ein verallgemeinertes lineares Poisson-Modell (GLM) verwenden, obwohl es in diesem Fall keinen großen Unterschied zu den statistischen Schlussfolgerungen macht.

Wenn Sie sich summary()eher ansehen als anova(), was marginale Effekte testet, werden Sie feststellen, dass nichts besonders bedeutsam aussieht (Sie müssen auch vorsichtig sein, wenn Sie die Haupteffekte bei Vorhandensein einer Interaktion testen: In diesem Fall wird die Wirkung der Ernährung bei a bewertet Gewicht Null : wahrscheinlich nicht sinnvoll, aber da die Wechselwirkung nicht signifikant ist (obwohl sie einen großen Effekt hat!), macht sie möglicherweise keinen großen Unterschied.

summary(fit.lm <- lm(rejections~diet*weight,data=dd2))
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)    0.3455     0.9119   0.379    0.710
## dietN          1.9557     1.4000   1.397    0.182
## weight         3.9573    21.6920   0.182    0.858
## dietN:weight  -5.7465    22.5013  -0.255    0.802

Zentrieren der Gewichtsvariablen:

dd2$cweight <- dd2$weight-mean(dd2$weight)
summary(fit.clm <- update(fit.lm,rejections~diet*cweight))
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)     0.7559     1.4429   0.524    0.608
## dietN           1.3598     1.5318   0.888    0.388
## cweight         3.9573    21.6920   0.182    0.858
## dietN:cweight  -5.7465    22.5013  -0.255    0.802

Keine großen Änderungen in der Geschichte hier.

car::Anova(fit.clm,type="3")
## Response: rejections
##               Sum Sq Df F value Pr(>F)
## (Intercept)   0.3149  1  0.2744 0.6076
## diet          0.9043  1  0.7881 0.3878
## cweight       0.0382  1  0.0333 0.8575
## diet:cweight  0.0748  1  0.0652 0.8017
## Residuals    18.3591 16               

Es gibt einige Argumente dafür, ob sogenannte "Typ 3" -Tests sinnvoll sind; Sie tun es nicht immer, obwohl das Zentrieren des Gewichts hilft. Eine Typ-2-Analyse, bei der die Haupteffekte nach dem Herausnehmen der Interaktion aus dem Modell getestet werden, ist möglicherweise vertretbarer. In diesem Fall werden Ernährung und Gewicht in Gegenwart voneinander getestet, jedoch ohne die enthaltenen Wechselwirkungen.

car::Anova(fit.clm,type="2")
## Response: rejections
##               Sum Sq Df F value  Pr(>F)  
## diet          4.1179  1  3.5888 0.07639 .
## cweight       0.0661  1  0.0576 0.81343  
## diet:cweight  0.0748  1  0.0652 0.80168  
## Residuals    18.3591 16                  

Wir können sehen, dass wir, wenn wir die Ernährung ohne Berücksichtigung der Auswirkungen des Gewichts analysieren , ein hoch signifikantes Ergebnis erzielen würden - dies ist im Wesentlichen das, was Sie aufgrund der sequentiellen ANOVA in Ihrer Analyse gefunden haben.

fit.lm_diet <- update(fit.clm,. ~ diet)
car::Anova(fit.lm_diet)
## Response: rejections
##           Sum Sq Df F value   Pr(>F)   
## diet       11.25  1  10.946 0.003908 **
## Residuals  18.50 18                    

Es wäre üblicher, diese Art von Daten an ein Poisson GLM ( glm(rejections~diet*cweight,data=dd2,family=poisson)) anzupassen , aber in diesem Fall macht es keinen großen Unterschied zu den Schlussfolgerungen.

Übrigens ist es besser, Ihre Daten programmgesteuert neu zu ordnen, als von Hand, wenn Sie können. Als Referenz habe ich das so gemacht (Entschuldigung für die Menge an "Magie", die ich verwendet habe):

library(tidyr)
library(dplyr)

dd <- read.table(header=TRUE,text=
"Trial A_Weight   N_Weight   A_Rejections N_Rejections
1     0.0496     0.1857     0     1 
2     0.0324     0.1112     0     2
3     0.0291     0.3011     0     2
4     0.0247     0.2066     0     3
5     0.0394     0.1448     3     1
6     0.0641     0.0838     1     3
7     0.0360     0.1963     0     2
8     0.0243     0.145      0     3
9     0.0682     0.1519     0     3
10    0.0225     0.1571     1     0
")

## pick out weight and rearrange to long format
dwt <- dd %>% select(Trial,A_Weight,N_Weight) %>%
    gather(diet,weight,-Trial) %>%
    mutate(diet=gsub("_.*","",diet))
## ditto, rejections
drej <- dd %>% select(Trial,A_Rejections,N_Rejections) %>%
    gather(diet,rejections,-Trial) %>%
    mutate(diet=gsub("_.*","",diet))
## put them back together
dd2 <- full_join(dwt,drej,by=c("Trial","diet"))

Plotcode:

dd_sum <- dd2 %>% group_by(diet) %>%
   do(data.frame(weight=mean(.$weight),
              rbind(mean_cl_boot(.$rejections))))

library(ggplot2); theme_set(theme_bw())
ggplot(dd2,aes(weight,rejections,colour=diet))+
geom_point()+
scale_colour_brewer(palette="Set1")+
scale_fill_brewer(palette="Set1")+
geom_pointrange(data=dd_sum,aes(y=y,ymin=ymin,ymax=ymax),
                size=4,alpha=0.5,show.legend=FALSE)+
geom_smooth(method="lm",aes(fill=diet),alpha=0.1)+
geom_smooth(method="lm",aes(group=1),colour="darkgray",
            alpha=0.1)+
scale_y_continuous(limits=c(0,3),oob=scales::squish)