在毛毛虫中，在抵抗捕食者方面，饮食比体型更重要吗？

机器算法验证 r 方差分析
2022-03-20 16:09:16

我试图确定吃天然食物（猴花）的毛虫是否比吃人造食物（小麦胚芽和维生素的混合物）的毛虫更能抵抗捕食者（蚂蚁）。我做了一项小样本的试验研究（20 只毛毛虫；每种饮食 10 只）。我在实验前称过每只毛毛虫的重量。我向一群蚂蚁提供了一对毛毛虫（每种食物一个），持续五分钟，并计算每只毛毛虫被拒绝的次数。我重复了这个过程十次。

这是我的数据的样子（A = 人工饮食，N = 自然饮食）：

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

我正在尝试在 R 中运行 ANOVA。这就是我的代码的样子（0 = 人工饮食，1 = 自然饮食；所有向量首先由十种人工饮食毛虫的数据组织，然后是十种自然饮食的数据毛毛虫）：

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)

这些是我的结果：

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

我的问题是：

  • ANOVA 在这里合适吗？我意识到小样本量对于任何统计测试都是一个问题。这只是一项试验性研究，我想在课堂演示中运行统计数据。我确实计划用更大的样本量重新进行这项研究。
  • 我是否将数据正确输入到 R 中？
  • 这是否告诉我饮食很重要，但体重不重要？
1个回答

tl； @whuber 博士是对的，您的分析中混淆了饮食和体重：这就是图片的样子。

在此处输入图像描述

脂肪点 + 范围显示仅饮食的平均置信区间和引导置信区间；灰线+置信区间表示与权重的整体关系；单独的线 + CI 显示了每组的权重关系。饮食=N 的拒绝更多，但这些人的体重也更高。

进入血腥的机械细节：你的分析是正确的，但是（1）当你测试饮食的影响时，你必须考虑体重的影响，反之亦然；默认情况下，R 执行顺序方差分析，仅测试饮食的影响；(2) 对于这样的数据，您可能应该使用泊松广义线性模型 (GLM)，尽管在这种情况下它对统计结论没有太大影响。

如果您查看测试边际效应的summary()而不是anova()，您会发现没有什么看起来特别重要（在存在交互作用的情况下，您还必须小心测试主效应：在这种情况下，饮食的影响被评估为权重为零：可能不明智，但由于交互不显着（尽管它有很大的影响！）它可能没有太大区别。

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

居中权重变量：

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

这里的故事没有太大的变化。

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

关于所谓的“类型 3”测试是否有意义存在一些争论。它们并不总是如此，尽管将重量居中会有所帮助。将交互作用从模型中剔除后测试主效应的类型 2 分析可能更有说服力。在这种情况下，饮食和体重在彼此存在的情况下进行测试，但不包括相互作用。

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

我们可以看到，如果我们在分析饮食时忽略体重的影响，我们会得到一个非常显着的结果——这基本上就是您在分析中发现的，因为顺序方差分析。

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

将这种数据拟合到 Poisson GLM ( glm(rejections~diet*cweight,data=dd2,family=poisson)) 会更标准，但在这种情况下，它对结论并没有太大的影响。

顺便说一句，如果可以的话，最好以编程方式而不是手动重新排列数据。作为参考，这就是我的做法（对不起我使用的“魔法”数量）：

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

绘图代码：

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)
其它你可能感兴趣的问题
上一篇如何决定将 REML 设置为 True 还是 False？ 下一篇逻辑回归何时不能正常工作？