# irisのデータを使う

library(car)  #Anova関数を使うため

model1 <- lm(Petal.Length ~ Sepal.Length * Species, data=iris)

summary(model1)

Anova(model1, type=3) # 平方和のタイプ選択ができる（1から3）

# 基本関数のanovaでも計算できるが、これはタイプ1の平方和（逐次平方和）なので、変数を入れる順序によって結果が変わる

anova(model1)

# Welchの分散分析

# Welchのt検定（等分散性を仮定しない）は有名（t.test関数でも引数var.equal=Fがデフォルト）だが、

# 一元配置ならoneway.test関数で等分散性を仮定しない分散分析が可能（引数var.equalの設定が可能で、デフォルトがF）

oneway.test(Petal.Length ~ Species, data=iris)

---------------------------------------------

★参考

守屋・広岡（2018）日畜会報 89(1)：1-6  https://doi.org/10.2508/chikusan.89.1

######################################

# Wilcoxonの順位和検定（Mann-WhitneyのU検定と実質同じ）

######################################

# 標準で入っているwilcox.test関数は、順位にタイ（同じ値）がある場合に不正確になる

iris2sp <- iris[1:100,]  # 2群の比較のため、2種のデータを用意する

# coinパッケージを使う方法

library(coin)

wilcox_test(Sepal.Length~Species, data=iris2sp)

# exactRankTestsパッケージを使う方法

library(exactRankTests)

wilcox.exact(Sepal.Length~Species, data=iris2sp)

######################################

# Kruskal-Wallis検定

######################################

# 3群以上の比較。coinパッケージを使う

kruskal_test(Sepal.Length~Species, data=iris)

# 例えば、以下の要素を持つDataを考える
# L：ある生物の計測値（長さとか）
# Exp：実験区（exp）か対照区（ctrl）か
# ID：ある生物の株や巣の番号（同じ株・巣から複数の個体を実験に使い、計測をしている）
# 実験によって計測値が変わるか（実験の効果があるか）を調べる

# IDをランダム効果として入れる

library(lme4)

library(lmerTest) # 係数の有意性検定を実施したい場合は入れておく

Data$ID <- as.factor(Data$ID)

Data$Exp <- as.factor(Data$Exp)

L.model <- lmer(L ~ Exp + (1 | ID), Data)

summary(L.model)

# 尤度比検定（likelihood ratio test: LRT）を実行

L.model.null <- lmer(L ~ (1 | ID), Data)

anova(SCL.model, SCL.model.null)

# lmerのデフォルトはREMLなのに対して、尤度比検定時にはMLで計算する必要があるが、自動で再計算してくれる

# "refitting model(s) with ML (instead of REML)" の注釈が出るはず

# Wald検定

library(car)

Anova(L.model)

# AIC

L.model.ML <- lmer(L ~ Exp + (1 | ID), Data, REML=F) # MLで計算する必要がある

L.model.null.ML <- lmer(L ~ (1 | ID), Data, REML=F)

AIC(L.model.ML, L.model.null.ML)

library(mgcv)

model <- gam(V1 ~ s(V2), data=Data, method="REML")  # 基本

# 周期的な変数を入れる場合（24時間周期、12か月周期など）

# 例えば、変数 Hour が何時のデータか（0～23）を表すとする

model1 <- gam(V1 ~ s(hour, bs="cc"), data=Data, method="REML")   # ccはcyclic cubic regression splineで端が合うようになる

# 場所等のカテゴリカルな変数による違いを追加。カテゴリカルな変数 Loc を考える

# Hour の効果は場所（Loc）によって変わらないことを仮定する場合

model2 <- gam(V1 ~ s(Hour, bs="cc") + Loc, data=Data, method="REML")  

# Hour の効果が場所（Loc）によって変わることを仮定する場合。Hour と Loc の交互作用を考慮するイメージ

model3 <- gam(V1 ~ s(Hour, bs="cc", by=Loc) + Loc, data=Data, method="REML")   

# Hour と連続量V2との交互作用をテンソル積スムーズで考慮

model4 <- gam(V1 ~te(Hour, V2, bs=c("cc", "cr")), data=Data, method="REML")

