Return to README.tex CVS log

Up to [local] / ext / README

Initial revision

\documentclass[a4j]{jsarticle}

% typeset: latex, bibtex, latex, latex

\usepackage{natbib} % interface for bibtex
%\bibliographystyle{jplain} %plain, alpha, abbrv, unsrt, jplain, jalpha, jabbrv, junsrt
%\bibliographystyle{plainnat} %plainnat, abbrvnat, unsrtnat for natbib
\bibliographystyle{jecon} %natbibと組み合わせて、日本語（要インストール）
%\bibliographystyle{ecology_letters2.bst}

\usepackage[deluxe, expert, multi,jis2004]{otf} %特にMacでおすすめ

\usepackage[fleqn]{amsmath} % fleqn左揃え
\usepackage{amssymb, mathrsfs, bm}
\usepackage{hyperref} 

%\usepackage[top=30mm, bottom=30mm, left=30mm, right=30mm]{geometry}
%\usepackage[a4paper]{geometry}

\usepackage[dvipdfmx]{graphicx}
\usepackage{wrapfig} % 図の回り込み
\graphicspath{{img/}} % 図のあるディレクトリを指定 ./img/

\usepackage{tabularx} % 表のパッケージ
\usepackage{booktabs} % 表のパッケージ

\begin{document}

\title{絶滅確率計算プログラムの使い方}
\author{箱山　洋\thanks{水産総合研究センター・上田庁舎}}
\date{\today}
\maketitle

\begin{abstract}
生物種の個体群が$t$年後までに絶滅する確率を計算するプログラムの使い方を説明します。このプログラムでは2つのモデルで絶滅リスク評価を行うことができます: 個体数の増加率にこみあい効果なしのモデル(\citealp{lande1988extinction,DENNIS:1991aa})、こみあい効果ありのモデル(\citealp{Hakoyama:2000aa})。モデルについては、原著論文や「水産資源の希少性評価手法について:定量基準E評価の概要」を参照してください。

\end{abstract}

\section{Rのダウンロード}

以下のサイトから辿って、Rをダウンロードします：
{\small
\begin{verbatim}
http://www.r-project.org/
\end{verbatim}}

WindowsやOSXなど、OSに合わせたbinary実行ファイルをダウンロードしてください。

Rを起動すると、図\ref{fig0}のような画面となってプロンプト$>$のあとに入力待ちとなります。

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=35zw]{0.png}
\caption{Rの起動画面}
\label{fig0}
\end{center}
\end{figure}

Rはインタプリンタ言語ですので、プロンプト$>$のあとにコマンドを入力すれば、すぐに結果が出力されます。例えば、$\log{123}$の計算をするには次のように入力します：

{\small
\begin{verbatim}
> log(123)
[1] 4.812184
\end{verbatim}}

\section{データ入力}
\subsection{Excelでデータファイルを作る場合}
個体数の時系列データをCSV形式で保存します。Excelでデータを入力し（図\ref{fig1}）、保存するときに「名前をつけて保存...」からフォーマットをCSVにすれば（図\ref{fig2}）、CSV形式のデータファイルを作ることができます。第一列には年を、第二列には個体数もしくは個体数の相対値を入力して下さい（図\ref{fig1}）。第一行にはヘッダーを入力して下さい（図\ref{fig1}）。ヘッダー名は任意です。ファイルの拡張子は.csvとなります。

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=35zw]{1.png}
\caption{Excelでのデータ入力}
\label{fig1}
\end{center}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=35zw]{2.png}
\caption{ExcelでCSV形式保存}
\label{fig2}
\end{center}
\end{figure}

参考となるテスト用のcsvデータファイルは以下からダウンロードできます：
{\small
\begin{verbatim}
http://cse.fra.affrc.go.jp/hako/extinction/Grizzlies.csv
http://cse.fra.affrc.go.jp/hako/extinction/Funa_biwa.csv
\end{verbatim}}

