# 羅吉斯回歸實例：蟑螂感染率
###### tags: `R教學` `基礎統計`

**練習檔案：**[蟑螂蟯蟲感染](https://drive.google.com/file/d/1sGcJevEKh7cPdlV3KBgAJYrM9ypcy7Cu/view?usp=sharing)

 在蟑螂的實驗中，我們希望了解提供不同食物是否會影響蟑螂體內的蟯蟲感染狀況(感染率)，而感染狀況的變化會不會進一步備蟑螂的發育零期與性別影響。因此我們做了三組實驗，分別使用蟑螂的雌成蟲、雄成蟲，以及若蟲，每組進一步分成兩群，分別給予高蛋白(狗飼料)與高纖維的食物(苜蓿粒)。在一個月後解剖看蟑螂體內的蟯蟲感染狀況。其中感染率的部分藉由羅吉斯回歸做模型，並用ANOVA做檢定。檢定蟑螂的感染率是否受到發育零期與性別的影響。而鑒於蟯蟲的雌蟲、雄蟲以及幼蟲的感染率分別代表不同的生物意義，我們在分析過程中也把不同的蟯蟲齡期分開來計算。

rm(list=ls())
setwd(choose.dir())

##### 數據輸入與資料整理 #####
dat <- read.csv("food.csv", head = T)
dat <- na.omit(dat)
dat$F.ratio <- replace(dat$F, which(dat$F>0), 1)
dat$M.ratio <- replace(dat$M, which(dat$M>0), 1)
dat$L.ratio <- replace(dat$L, which(dat$L>0), 1)

##### 模型建立與統計分析 #####
#雌蟯蟲
F.model <- glm (F.ratio ~ food * stage, 
                data = dat, family = "binomial")
anova(F.model, test = "Chisq")
#雄蟯蟲
M.model <- glm (M.ratio ~ food * stage, 
                data = dat, family = "binomial")
anova(M.model, test = "Chisq")
#幼期蟯蟲
L.model <- glm (L.ratio ~ food * stage, 
                data = dat, family = "binomial")
anova(L.model, test = "Chisq")

##### 繪圖資料整理 #####
Ratio <- aggregate(cbind(F.ratio, M.ratio, L.ratio) ~ food + stage, 
          data = dat, 
          FUN = function(x){sum(x)/length(x)})

##### 繪圖 #####
tiff(filename = "F.infection.tif", res = 300,
     width = 450*7, height = 400*7, pointsize = 20)
plot(Ratio$F.ratio ~ as.factor(rep(c(1,2), 3 )), 
     xlim = c(0.5,2.5), ylim = c(0,1),
     xaxt = "n", las = 1, boxwex = 0.5,
     ylab = "prevalence", xlab = "diet")
points (rep(c(1,2), 3 ), Ratio$F.ratio,
        pch = 19, col = rep(1:3, each = 2))
axis (1, at = 1:2, labels = c("Fiber", "Protein"))
segments (rep(1,3), Ratio$F.ratio[c(1,3,5)],
          rep(2,3),Ratio$F.ratio[c(2,4,6)], 
          lty = 2, col = 1:3, lwd = 2)
text(rep(0.5, 3), Ratio$F.ratio[c(1,3,5)],
     labels = c("Female", "Male", "Nymph"), adj = c(0,0),
     col = 1:3)
dev.off()


##### additional figure #####
tiff(filename = "M.infection.tif", res = 300,
     width = 450*7, height = 400*7, pointsize = 20)
plot(rep(c(1,2), 3 ), Ratio$M.ratio,
     xlim = c(0.8,2.2), ylim = c(0,1),
     xaxt = "n", las = 1,
     ylab = "prevalence", xlab = "diet",
     pch = 19, col = rep(1:3, each = 2))
axis (1, at = 1:2, labels = c("Fiber", "Protein"))
segments (rep(1,3), Ratio$M.ratio[c(1,3,5)],
          rep(2,3),Ratio$M.ratio[c(2,4,6)], 
          lty = 2, col = 1:3, lwd = 2)
text(rep(0.8, 3), Ratio$M.ratio[c(1,3,5)],
     labels = c("Female", "Male", "Nymph"), adj = c(0,0))
dev.off()

tiff(filename = "L.infection.tif", res = 300,
     width = 450*7, height = 400*7, pointsize = 20)
plot(rep(c(1,2), 3 ), Ratio$L.ratio,
     xlim = c(0.8,2.2), ylim = c(0,1),
     xaxt = "n", las = 1,
     ylab = "prevalence", xlab = "diet",
     pch = 19, col = rep(1:3, each = 2))
axis (1, at = 1:2, labels = c("Fiber", "Protein"))
segments (rep(1,3), Ratio$L.ratio[c(1,3,5)],
          rep(2,3),Ratio$L.ratio[c(2,4,6)], 
          lty = 2, col = 1:3, lwd = 2)
text(rep(0.8, 3), Ratio$L.ratio[c(1,3,5)],
     labels = c("Female", "Male", "Nymph"), adj = c(0,0))
dev.off()