# 引数methodのデフォルトは"GCV.Cp"だが"REML"または"ML"を使うことが推奨されている（Pedersen et al. 2019）

---------------------------------------------

★参考

Pedersen et al. (2019）PeerJ 7:e6876 https://doi.org/10.7717/peerj.6876

#################################################################

# 生物種の存在データ収集：spoccパッケージを用いてデータベース検索

#################################################################

library(spocc)

library(mapr)

library(scrubr)

spnames <- c("Mauremys japonica") # 二ホンイシガメの場合

# GBIF, iNaturalistを検索

Mj <- occ(query = spnames, from=c("gbif", "inat"), gbifopts=list(country="JP"), inatopts=list(bounds=c(22.5, 122.5, 46.0, 147.0)),limit=500) # GBIFはcountryで国指定が可能。iNaturalistは緯度経度を指定

Mj$gbif # GBIFのレコード

Mj$inat # iNaturalistのレコード

#data frame形式に変換

df_Mj <- occ2df(Mj)

# 緯度経度が文字列として入っているので数値に変換

df_Mj$longitude <- as.numeric(df_Mj$longitude)

df_Mj$latitude <- as.numeric(df_Mj$latitude)

df_Mj <- na.omit(df_Mj)

# データのクリーニング

dframe(df_Mj) %>%

  coord_impossible() %>%

  coord_incomplete() %>%

  coord_unlikely()

# 重複レコードの削除

df_Mj2 <- dframe(df_Mj) %>% dedup()

NROW(df_Mj2) # レコード数確認

# 図示

map_leaflet(df_Mj2)

map_plot(df_Mj2)

# 在情報を出力しておく

write.table(df_Mj2, "Mjaponica_occurence.txt", sep="\t", col.names=TRUE, 

  row.names=FALSE, quote=FALSE, na="NA")

#################################################################

# MaxEntを用いたモデリング

#################################################################

library(gridExtra)

library(rasterVis)

library(raster)

library(rgdal)

library(ENMeval)

## 生物気候データのダウンロードと準備

# GUI上で準備しておくことも可

dir.create("WorldClim_data", showWarnings=F)

# 現在の生物気候：10minのデータ

download.file(url="https://data.biogeo.ucdavis.edu/data/climate/worldclim/1_4/grid/cur/bio_10m_esri.zip",

  destfile="WorldClim_data/current_bioclim_10min.zip", method="auto")

# 解凍

unzip(zipfile="WorldClim_data/current_bioclim_10min.zip", exdir="WorldClim_data/current", overwrite=T)

list.files("WorldClim_data/current/bio")

## 準備終わり

jp_ext <- extent(122, 147, 22, 46)

stk_current <- stack(list.files(path="WorldClim_data/current/bio/", pattern="bio_", full.names=T), RAT=F)

stk_current <- crop(stk_current, jp_ext)

plot(stk_current)

# bio_1, bio_12, bio_8を仮に選択する（レイヤ名とスタックの順番を統一）

selected_vars <- c("bio_1", "bio_12", "bio_8")

stk_select <- stack(subset(stk_current, selected_vars))

plot(stk_select)

# モデリング

occ <- cbind(df_Mj2[,2], df_Mj2[,3])

result1 <- ENMevaluate(occ, stk_select, fc="H", algorithm="maxnet", method="randomkfold", kfolds=10)

# bg.coords: バックグラウンドの緯度経度データフレーム

# method: 検証方法の指定。ENMevalの論文を参照

# algorithm: defaultはRのmaxent使用、java版の場合はmaxent.jar指定

# 結構時間はかかる

result1 

result1@results # 評価指標

# AICcが最小のもので予測プロット

plot(result1@predictions[[which (result1@results$delta.AICc == 0) ]])

points(result1@occ.pts, pch=21, bg=result1@occ.grp)

---------------------------------------------

★MaxEntのモデリングはこの教科書のサンプルコードを参考にしています

「野生生物の生息適地と分布モデリング」Antoine Guisan・Wilfried Thuiller・Niklaus E. Zimmermann 著・久保田 康裕 監訳、共立出版、ISBN 9784320057906