# 事前の準備
source("http://www.rforge.net/FSA/InstallFSA.R")
install.packages("FSAdata",repos="http://www.rforge.net/",type="source")


# ライブラリの読み込み
library(FSA)
library(nlstools)
library(FSAdata)

# データの読み込み
data(Croaker2)
Croaker2

# データの表示
Croaker2
head(Croaker2)

#　性別ごとの体長の期待値
tapply(Croaker2$tl,Croaker2$sex,summary)

#　性別ごとに分けた年齢体長関係
plot(tl~age,pch=c(21,19)[Croaker2 $sex],col=c(1,2)[Croaker2$sex],data=Croaker2,main="雌雄別年齢体長プロット")　
legend("topleft",legend=c("F","M"),pch=c(21,19),col=c(1,2))


#　データを性別で分ける
Croaker.M <- subset(Croaker2,sex=="M")
Croaker.F <- subset(Croaker2,sex=="F")


# 初期値の設定　FSAパッケージより
start.full <- vbStarts(tl~age,data=Croaker2,type="typical")    
start.full

# 雌雄まとめたモデル
Bertalanffy.full<-nls(tl~Linf*(1-exp(-K*(age-t0))),data=Croaker2,start=start.full)
overview(Bertalanffy.full)

#雌雄別モデル
Bertalanffy.M <- nls(tl ~ Linf * (1 - exp(-K * (age - t0))), data=Croaker.M, start=start.full)
Bertalanffy.F <- nls(tl ~ Linf * (1 - exp(-K * (age - t0))), data=Croaker.F, start=start.full)

#　プロット
plot(tl ~ age, xlab="年齢", ylab="全長 (mm)", pch=16, col=c(4, 2)[Croaker2$sex]
　,data=Croaker2, main="雌雄別年齢 体長プロット")  
lines(c(1,2,3,4,5,6,7,8,9,10), predict(Bertalanffy.full,
　　new=data.frame(age=c(1,2,3,4,5,6,7,8,9,10))), lwd=2)
lines(c(1,2,3,4,5,6,7,8,9,10), col=2, predict(Bertalanffy.M,
　　new=data.frame(age=c(1,2,3,4,5,6,7,8,9,10))), lwd=2)
lines(c(1,2,3,4,5,6,7,8,9,10), col=4, predict(Bertalanffy.F,
　　new=data.frame(age=c(1,2,3,4,5,6,7,8,9,10))), lwd=2)
legend("bottomright", legend=c("合わせた","雄","雌"), lwd=2, col=c(1,2,4))



# F検定
Smf <- deviance(Bertalanffy.full)
Sf <- deviance(Bertalanffy.M)
Sm <- deviance(Bertalanffy.F)

# F比の計算
F <- (Smf - Sm - Sf)/3  / ( (Sm + Sf) / (length(Croaker2$age) - 3 * 2) )
F

# p値
1 - pf(F, 3, length(Croaker2$age) - 3 * 2)