Grizzlies.csvはイエローストーン国立公園のグレズリーの個体数データ\citep{DENNIS:1991aa}、Funa$\_$biwa.csvは琵琶湖のフナのCPUEデータ\citep{Hakoyama:2000aa}です。

\subsection{CSVファイルのデータをRに読み込ませる}
次に、作成したCSVファイルをRに読み込ませる方法を説明します。まず、Rのコマンドラインにgetwd()と入力して、現在の作業ディレクトリを調べます。例えばOSXの場合、

{\small
\begin{verbatim}
> getwd()
[1] "/Users/hako/Temp/R_workspace"
\end{verbatim}}

などとパスが表示されます。確認した作業ディレクトリに、先に作成したCSV形式のデータファイルを置きます。その上で、次のコマンドを入力してデータを読み込んで下さい：

{\small
\begin{verbatim}
> dat <- read.csv("Grizzlies.csv", header = TRUE)
\end{verbatim}}

``dat''はデータの名前ですので、任意の名前に変更してかまいません。コマンドラインにdatと入力すれば、データがうまく読み込めたかを確認できます：

{\small
\begin{verbatim}
> dat
   Year Population
1  1959         44
2  1960         47
3  1961         46
4  1962         44
5  1963         46
6  1964         45
7  1965         46
8  1966         40
9  1967         39
10 1968         39
11 1969         42
12 1970         44
13 1971         41
14 1972         40
15 1973         33
16 1974         36
17 1975         34
18 1976         39
19 1977         35
20 1978         34
21 1979         38
22 1980         36
23 1981         37
24 1982         41
25 1983         39
26 1984         51
27 1985         47
28 1986         57
29 1987         47
\end{verbatim}}

後述する絶滅計算プログラムにはテストデータがついていますので、データを入力しなくてもテストすることができます。

\subsection{欠測年(不等間隔での測定)がある場合}

時系列データに欠測年(不等間隔での測定)がある場合や、外れ値を除く場合でも、データから絶滅リスクを最尤推定することができます。欠測年がある場合には、その年の個体数にNAを入れるか：

{\small
\begin{verbatim}
> dat
   Year Population
1  2001  400000.00
2  2002         NA
3  2003 1003508.77
\end{verbatim}}

欠測年の行自体を削除して入力してください：

{\small
\begin{verbatim}
> dat
   Year Population
1  2001  400000.00
2  2003 1003508.77
\end{verbatim}}

\subsection{Rに直接データ入力する場合}

直接Rにデータを入力する場合は次のように、コマンドラインに書きます：

{\small
\begin{verbatim}
> dat <- data.frame(Year=c(1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 
1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979,
 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987), 
Population=c(44, 47, 46, 44, 46, 45, 46, 40, 39, 39, 42, 44, 41, 40, 33, 36,
 34, 39, 35, 34, 38, 36, 37, 41, 39, 51, 47, 57, 47))
\end{verbatim}}

テキストファイルに保存しておいて、Rのプロンプトにコピー\&ペーストしてもよいと思います。

\section{絶滅プログラムのロード}
\subsection{インターネットからダウンロードする場合}
コマンドラインに次のように入力します：

{\small
\begin{verbatim}
> source("http://www5.ocn.ne.jp/~hakoyama/hiroshi/extinction_risk_1.R", encoding = "UTF-8")
> source("http://www5.ocn.ne.jp/~hakoyama/hiroshi/extinction_risk_2.R", encoding = "UTF-8")
\end{verbatim}}

extinction\_risk\_1.Rは、個体数の増加率にこみあい密度効果のないモデル(\citealp{lande1988extinction,DENNIS:1991aa})、extinction\_risk\_2.Rはこみあい効果のあるモデル(\citealp{Hakoyama:2000aa})の絶滅確率の計算プログラムです。

