# 生長積溫 (Thermal-time)
###### tags: `R教學` `生態統計`

**參考資料**：
1. [Model evaluation](https://cran.r-project.org/web/packages/devRate/vignettes/modelEvaluation.html)
2. [Using devRate package to fit development rate models to an empirical dataset](https://cran.r-project.org/web/packages/devRate/vignettes/quickUserGuide.html)

練習檔案：[真菌溫度發芽率](https://drive.google.com/file/d/16o3_ZwpbqaX2VJy_Au1UcQ3YW4ZFZHQP/view?usp=sharing)

setwd(choose.dir())

require("devRate")

###數據輸入
dat <- read.csv("Bb germination.csv", head = T)
dat$temp <- as.factor(dat$temp)
dat$isolate <- as.factor(dat$isolate)

### Glm 模型檢驗實驗數據
G <- glm(germination~isolate + temp, 
         weight = n, 
         data = dat,
         family = binomial(link = "logit"))

### 因子顯著性比較 (likelihood ratio test)
anova(G, test = "Chisq")

### 生長曲線模擬 
isolate <- levels(dat$isolate)
isolates.Tm <- rep(NA, length(isolate))
isolates.Tm.se <- rep(NA, length(isolate))

for (i in 1:length(isolate)){
  #i = 6
  temp.dat <- dat[dat$isolate == isolate[i], c(2, 3)]
  colnames(temp.dat) <- c("temp", "devRate")
  temp.dat$temp <- as.numeric(as.character(temp.dat$temp))
  
  Temp.model <- devRateModel(eq = taylor_81, 
                             dfData = temp.dat,
                             startValues = list(Rm = 1, Tm = 21, To = 10))
  
  isolates.Tm[i] <- summary(Temp.model)$coefficients[2,1]
  isolates.Tm.se[i] <- summary(Temp.model)$coefficients[2,2]
}

### Tm estimation ###
plot (isolates.Tm, 
      1:length(isolates.Tm), 
      pch = 16, yaxt = "n", 
      ylab = "", xlab = "", 
      xlim = c(20,28))
axis(2, at= 1:10, isolate, las=1)
arrows(isolates.Tm - isolates.Tm.se,
       1:10,
       isolates.Tm + isolates.Tm.se,
       1:10,
       angle = 90, length = 0.05, code = 3, lty = 3)

### curve ###
i = 1
temp.dat <- dat[dat$isolate == isolate[i], c(2, 3)]
colnames(temp.dat) <- c("temp", "devRate")
temp.dat$temp <- as.numeric(as.character(temp.dat$temp))
Temp.model <- devRateModel(eq = taylor_81, 
                           dfData = temp.dat,
                           startValues = list(Rm = 1, Tm = 21, To = 10))
devRatePlot(eq = taylor_81, 
            nlsDR = Temp.model, 
            pch = 16, ylim = c(0, 1.5))