在毛毛蟲中，在抵抗捕食者方面，飲食比體型更重要嗎？
我試圖確定吃天然食物（猴花）的毛蟲是否比吃人造食物（小麥胚芽和維生素的混合物）的毛蟲更能抵抗捕食者（螞蟻）。我做了一項小樣本的試驗研究（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 中？
- 這是否告訴我飲食很重要，但體重不重要？
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)