これらのファイルを読み込ませることで、絶滅確率を計算する関数ext1とext2が使えるようになります。また、テストデータdatとdat2(dat, dat2：それぞれ、グレズリーとフナのデータ)もロードされます。

\subsection{コンピュータ上にあるプログラムファイルをRに読み込ませる場合}

配布されたextinction\_risk\_1.Rとextinction\_risk\_2.Rがコンピュータに保存してある場合、前述のようにgetwd()で作業ディレクトリを調べて、そのディレクトリ（フォルダ）にext1.Rとext2.Rをコピーします。その上で、次のように入力してプログラムをロードして下さい：

{\small
\begin{verbatim}
> source("extinction_risk_1.R", encoding = "UTF-8")
> source("extinction_risk_2.R", encoding = "UTF-8")
\end{verbatim}}

Rのコンソールにextinction\_risk\_1.Rやextinction\_risk\_2.Rをドラッグ\&ドロップしてもロードできます。

\section{密度効果のないモデル（環境収容力なし）で絶滅確率を計算する場合(extinction\_risk\_1.R)}

環境収容力なしの密度効果のないモデル(\citealp{lande1988extinction,DENNIS:1991aa})で絶滅確率を計算する関数はext1です。減少している個体群にはこちらのモデルが適しています。ext1は次のような引数をとります：

{\small
\begin{verbatim}
Usage:
ext1(dat, t = 100, ne = 1, verbose = FALSE)
\end{verbatim}}

ただし、datはデータファイルの名前、$t$は絶滅確率を計算したい期間（$t$年後までに滅びている確率）、$ne$は絶滅と判断する個体数の閾値、verboseは詳しい結果を表示するかのオプションです。初期値は$t = 100$、$ne = 1$、verbose = FALSEとなっています。初期値のある引数は省略できます。モデルの性質上$ne$の最小値は1です（0を入力しないでください）。

イエローストーンのグレズリー(dat)について、100年後までに個体数が1になる確率は、ext1(dat, t = 100, ne = 1)と入力して求めることができます：

{\small
\begin{verbatim}
> ext1(dat, t = 100, ne = 1)
                                                Estimate
Probability of decline to 1 within 100 years: 9.4128e-05
\end{verbatim}}

$t=10, 20, 100$年でExtinction probabilityを計算すれば（もしくは、3世代、5世代の時間が$t=10, 20$より長ければ、その時間）、基準Eの評価となります。

詳しい出力をしたい場合は、verbose=TRUEとしてext1(dat, t = 100, ne = 1, verbose = TRUE)のように入力します：

{\small
\begin{verbatim}
> ext1(dat, t = 100, ne = 1, verbose = TRUE)
                                                Estimate
Growth rate:                                   0.0023556
Variance:                                        0.01087
Unbiased variance:                              0.011273
Probability attaining the threshold:             0.18849
Conditional mean extinction time:                 1634.4
Probability of decline to 1 within 100 years: 9.4128e-05
\end{verbatim}}

ただし、Growth rate（増加率）とVariance（環境変動の分散）は、パラメーターの最尤推定値、Unbiased varianceは環境変動の分散の不偏推定値、Probability attaining the thresholdは、個体数の閾値(ne)に達する確率（増加率が正の場合、このモデルではサンプルパスが無限に発散し、個体数の閾値に達しないことがある）、Conditional mean extinction timeは、個体数の閾値に達するという条件付きでの平均絶滅時間、Probability of decline to $ne$ within $t$ yearsは、$t$年後までに個体数が閾値$ne$に達する確率を表します。

同様に、100年後までに個体数が10になる確率は、ext1(dat, t = 100, ne = 10, verbose = TRUE)と入力します：

{\small
\begin{verbatim}
> ext1(dat, t = 100, ne = 10, verbose = TRUE)
                                                Estimate
Growth rate:                                   0.0023556
Variance:                                        0.01087
Unbiased variance:                              0.011273
Probability attaining the threshold:             0.51133
Conditional mean extinction time:                 656.96
Probability of decline to 10 within 100 years:  0.096852
\end{verbatim}}

