Nas lagartas, a dieta é mais importante do que o tamanho na resistência aos predadores?

Estou tentando determinar se as lagartas que comem uma dieta natural (monkeyflower) são mais resistentes aos predadores (formigas) do que as lagartas que comem uma dieta artificial (uma mistura de gérmen de trigo e vitaminas). Fiz um estudo experimental com uma amostra pequena (20 lagartas; 10 por dieta). Eu pesava cada lagarta antes do experimento. Ofereci um par de lagartas (uma por dieta) a um grupo de formigas por um período de cinco minutos e contei o número de vezes que cada lagarta foi rejeitada. Repeti esse processo dez vezes.

É assim que meus dados são (A = dieta artificial, N = dieta natural):

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

Estou tentando executar uma ANOVA em R. É assim que meu código se parece (0 = Dieta artificial, 1 = Dieta natural; todos os vetores são organizados com dados para as dez lagartas artificiais da dieta primeiro, seguidos pelos dados da dieta dez natural lagartas):

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)  

E são estes os meus resultados:

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 

Minhas perguntas são:

• A ANOVA é apropriada aqui? Sei que o pequeno tamanho da amostra seria um problema em qualquer teste estatístico; este é apenas um estudo de avaliação no qual eu gostaria de exibir estatísticas para uma apresentação em classe. Eu pretendo refazer este estudo com um tamanho de amostra maior.
• Eu inseri meus dados no R corretamente?
• Isso está me dizendo que a dieta é significativa, mas o peso não?

tl; dr @whuber tem razão em que a dieta e o peso são confundidos em sua análise: é assim que a imagem se parece.

Os pontos de gordura + os intervalos mostram os intervalos de confiança média e de autoinicialização apenas para a dieta; a linha cinza + intervalo de confiança mostra a relação geral com o peso; as linhas individuais + IC mostram as relações com o peso de cada grupo. Há mais rejeição à dieta = N, mas esses indivíduos também têm pesos mais altos.

Indo para os detalhes mecânicos sangrentos: você está no caminho certo com sua análise, mas (1) quando você testa o efeito da dieta, precisa levar em consideração o efeito do peso e vice-versa ; por padrão, R faz uma ANOVA seqüencial , que testa o efeito da dieta sozinha; (2) para dados como esse, você provavelmente deveria estar usando um modelo linear generalizado de Poisson (GLM), embora não faça muita diferença nas conclusões estatísticas neste caso.

Se você observar, em summary()vez de anova()testar os efeitos marginais, verá que nada parece particularmente significativo (você também deve ter cuidado ao testar os principais efeitos na presença de uma interação: neste caso, o efeito da dieta é avaliado peso zero : provavelmente não é sensível, mas como a interação não é significativa (embora tenha um grande efeito!), pode não fazer muita diferença.

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

Centralizando a variável de peso:

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

Não há grandes mudanças na história aqui.

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               

Há algum argumento sobre se os chamados testes "tipo 3" fazem sentido; nem sempre, embora seja útil centralizar o peso. A análise do tipo 2, que testa os principais efeitos após remover a interação do modelo, pode ser mais defensável. Nesse caso, a dieta e o peso são testados na presença um do outro, mas sem as interações incluídas.

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                  

Podemos ver que, se analisássemos a dieta ignorando os efeitos do peso , obteríamos um resultado altamente significativo - isso é essencialmente o que você encontrou em sua análise, devido à ANOVA seqüencial.

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                    

Seria mais padrão ajustar esse tipo de dados a um Poisson GLM ( glm(rejections~diet*cweight,data=dd2,family=poisson)), mas, neste caso, não faz muita diferença para as conclusões.

A propósito, é melhor reorganizar seus dados programaticamente em vez de manualmente, se você puder. Para referência, foi assim que fiz (desculpe pela quantidade de "mágica" que usei):

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"))

Código de plotagem:

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)