フィールドデータの解析覚書3
unmarkedと同様にdynamic occupancy modelを解析できるソフトとして、USGS（合衆国地質調査所）のウェブサイトでPRESENCEというソフトが公開されています。SRTM解析のときといいほんとうにUSGSにはお世話になっています。このソフトでは、閉鎖系のとき（つまり季節による変化は想定しない＝single season)と、季節変化も想定したmultiple seasonのときの両方で解析可能です。上記のウェブサイトに”Worked Examples and Exercies in PRESENCE"という練習メモがあります。どこまでも親切です。これとユーザーマニュアルをみながらこのソフトを使えます。
なので、詳しくはそちらを参照したほうがよいのですが、ちょっと引っかかりやすい場所だけメモしておきます。
まず、File＞New projectとして新しい作業を開始します。このとき右下にでてくる"Input Data Form"をクリックするとデータシートがでてきます。Excelデータからコピペできるのですが、まずデータはサンプルの一つNSO_psg209.xlsを示すと下記のようなものです。1997,1998は年（マニュアルなどでは年を季節=seasonとよみかえて表現しています。その下の1,2,3,4,）は調査回です。Territoryは調査地です。あとはいる（1)、いない(0)、と欠測(-)があります。このデータをコピペするときには0,1部分だけを選択してコピーします。
そしてPRESENCEに貼り付けるときに、まず左上のセルを選択し、下記の場所から"Pate values"をつかって貼り付けてください。
そして、下記の#srvy/seasn欄に各季節の調査回数を記入します。これをしないとPRESENCEにデータの階層を理解してもらえません。
ここでFile>save asを押します。"Use last col of data as frequency"というダイアログではNoを選択、これでpaoというファイルをつくらないとそのあとの解析がすすみません。
トップ画面に戻ってOKを押すと、プロジェクトのフォルダが作られます。こうしてファイルを整理してからRus>analyzeと分析を進めてゆきます。
フィールドデータの解析覚書2
今回はロシアの共同研究者から頼まれて、「ある／なし」データの解析をすることになりました。彼も数学が苦手で私に頼んできたのですが、やはり数学な苦手な私になにができるか・・・。とりあえず情報集めは得意なのでとにかく調べまくりました。ある野生動物Aが、いる・いないというデータを彼は10年近く、3か所の調査地で集めていました。各調査地にはさらに３～４の区画があります。各年に１～３回ずつの調査回数です。とにかく体系だっていないデータ！欠測も多い…となるとANOVAなど旧来の方法は使いにくいです。そこで階層ベイズを思いつきました。
ある・なしデータはベルヌーイ分布へのフィットが基本になります。このベルヌーイ分布を主軸としつつ、そこに特定のハビタット要素の影響とか季節性、測定誤差などいろいろなことを考えなければなりません。欠測値を補う必要もあります・・・
あまりの複雑さに泥沼に陥りました。そこで、階層ベイズはとりあえずあきらめて、可能な方法を調べました。するとDynamic occupancy modelというのにあたりました。フィールドでいる、いない（＝ある場所を占有した・しない）データをとったときに、それをフィットさせるというものです（MacKenzie et al.2003)。
いま、調査地 i に、季節 tのあいだにj 回調査におとずれたとき、ある種を確認したかどうかをYijt で表すとします。つまりYijt は1か0です。モデルは次の仮定をおいています。
＊ある季節のなかでの各調査は独立であること
＊調査地 iにその種が実際にいるかいないかの状態Zit（存在状態でもいいますか）は、季節 tの間は
変化しない。
＊実際にはいないのに確認できたというような間違え方（勘違いタイプfalse-positibe errors）はない。
このとき最初の季節（つまりt=1)の調査地iの存在状態Zi1は次のようになります。もちろん0か1です。
Z i1＝Bernoulli (Ψi1)
このときΨ（プサイとよむギリシャ語です）は、存在可能性を示すパラメータであり、これは実際に観測できるデータではありません。同じくZも観測はできません。t=2,3,…としたとき、
Zit～Bernoulli (Zi, t-1×Φit +(1-Zi,t-1)Γit）となります。
ひとつ前の季節t-1にiにその種が存在していなければZ i,t-1=0となります。そしてZitは０×Φit ＋（１－０）×ΓitをパラメータとしたBernoulliとなります。
ΦitのΦはファイ、ΓitのΓはガンマです。それぞれ局所的な生存確率、移住確率です。
このようなバックグラウンドをもとに観察失敗も加味するとある種を確認したかどうかを示すYijtは
Yijt～Bernoulli (Zit×Pijt）となります。
さらにこのモデルに環境変数を加えてその効果をモデリングしてゆくというわけです。この詳細はunmarkedというRのパッケージを参照ください。
（文献）
MacKenzie DI, Nichols JD, Hines JE, Knutson MG, Franklin AB. 2003. Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology. 84:2200–2207
フィールドデータの解析覚書１
（私は統計や数学にとんと弱い研究者です。なので、ざくっとした理解をしてデータ解析をしています。忘れやすい自分への覚書ですので、あまり信用しないで参考程度にしてください。）
かつて野生動物学や生態学の分野では、フィールドでとったデータを解析すると、x-y相関を求めるか有意差（5%とか1%とか）を求めることが多かったです。有意差はカイ二乗検定（データが整数値の場合に使用）とかt検定（データの正規分布を仮定）、またはホイットニーのU検定（これはデータを小さい順とか順位に直して検定）などを利用しました。変数が多くなる場合はANOVA（変数の正規分布を仮定した分散分析）を用いるか、複数の変数を含めて差があるかないかを判断する（クラスカルワリスとか）方法がとられました。しかし、フィールドで得られるデータはだいたいバラつきが大きく、有意差5%というハードルで涙をのんだケースも多かったです。しかし、それは特定の統計モデルに対して、苦労して集めたフィールドデータをあてはめようとしたらあわなかったということです。その統計モデルが適当でなかったということは考えなかったのです。順位に直して検定する方法ではデータのバラつきは考えなくてよかったのですが、逆にせっかくとったデータに含まれている情報を無駄にしているともいえます。
しかし、数年前からより柔軟な解析方法が利用されるようになってきています。まずは一般化線形モデル（GLM）や一般化線形混合モデル（GLMM）です。前者は正規分布以外の分布（たとえば二項分布やベルヌーイ分布）を想定してフィールドデータ（観測データ）をあてはめ、考えられる統計モデルのうち良いモデルを選択するというものです。後者はさらに個体差などもモデルに組み込むことができます。たとえば個体差がランダムに生じるとか、個体差は正規分布に従うといったケースバイケースの柔軟な仮定をモデルにしてゆける便利さがあります。そしてモデルのあてはまりの良さをAICという値で評価して、複数のモデルからもっともよいものを選ぶという方法です。AICはほんとうによく見かけます。
さらに変数が増えて複雑なモデルが必要な場合、階層ベイズ推定が利用されます。ベイズ推定が私には長らくわかりませんでしたが、要は得られた結果から元の状態を推定するということです。フィールドでは「得られた結果」にいろいろな変数や単なる誤差が含まれています。たとえばライトセンサスによるシカの個体数調査なら、サンプリングをしたときの天候とか季節、調査コースが変数になり、誤差は調査者によるものとか、その日のシカの動きによるものとかが考えられます。ベイズ推定では、これらの変数や誤差をひっくるめて観測データができているという前提で解析します。そして、使える情報はなるべく使いたいという発想に立ちます。個体数なら前年の個体数はたぶん影響することでしょう。季節による差も出産期の前後なら、あとのほうが増えると予測が立ちます。それらを推定の柱として、あとはバラつきなどを加味して統計モデルをつくるのです。このとき、さきほどのGLMやGLMMなどをモデルの柱にすると柔軟にデータに対応できるということです。
詳しくは北大の久保さんのウェブサイトや御著書（下記）を参照ください。
Amazon.co.jp： データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学): 久保 拓弥: 本
「パラメータ推定」という言葉が頻繁に出てきますが、これがベイズ推定のミソです。統計モデルをつくるときに変数にかかる係数とか想定した分布の形状を決める係数とか定数項とかがパラメータになります。観測データをもとにモデルをつくり、そのモデルから数千回とか数万回シミュレーション（MCMCとか）を行ってみたとき、それぞれのパラメータが推測できるということです。推測結果は2.5%,5%,95%,97.5%でのパラメータの推定値が示されます。たとえばある係数パラメータが2.5%や5%で負の値、95%や97.5%で正の値ということだと、その係数はあまり影響してないとわかります。逆にいずれも正の値だと推定にプラスの影響を与えていることがわかります。
備忘録・筆ぐるめ18で印刷指定、喪中の人に年賀状を送らない方法
年末になると悩ましい年賀状の季節。数年前から「筆ぐるめ」を使っており、今はバージョン18にしています。しかしこのソフトは少々癖があり、感覚的には操作できないことがあります。今回のテーマの印刷指定、喪中の人に年賀状を送らない方法は毎年、思い出しながら苦労しています。そこで、自分のために方法をメモしておくことにしました。
１）「宛て名」メニューにして、それぞれの人のマーク欄にチェックを入れる。下の例では、マーク１と送付せずにマークを入れている。「送付せず」は自分で入力しt文字、これを「喪中」とかにしてもよい。
２）印刷・メールメニューで下図の「印刷指定チェック」の項を開く。その中の「マークでの印刷指定をチェック」をえらぶ。ここがまずわかりにくい。
３）立ち上がった画面は下図。ここで送付せずを「チェックする」にすると、送付せずにレ点を入れた人たちに送られる。「チェックしない」だとレ点を入れなかった人たちに送られる。これも表現を変えてほしいですね！わかりにくい。
３）印刷・メールメニューで「印刷指定を使用する」をクリック。そのあと、下記のとおり全て印刷のままにしておく。右側にプレビューがあって、そこには出さない人の年賀状も見えるし、枚数も全員分のままなので、どうしても不安になります。でもこういう設定なのです。
３）印刷実行をクリック
上の状態で不安をかかえたまま印刷実行をクリックします。
しかし、すると枚数が指定した枚数になるのでした。
わかりにくいぞ、筆ぐるめ！
※2015年12月追記：住所録を開いたときマーク１～５とか表示項目の設定をするには、拡張タブの中の一覧表示項目設定でチェックを入れます。今回はここで悩みました。
※2020年12月追記：バージョン24では、プリンターで印刷→印刷する宛名を履歴から選択というボタンをクリック→マークのついた人を除外という手順に代わっています。わかりにくいのは相変わらずです。
QGISで世界のネコの分布図をつくる
必要が生じて、世界のネコ科の分布図を描こうと思いました。さて、GIS初心者の私はまた悩みました。そういえばGISでは座標系をあわせるのが基本だったな～。でも全世界なら範囲が広すぎて座標系がいっぱいあります。もしかして、UTMゾーンごとに何十枚と図を作成してから一括してはりつける・・・！？そんなことしていたらたぶん1年かかります。でも、きっとそういうニーズは多いはずだし、GISセミナーで世界地図の上で航空路線を描いたかっこいい図もみたことがありました。そこであれこれ調べたところ、世界全体を図化するときにはロビンソン（Robinson）図法というのがあるようです。もともと丸い地球なので、ゆがみは生じているはずですが、とりあえず図化することはできそうです。次に世界地図を探します。世界の分布図なので、あまり細かい情報は不要です。これにあうものとして、thematicmappingというサイトで、ダウンロード可能でした。シェイプファイル（shp形式）では、simpleとそうでないものがありましたが、今回はシンプルのほうを利用しました。（英語サイトです）
http://thematicmapping.org/downloads/world_borders.php
そして、ネコ科の分布図はIUCNのレッドリストページからダウンロードします。こちらも各種の分布シェイプファイル（shp形式）が提供されています。ログイン登録が必要かもしれません。私は依然にログインidを取得していたので、これはすぐクリア。さらに、各種ごとに使用目的を確認してからの提供となっています。ポップアップ表示されるアンケートに営利目的か・・・などの項目があり、それぞれにYes、Noで答えます。すべてNoだとすぐにダウンロード可能に。ひとつでもYesがあるときは、メールが送信されて数日後にOKかが届きます。出版物での使用などでは注意が必要ですね。それをクリアするとめでたく各種の分布シェイプファイルをゲット。
分布図とともに一つのフォルダにまとめておきます。このとき、フォルダ名を日本語にしないことがポイントです。日本語だとQGISでエラーがでます。
QGIS（QGIS2.2を利用しています）でベクタレイヤを追加と指定して、それぞれのshpファイルを選択して重ねられます。SRSはデフォルトで指定されているように思いますが、必要なら54030（World Robinson）です。
QGISで解析・積雪深マップの作成
北海道内の積雪深マップを作成しようとしました。どの地域が雪深いかを地図上で示すのです。データは国土交通省が提供している国土数値情報の気候データをダウンロードして使用しました。ダウンロードサイトは下記で、統一フォーマット（SHP・GML）のタブから気候値メッシュを選択（ボタンはずっと下のほう）します。
http://nlftp.mlit.go.jp/ksj/gmlold/index.html
このデータには、気温や降水量のデータも含まれています。今回は最深積雪深だけ使えればよいのですが、QGISで解析するときに使用データを絞り込みます。
SHPファイル（シェープファイル、シェイプファイル）とは、各メッシュの測定値がずらっと入っている表のようなものです。詳しくはわかりませんが、拡張子がshpというファイルをシェイプファイルといいつつ、結局はセットでdbfやshxというファイルもついてきます。だからシェイプファイルといっても、結局３つずつのファイルがセットになっているわけです。
もう一つ、ベクタ形式とラスタ形式という言葉もGIS解析のはじめにつまずきました。これらはGISで使用するデータの２つの形式を示しています。ベクタ形式というのは上のシェイプファイルが主ですが、特定の区域にデータが伴っているものです。たとえば今回使用する気候値メッシュは、全国をたくさんのメッシュに分けて、それぞれのメッシュごとに気候データが入ったものです。一方、ラスタ形式では、連続的な変化をもつデータ形式です。衛星画像とか、衛星から観測した標高データなどがあてはまります。GPSで観測したポイントとか、エクセルに入力した観測データはふつう、ベクタ形式のデータとして扱います。そんな風にベクタ形式とラスタ形式の分け方は重要ですが、解析の都合によって両者の形式はQGISの中で互いに変換もできます。
＊気候値メッシュもそうですが、国土数値情報はJDS2000という日本独自の規格（地図作成上の規格）でつくられています。ファイル読み込みの際にはJDS2000のデータとして読み込まなければなりません。
＊北海道全体の気候値データをダウンロードすると30くらいのフォルダにshpファイル、dbf,shx,xmlというファイルができます。このうちshp,dbf,shxというファイルを30セット、つまり90くらいのファイルを一つのフォルダにまとめて入れておきます。QGISの基本として、フォルダ名やファイル名は英小文字と数字のみにしておくと無難です。
＊「ベクタ」→「データマネジメントツール」→「複数のシェープファイルを１つに結合する」と選択します。上記をひとつひとつ開いてあとで結合するよりこのほうが楽です。ブラウズでさきほどのフォルダを選択します。ふつう、選択したら選択したフォルダ内のファイルが画面に表示されると思いますが、なぜか出てきません。でも気にしないで続けます。するとフォルダ名は入力ディレクトリに表示されます。出力シェープファイルは何か名前をつけてください。たとえばsnow.shpとか。shpとつけておきましょう。
＊「マップキャンパスに結果を追加する」の項にはチェックを。これでマップ画面に表示できます。ＯＫを押すと空間参照システムの選択をしなければなりません。JDS2000/UTMzone54です。日本はだいたいこの範囲のようですね。ファイルの数分、これを繰り返し、最後に画面にギザギザの北海道がでてくればまずは成功です。
＊北海道地図の左側に自分でつけたシェープファイル名があるはずです。これを選択して右クリックすると、そのなかに「属性テーブルを開く」という項目があります。「レイヤ」タブ内からでも選択できます。これで情報をみることができます。最深積雪深はG02_058の列です。
＊この時点では北海道は真っ黒です。次に表示させるデータを選ばなければいけません。しかし、ここからトラップが用意されています。「数値情報なのにQGISに来たら、情報はみんな文字列になっているのです！！」これは3日間悩みました。このまま解析をしようとしても、一向に積雪深地図はできません。（だって文字列ですから！）文字列を数値にしなければなりません。
＊この作業にはまず編集モードに切り替えます。属性テーブルを開いたときに画面上の左のほうにある鉛筆アイコンをクリックします。今度は属性テーブル上の右のほうにあるそろばんのアイコンをクリック。これで「フィールド計算機」というのが立ち上がります。新しいフィールドをつくると選択してsnowdepthとか適当な名前をつけます。このsnowdepthというフィールド（つまり列）にG02_058の情報を数値にして入れたいわけです。
＊フィールド計算機は条件タブのcase（条件式をつくる）、変換タブのtoreal（実数にするもの）、フィールドと値タブからG02_058を使って次のような演算式をつくります。実はこのプロセスにもトラップがあって、文字列になった情報の中にはもともと文字列だった情報（たとえば”unknown”があり、それは変換しようとしても数値にできないのです。ですから、それは条件式で処理しなければなりません。よって・・・
CASE WHEN “G02_058” != 'unknown' THEN toreal( “G02_058” ) END
とします。
＊演算が成功して数値に変換されたら、ふたたびsnowdepthを選択して「レイヤ」→「プロパティ」としてシンボルを「段階に分けられた」とすると5段階くらいに分けてくれます。画面上に表示できれば成功です。
ＱＧＩＳによる解析事始め
QGISを使いたい、衛星画像解析をしたい・・・という前に、これらはツールでしかないので、解析方法を吟味したりデータをとる作業は必要です。よく、目的と方法がひっくり返ってしまうことがあるので要注意です。目的が定まらないでデータも不十分だと、結局ただPCに向かっているだけで時間が過ぎてゆきます。あたりまえのことなのですが、これはスタート地点の大きなトラップですので一応記しておきます。「とにかく解析をしたい」ではなく、「何を明らかにしたいか？」を考えましょう。
では、次のステップで、QGISを使うメリットです。最大のメリットはソフトが無料、ライセンスが不要ということでしょう。ほかには、開発プロセスがオープンになっているので、多くの人による改良が日々加えられていることもあげられます。一方、デメリットはまだ完成度の低い部分があり、初心者はトラップにかかりやすいということでしょう。完全なマニュアルがないので、自身の解析に必要な情報をネットで探しながら、答えを探すことになります。専門家と話をしたら数秒で済むことに1か月悩み続けることもありました。ですからあれこれ悩まずに専門家に聞くのも良い（そういうための掲示板もあります）ですが、私自身がそうだったように初心者は質問もしにくいのが実情です。このブログを立ち上げた理由の一つです。
それでは、解析のスタートですが、まず解析用のパソコンとQGISをそろえる必要があります。QGISはGISソフトの中では比較的軽いソフトですが、そこそこのスペックのパソコンでないと、あとの解析でつまります。ネット環境（光回線）も必要です。画像や地図情報をダウンロードしなければならないからです。QGISは
こちらからダウンロードできます。windows版とMac版、Linax版があります。Linax版はわかりませんが、windows版とMac版とではバージョンも挙動も違いますので要注意です。つまずいている原因がソフトのバージョンによる不具合ということがしばしばありますので、そのことも念頭においておきましょう。また、最新バージョンが良いかというと、そうでもありません。往々にして新しいバージョンでは新たな不具合が発生しています。私は最近windows版のQGIS2.2を使いはじめましたが、少し前の1.8と比べつつ、不具合が多ければバージョンを戻そうかと思います。
