1 はじめに
東京海洋大学資源解析学研究室で作成した， 海ごみ解析パッケージ(malia)内で使用できる関数マニュアルです．
1.1 前提条件
R (>= 3.5.2)のバージョンであること， Rtoolsがインストールされている必要があります．
1.2 インストール
#install.packages("devtools") #devtoolsをインストールしていない場合
#devtools::install_github("y-yasutomo/malia") #パッケージのインストール
library(malia) #呼び出し##
## Attaching package: 'malia'## The following object is masked from 'package:stats':
##
## model.extract##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union一連の解析手順マニュアル簡易版
1.3 参考文献
ライントランセクト法
> Buckland, S.T. et al. (2015) Distance sampling: Methods and applications. Springer International Publishing, Switzerland.
> Buckland, S.T. et al. (2001) Introduction to Distance Sampling. Oxford University Press, Oxford, United Kingdom.
2 解析手順
解析の手順としては次の段階がある．
1. データ整形
2. 推定
3. 結果の統合，作図・作表
解析手順と関数関係:
関数の依存関係
3 データ整形
タブレットデータは解析に用いる形に整形する必要があります.
解析には目視観測データと航海の努力量データが必要となります．それぞれの整形にはSight.hand(),Effort.hand()関数を使います．
3.1 目視観測データ
Sight.hand()関数で目視観測データを解析に用いる形に整形できます.
問題なく実行されると，作業ディレクトリにRev.name.debrisという名称のcsvファイルが作成されます.
#データ(csvファイル)の名前を指定
Data.name<-"まとめ全データ"
#R内での名称(なんでもいい)
Rev.name<-"S19y1"
Sight.hand(Data.name,Rev.name)## Warning in Sight.hand(Data.name, Rev.name): Since S19y1.debris.csv already exists,
## S19y12.debris.csv was created.データの確認
3.2 努力量データ
努力量データにはEffort.hand()関数を使います．
こちらは，作業ディレクトリにRev.name.effortという名称のcsvファイルが作成されます.
## Warning in Effort.hand(Data.name, Rev.name): Since S19y1.effort.csv already exists,
## S19y12.effort.csv was created.データの確認
4 解析
maliaパッケージでは，ライントランセクト法に基づく密度推定を行えます．
発見関数は次の4種類のオプションがあります．
4.1 発見関数
4.1.1 Half normal
\[g(x)=exp(-\frac{x^2}{2\sigma^2})\]
## Warning: package 'ggplot2' was built under R version 3.6.3
Figure 1: Half normal detection function
4.1.2 Hazard rate
\[g(x)=1-exp\left\{ \displaystyle -\left(\frac{x}{\sigma}\right)^{-b} \right\} \quad\]
Figure 2: Hazard rate detection function
4.1.3 Half half normal
\[\begin{equation} g(x) = \begin{cases} exp(-\frac{(x-cutpoint)^2}{2\delta^2}) \quad, & x \leq cutpoint \\ exp(-\frac{(x-cutpoint)^2}{2\sigma^2}) \quad, & x > cutpoint \end{cases} \end{equation}\]
Figure 3: Half half normal detection function
4.1.4 Half hazard rate
\[\begin{equation} g(x) = \begin{cases} exp(-\frac{(x-cutpoint)^2}{2\delta^2}) \quad, & x \leq cutpoint \\ 1-exp\left\{ \displaystyle -\left(\frac{x-cutpoint}{\sigma}\right)^{-b} \right\} \quad, & x > cutpoint \end{cases} \end{equation}\]
Figure 4: Half hazard rate detection function
4.2 発見関数とlegごとの密度推定
推定には，MALIA()関数が使えます.
引数に観測データと努力量データ,
keyで発見関数の型，tdで最大目視距離を与えることで発見関数の推定と，leg,grid,areaごとの密度推定を行います.
発見関数は
key="hn" : Harf normal
key="hr" : Hazard rate
key="hhn" : Harf half normal
key="hhr" : Harf hazard rate
となります.
Sight.Data<-malia::debris
Effort.Data<-malia::effort
td<-200
res<-MALIA(Sight.Data = Sight.Data,
Effort.Data = Effort.Data,
key = "hn",
td = td)## Start from 3973.39950231642
## ite 1 value = 2660.64454140305
## ite 2 value = 2660.64454140305
## generated initial value!
## initial value 2660.644541
## final value 2658.962182
## converged## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました
## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました推定された発見関数の確認
hist(Sight.Data$distance,breaks = c(0, 5, 10, 15, 20, 25,
30, 35, 40, 45, 50, 75, 100, 150, 200),xlim = c(0,70),col="skyblue",main="",xlab = "Distance")
fx<-function(x,para){
integrate(HN,x[1],x[2],para=para)$value/
integrate(HN,0,200,para=para)$value
}
Breaks<-c(0, 5, 10, 15, 20, 25,
30, 35, 40, 45, 50, 75, 100, 150, 200)
ans<-numeric()
for(i in 0:80){
ind<-sum(Breaks<=i)
ans[i+1]<-fx(c(Breaks[ind],Breaks[ind+1]),res$est$est.val)/
integrate(HN,0,200,para=res$est$est.val)$value
}
points(0:80,ans,type="l",lwd=2)
4.3 gridの密度推定
grid.D()関数を使います.
引数には，MALIA()関数の結果オブジェクトに格納されている，legごとの密度推定結果を使います.
## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました
## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました4.4 areaの密度推定
area.D()関数を使います.
引数には，MALIA()関数の結果オブジェクトに格納されている，legごとの密度推定結果を使います.
## $area.density
## Area Density Total.Length
## 1 1 NaN 0.000
## 2 2 NaN 0.000
## 3 3 NaN 0.000
## 4 4 NaN 0.000
## 5 5 NaN 0.000
## 6 6 NaN 0.000
## 7 7 NaN 0.000
## 8 8 NaN 0.000
## 9 9 149.07025 225.183
## 10 10 NaN 0.000
## 11 11 NaN 0.000
## 12 12 40.75428 163.766
## 13 13 NaN 0.000
## 14 14 NaN 0.000
## 15 15 NaN 0.000
## 16 16 NaN 0.000
## 17 17 NaN 0.000
##
## $area.density.leg
## Leg.No. Density Abundance Detected.Number Leg.Length Lat.Start
## 1 1 46.13634 821.9413 25 17.81548 34.22011
## 2 2 80.68412 821.9413 25 10.18715 34.16826
## 3 3 97.16820 1709.6380 52 17.59462 34.12495
## 4 4 48.51237 854.8190 26 17.62064 34.04033
## 5 5 23.47864 361.6542 11 15.40354 33.91543
## 6 6 127.62128 2728.8452 83 21.38237 34.07920
## 7 7 230.07250 3879.5630 118 16.86235 34.61295
## 8 8 234.89176 5391.9351 164 22.95498 34.55104
## 9 9 228.90738 5359.0574 163 23.41147 34.47385
## Lon.Start Lat.End Lon.End area
## 1 137.2803 34.16826 137.0974 12
## 2 137.0974 34.10973 137.0123 12
## 3 136.7982 34.04033 136.6369 12
## 4 136.6369 34.10243 136.4612 12
## 5 134.8557 34.03701 134.7751 9
## 6 134.7292 34.18823 134.9203 9
## 7 135.0398 34.58147 134.8600 9
## 8 134.7499 34.47385 134.5180 9
## 9 134.5180 34.42164 134.2712 9
## [ reached 'max' / getOption("max.print") -- omitted 10 rows ]4.5 一航海分の解析
MALIA()関数では，与えたデータすべてを用いて発見関数を推定する．SDAM()関数は観測データのごみの種類，発見関数，共変量なし・共変量(occo,weather,size)ごと
つまり，一種類のごみに対して16通りの推定を行います．
計算には30-60分程度時間がかかります.
計算が終了すると,作業ディレクトリに航海名.result.objというファイルが作成されます．
以降は，既に作成された結果を用いて結果のハンドリングを行います.
5 モデル選択と結果の抽出
SDAM()関数の結果オブジェクトの読み込み.
aic.summary()関数で，ごみの種類ごとに一覧の結果を返します.
種類EPSの結果.
model.extract()関数を使えば，すべてのごみについてAICによって選択されたモデルを抽出できます. 引数には，aic.summary()関数の結果オブジェクトと航海の解析結果のオブジェクトを与えます．
6 作図
作図用には次の関数があります.
survey.plot() : 航跡図
leg.D.plot() : legごとの密度推定結果
grid.D.plot() : gridごとの密度推定結果
6.1 航跡図
航跡図の作図には，survey.plot()関数が使えます．
引数には，努力量データを与えます．
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
引数のxl,ylを変えれば，描画範囲を選択できます.
デフォルトはxl = c(120, 180),yl = c(10, 50)となっています.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
6.2 leg
leg.D.plot()関数が使えます.
引数には，MALIA()関数の結果からlegごとの推定値の結果を抽出したものを与えます．
## Start from 3973.39950231642
## ite 1 value = 2660.64454140305
## ite 2 value = 2660.64454140305
## generated initial value!
## initial value 2660.644541
## final value 2658.962182
## converged## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました
## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりましたlegごとの密度推定結果
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 5: legごとの密度推定結果
同様に引数のxl,ylを変えれば，描画範囲を選択できます.
デフォルトはxl = c(120, 180),yl = c(10, 50)となっています.
また，引数のTypeに文字列を与えると図中に表示することができます.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 6: legごとの密度推定結果2
6.3 grid
grid.D.plot()関数が使えます.
引数には，grid.D()関数の結果を与えます.
## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました
## Warning in sqrt(E(z)^2 - F(z)): 計算結果が NaN になりました## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 7: gridごとの密度推定結果
引数はleg.D.plot()関数と同様に変更できます.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 8: gridごとの密度推定結果2
7 複数航海の描画
7.1 航跡図
データの読み込み
複数航海を同時にプロットしたい場合，努力量データにvoyage_nameという列を追加します．
Data1<-read.csv("S19y1.effort.csv") %>%
mutate(.,voyage_name = "S19y1")
Data2<-read.csv("U19y1.effort.csv") %>%
mutate(.,voyage_name = "U19y1")データを縦に結合する
プロットには，survey.plot()関数が使えます．この時，引数のmulti=Tとします．(デフォルトではmulti=F)
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
multi=Fとすると，一つの航海として扱われます.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
7.2 leg
複数航海の結果を同時にプロットしたい場合でも，leg.D.plot()関数が使えます.
legごとの推定値の結果を縦に結合し，引数に与えます．
データ用意
#'S19y1で推定'
res1<-MALIA(Sight.Data = read.csv("S19y1.debris.csv"),
Effort.Data = read.csv("S19y1.effort.csv"),
key="hn",td=200)## Start from 8917.63166023761
## ite 1 value = 5597.06850886512
## ite 2 value = 5597.06850886512
## generated initial value!
## initial value 5597.068509
## final value 5556.106719
## converged#'U19y1で推定'
res2<-MALIA(Sight.Data = read.csv("U19y1.debris.csv"),
Effort.Data = read.csv("U19y1.effort.csv"),
key="hn",td=200)## Start from 9245.72635954634
## ite 1 value = 6535.91522660305
## ite 2 value = 6535.91522660305
## generated initial value!
## initial value 6535.915227
## final value 6473.753813
## convergedプロットする
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 9: 複数航海のlegごとの密度推定結果
7.3 grid
複数航海のlegごとの密度推定結果をもとに，gridごとの密度を計算します．計算にはgrid.D()関数が使えます.
この時，先ほどのlegごとの推定値の結果を縦に結合したものを引数に与えます．
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 10: 複数航海のgridごとの密度推定結果
8 作表
作表には
survey.inf()
leg.D.table()
grid.D.table()
area.D.table()
関数が使えます．
8.1 航海の情報
survey.inf()関数は引数に観測データと努力量データ,Artificial,Naturalにそれぞれ人工物・天然物の文字列を与えることで，航海の要約を行えます．
なお，Artificialには"FGN","FGF","FGO","EPS","PBA","PBO","FP","PC","G","M","W","UO"，
Naturalには"SW","DW","NO"がデフォルトで与えられています.
## function (Sight.Data, Effort.Data, Artificial = c("FGN", "FGF",
## "FGO", "EPS", "PBA", "PBO", "FP", "PC", "G", "M", "W", "UO"),
## Natural = c("SW", "DW", "NO"))
## NULLsurvey.info <- survey.inf(Sight.Data = Sight.Data,
Effort.Data = Effort.Data)
survey.info$voyage.inf8.2 leg
調査レグごとの密度推定結果
モデル選択と結果の抽出で用いたmodel.extract()関数の結果オブジェクトを引数に使います．
それぞれのごみについて，選択されたモデルで作表します.
9 人工物・天然物の合算
ls.()関数が使えます．
引数のleg.tableには,leg.D.table()関数の結果オブジェクト，Ar,Ntにそれぞれ人工物・天然物の文字列を与えます. なお，Arには"FGN","FGF","FGO","EPS","PBA","PBO","FP","PC","G","M","W","UO"，
Ntには"SW","DW","NO"がデフォルトで与えられています.
この表には，legごとの人工物・天然物の密度推定結果が格納されていることになるので，これまでのようにlegごとの密度推定結果の作図，gridごとの密度推定結果の作図を行えます．
9.1 人工物
9.1.1 leg
人工物の結果だけを取り出す
Artificial.leg <- An.table %>% select(Leg.No.,Leg.Length,Lat.Start,Lon.Start,Lat.End,Lon.End,Artificial_Density) %>% rename(Density=Artificial_Density)## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 11: legごとの密度推定結果(人工物)
9.1.2 grid
まず，grid.D()関数で，gridごとの密度を計算します．
引数には，人工物のlegごとの密度推定結果が入ったArtificial.legをつかっていることに注意してください．
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 12: gridごとの密度推定結果(人工物)
9.2 天然物
9.2.1 leg
人工物の時と同様の手順です． まず，天然物の結果だけを取り出します．
Natural.leg <- An.table %>% select(Leg.No.,Leg.Length,Lat.Start,Lon.Start,Lat.End,Lon.End,Natural_Density) %>% rename(Density=Natural_Density)## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 13: legごとの密度推定結果(天然物)
9.2.2 grid
grid.D()関数で，gridごとの密度を計算します．
今回，引数には天然物のlegごとの密度推定結果が入ったNatural.legをつかっています．
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.
Figure 14: gridごとの密度推定結果(天然物)
9.2.3 area
## $area.density
## Area Density Total.Length
## 1 1 NaN 0.00000
## 2 2 NaN 0.00000
## 3 3 NaN 0.00000
## 4 4 7.500576 32.70094
## 5 5 50.757654 254.49150
## 6 6 41.979273 183.07380
## 7 7 44.190659 234.70557
## 8 8 NaN 0.00000
## 9 9 147.810139 110.93061
## 10 10 NaN 0.00000
## 11 11 77.834004 104.52435
## 12 12 50.846878 170.78627
## 13 13 NaN 0.00000
## 14 14 NaN 0.00000
## 15 15 NaN 0.00000
## 16 16 NaN 0.00000
## 17 17 NaN 0.00000
##
## $area.density.leg
## Leg.No. Leg.Length Lat.Start Lon.Start Lat.End Lon.End Density area
## 1 1 15.531581 35.04131 139.5642 34.95652 139.4288 64.58159 12
## 2 2 12.434605 34.95652 139.4288 34.89214 139.3174 38.76693 12
## 3 3 23.582111 34.32673 137.5644 34.26223 137.3203 38.95511 12
## 4 4 24.042884 34.26223 137.3203 34.18926 137.0746 31.06039 12
## 5 5 24.422726 34.18926 137.0746 34.12244 136.8223 80.40501 12
## 6 6 11.220152 34.12244 136.8223 34.10081 136.7035 52.56989 12
## 7 7 16.878680 34.07303 136.4176 34.07811 136.2349 28.87378 12
## 8 8 11.821129 33.30119 134.9511 33.25234 134.8383 49.04212 11
## 9 9 16.488111 33.23014 134.7123 33.20872 134.5373 36.96090 11
## 10 10 8.628306 33.20872 134.5373 33.20441 134.4449 62.12303 11
## 11 11 21.778540 33.08532 133.9477 32.98587 133.7467 86.79143 11
## [ reached 'max' / getOption("max.print") -- omitted 53 rows ]## $area.density
## Area Density Total.Length
## 1 1 NaN 0.00000
## 2 2 NaN 0.00000
## 3 3 NaN 0.00000
## 4 4 15.26223 32.70094
## 5 5 86.05811 254.49150
## 6 6 78.91200 183.07380
## 7 7 47.27352 234.70557
## 8 8 NaN 0.00000
## 9 9 108.39523 110.93061
## 10 10 NaN 0.00000
## 11 11 148.49694 104.52435
## 12 12 42.43476 170.78627
## 13 13 NaN 0.00000
## 14 14 NaN 0.00000
## 15 15 NaN 0.00000
## 16 16 NaN 0.00000
## 17 17 NaN 0.00000
##
## $area.density.leg
## Leg.No. Leg.Length Lat.Start Lon.Start Lat.End Lon.End Density
## 1 1 15.531581 35.04131 139.5642 34.95652 139.4288 145.932626
## 2 2 12.434605 34.95652 139.4288 34.89214 139.3174 12.353263
## 3 3 23.582111 34.32673 137.5644 34.26223 137.3203 7.145487
## 4 4 24.042884 34.26223 137.3203 34.18926 137.0746 19.608936
## 5 5 24.422726 34.18926 137.0746 34.12244 136.8223 40.619616
## 6 6 11.220152 34.12244 136.8223 34.10081 136.7035 18.603033
## 7 7 16.878680 34.07303 136.4176 34.07811 136.2349 49.057632
## 8 8 11.821129 33.30119 134.9511 33.25234 134.8383 7.577291
## 9 9 16.488111 33.23014 134.7123 33.20872 134.5373 127.706371
## 10 10 8.628306 33.20872 134.5373 33.20441 134.4449 118.338281
## 11 11 21.778540 33.08532 133.9477 32.98587 133.7467 117.082519
## area
## 1 12
## 2 12
## 3 12
## 4 12
## 5 12
## 6 12
## 7 12
## 8 11
## 9 11
## 10 11
## 11 11
## [ reached 'max' / getOption("max.print") -- omitted 53 rows ]