\section{密度効果のあるモデル（環境収容力あり）で絶滅確率を計算する場合(extinction\_risk\_2.R)}

環境収容力ありの密度効果のあるモデル(\citealp{Hakoyama:2000aa})で絶滅確率を計算する関数はext2です。横ばいの個体群にはこちらのモデルが適しています。ext2は次のような引数をとります：

{\small
\begin{verbatim}
Usage:
ext2(dat, gt = 3, t = 100, INDEX = FALSE, PS = 10^6, verbose = FALSE, accuracy = 0.25)
\end{verbatim}}

ただし、datはデータファイルの名前、$gt$は世代時間(年)、$t$は絶滅確率を計算したい期間（$t$年後までに滅びている確率）、INDEXはデータが個体数の相対値の場合にTRUE、個体数そのものの場合にFALSE、$PS$は別の情報から推定した個体数（データが個体数の相対値の場合に絶滅確率を計算するのに必要）、verboseは詳しい結果を表示するかのオプション、accuracyは数値積分関数の精度に関するパラメーターです。初期値は$gt = 3$, $t = 100$, INDEX = FALSE, $PS = 10^6$, verbose = FALSE, accuracy = 0.25となっています。こちらのモデルでは個体数が0になることを絶滅とします。

イエローストーンのグレズリー(dat)について、世代時間を10年としたときの、100年後までに個体数が0になる確率は、つぎのようにして求められます：

{\small
\begin{verbatim}
> ext2(dat, gt = 10, t = 100)
                                           Estimate
Extinction probability within 100 years: 1.3414e-34
> ext2(dat, gt = 10, t = 100, verbose = TRUE)
                                           Estimate
Growth rate per generation:                  4.3849
Carrying capacity:                           41.621
Variance per generation:                    0.12365
Mean extinction time (year):             7.4551e+35
Extinction probability per year:         1.3414e-36
Extinction probability within 100 years: 1.3414e-34
\end{verbatim}}

ただし、Growth rate per generation（世代あたりの増加率）、Carrying capacity（環境収容力）、Variance per generation（世代当りの環境変動の分散）は、パラメーターの近似最尤推定値、Mean extinction time (year)は平均絶滅時間（年）の近似最尤推定値、Extinction probability per yearは一年当りの絶滅確率（平均絶滅時間の逆数）、Extinction probability within $t$ yearsは$t$年後までに個体数が絶滅する確率を表します。

個体数の相対値しか利用できない場合で、別の情報から個体数の推定値PSがある場合、次のようにコマンド入力して絶滅確率を計算できます。ただし、人口学的な確率性が無視できる程に個体数が十分大きいという条件を満たさなければいけません。

{\small
\begin{verbatim}
> ext2(dat2, gt = 3, t = 100, INDEX = TRUE, PS = 10^6)
                                            Estimate
Extinction probability within 100 years: 6.6912e-104
> ext2(dat2, gt = 3, t = 100, INDEX = TRUE, PS = 10^6, verbose = TRUE)
                                            Estimate
Growth rate per generation:                   1.3086
Variance per generation:                     0.11524
Extinction probability within 100 years: 6.6912e-104
\end{verbatim}}

\subsection{平均絶滅時間}

extinction\_risk\_2.Rをロードすると、平均絶滅時間のプログラムも同時に使えるようになります。例えば、増加率$r=0.1$、環境収容力$K=10^5$、環境分散$\sigma_e^2=0.4$、初期個体数$x_0=10^5$の場合の平均絶滅時間は、4906.355年と計算されます：

{\small
\begin{verbatim}
> Te(0.1, 10^5, 0.4, 10^5)
[1] 4906.355
\end{verbatim}}

\bibliography{ref} % ref.bibを引用

\end{document}