過去のトップページ日記2009.04.17
オオクチバス。ブラックバスというニックネームだけど、あまり黒いって感じではないですね。しかし、アゴは張っているので、確かにオオクチって感じはします。 特定外来生物の代表例みたいなブラックバスだが、日本全国にはバス釣り愛好家が存在し、その生息を望んでいるらしい。実際、ほとんどの河川や湖沼にブラックバスを放ったのはバス釣り愛好家のようだし、はては北海道で養殖していた業者もあるようだ(参考ページ)。 歴史的には、最初は食料として実験的に導入し、次にバス釣りブームにのって各地で勝手にバス釣り愛好家が放流したというのが定説になっている。存在自体が環境破壊なのだが、あまりそうは認識されていないようだ。バサーにとっては河川や湖沼の環境はどうでもいいのだろう。日本釣振興会や、麻生首相を構成員に含む釣魚議員連盟などが「バス釣りの印象が悪くなる」からという理由で、特定外来生物被害防止法の対象魚からブラックバスを外そうと運動していたし、環境よりも自分の趣味や商売という主張を堂々と展開しているのはやや気になる。環境破壊を省みずに違法放流を繰り返したからバス釣りができるわけで、バス釣りの印象は悪くて当然だという認識がまったく感じられない。そういえば、ブラックバス駆除の仕掛けを壊したことを2ch.netで吹聴して報道される事件もあった(参考ページ)。 秋田県の足田堤のような小規模なところでもブラックバスの駆除には大変な費用と手間隙がかかるようだし、琵琶湖に至って言えば底引き網を持ち出して駆除する方向になっている。最近でも違法放流は続いているようなので、完全駆除ができてもいたちごっこになるのだろう。ブラックバスを放流するほうは気軽なのだが、駆除する方は膨大な労力を費やす必要があるのだ。しかしながらバサーは、ブラックバスを駆除する意義も感じていないせいか、この点を認識していない。 こういった状態を改善するには、バサーに費用負担をしてもらえばいいと思う。方法としては、まずバス・フィッシング用の水域と、環境保全の水域を分離する。バス・フィッシング用水域の釣り券に、相当な額の税金をかける。そうやって資金を確保して、環境保全の水域のブラックバス駆除とパトロールをしていけばいいと思う。現実的には、モラルの低いバサーがいる限りブラック・バスを日本から完全に駆除はできないだろうから、管理された水域で、お金を取りつつバス・フィッシングをさせてやるサービスを提供するという算段だ。もちろん、環境保全水域で勝手にバスを放して、勝手にバス・フィッシングをはじめるバサーもいるだろうが、放流禁止地域での釣りに罰則規定を設ければいいと思う。放流は監視できないであろうという技術的な問題と、釣りができなきゃ放流する意欲もわかない事を期待できる。 まぁ、こういう計画を立てるのは簡単なのだが、環境保全派と、バサーの両方の批判を受けることは間違いない。でも環境保全の立場からは、お互いの場所を確保した方が建設的だと思う。釣魚議員連盟ぐらいに、頑張って現実的な制度を作ってもらいたいと思うこの頃です。
画像元ページ:オオクチバス - 水生生物センター
Trackback Ping-URL: http://uncorrelated.no-ip.com/cgi-bin/trackback/20090417
あらしの為のTips(;゚д゚) R言語で質的選択モデルを推計せよ! 住居用の不動産(分譲住宅やマンション)を購入する人は、世の中にたくさんいると思うが、それを何軒も購入する人はほとんどいない。ほとんどの人々にとって、所有する住居は0か1といった選択になる。基本的な計量分析であるOLSでは、説明変数に比例して被説明変数が増加する事を仮定しているので、購入するか否か、つまり1か0かといった『質的選択』で推計を行うと、推計結果にバイアスが入ることになる。 質的選択モデルに関しては、次のような種類に分類されて考えられている。
R言語でこれらの質的選択モデルをどのように推計をするか、順に確認していこう。 1. プロビット・モデル淀・井口(2004)によれば、ブラックバスは1920年〜2000年の間に広まったが、特に1970年代以降はバス・フィッシングのブームに乗って、違法放流者が放流を行い、生息域を拡大したと考えられている。2000年には、富山県で違法放流者が逮捕・起訴される事件が発生しているが、急速な拡大の主要因が違法放流者である明確な証拠は無い。 しかし、ブラックバスは歩行ができないし、流速の速い河川にも生息しないため河川間を移動する可能性も相対的に低い。つまり生息域を拡大には、人為的な操作が必要である。初期のブラックバスの導入は、養殖魚導入されたことが知られているが、その回数は限られている。初期を除いては、養殖用種苗(稚魚)にブラックバスが混じっていたか、違法放流者が放流を行ったかに要因は限られるであろう(もっとも養殖用種苗にオオクチバスが混入していた例は観察されていないらしい)。 ブラックバスの違法放流者の行動を考えよう。ブラックバスがいない内水面(河川や湖沼)に放流すれば、そこでバス・フィッシングを愉しむことができようになる。しかし、ブラックバスを放流するにはバスを捕獲し運送する必要があるので、近い内水面間の違法放流に限られるはずだ。逆に、養殖用種苗の売買は遠隔地間でも、国際取引もあるぐらいなので、ブラックバスの生息域の拡大に地域性は無いと考えられる。 つまり、生息域の拡大パターンに地域性が見られれば、ブラックバスの生息域の拡大には違法放流者が大きな役割を果たしてきたと考える事ができるであろう。 1.1. プロビット・モデル式この仮説を検定するために、以下のプロビット・モデルを設定する。 yは1か0の二項で、1が発生する確率が説明変数xとその係数βで変化する。Φは標準正規分布(小文字が密度関数、大文字が累積分布関数)。 なぜこういうモデルを設定するかというのは、次のような理由である。つまり、y = x'β + εのような線形モデルで被説明変数が1か0かの二項だと、yが1のときにεは1-x'β、yが0のときに-x'βになってしまい、yの分散が(x'β)(1-x'β)となり、誤差項εと説明変数xが相関することになる。つまり、不均一分散になるため、正しく推計できない。 プロビット・モデルは線形結合でないので、OLSなどは使えない。例えば最尤法で推計する必要があり、次のような尤度関数が導出されている(級数ではなくて、行列式で書いた方がRっぽいが、ここは計量のテキストにあわせている。)。
もっとも最尤法はパラメーター初期値に依存する上に計算量が多いので、実際には使われる事は少ない。 1.2. データ・セットデータセットは、淀・井口(2004)のオオクチバスの分布の日本地図を参考に作成した(表「オオクチバスの分布(抜粋)」を参照)。1964年、1969年、1974年の各県の分布から、1969年と1974年の2期を含むサンプルを作成している。なおサンプル中、1期前に既にオオクチバスの生息域になっている県は分析から除外した。違法放流者としては、既にブラックバスがいれば放流する/しないの選択はする必要が無い。
実際のファイルは以下となるので、右クリックで「名前を付けてリンク先を保存」をして利用すること。 1.3. glm()関数による推定R言語でプロビット・モデルを推定する場合、実はnlm()関数ではなくて、glm()関数を使う方が一般的である。どういう仕掛けでGLSで推計できるのかが不明なのだが、推計結果が初期値に依存する可能性がある最尤法よりは扱いやすい。
# ファイルをロード
bbdist <- read.table("オオクチバスの分布.txt", header=TRUE, sep="\t") # 1974年ダミーを作成 年1974 <- ifelse(bbdist$年==1974, 1, 0) # 分布・前期・過去を展開 分布 <- bbdist$分布 前期 <- bbdist$前期 過去 <- bbdist$過去 # 隣県は1以上は1、その他は0にする 隣県 <- ifelse(bbdist$隣県>0, 1, 0) # 前期に既にオオクチバスが生息しているサンプルは除外 年1974 <- 年1974[過去==0] 分布 <- 分布[過去==0] 隣県 <- 隣県[過去==0] # モデル指揮を設定 frml <- 分布 ~ 隣県 + 年1974 # GLSでprobitモデルを推計 pr <- glm(frml, family=binomial(link="probit")) summary(pr)として結果を出力すると、以下のようになる。
> summary(pr)
Call: glm(formula = frml, family = binomial(link = "probit")) Deviance Residuals: Min 1Q Median 3Q Max -1.0687 -0.6766 -0.4992 -0.4992 2.0709 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.1893 0.2607 -4.563 5.05e-06 *** 隣県 0.6619 0.3333 1.986 0.0471 * 年1974 0.3641 0.3338 1.091 0.2753 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 87.163 on 78 degrees of freedom Residual deviance: 80.075 on 76 degrees of freedom AIC: 86.075 Number of Fisher Scoring iterations: 4 隣県の係数が正で5%有意なので、前期に隣県にブラックバスがいれば、その県にブラックバスが侵入してくる可能性が高いことが分かる。つまり、違法放流者の可能性が高いと言えるだろう。 なお、R2は出ないので、たぶん(Residual deviance)/(Null deviance)をして、0.918681を得る必要があるだろうヽ(´д`)ノ 1.4. 5年前に隣県にオオクチバスが生息するときの、オオクチバスの生息確率プロビット・モデルの切片項、説明変数、その係数は、y=1である確率P(y|x'β)を直接表すわけではなく、累積分布関数のパラメーターにしか過ぎない。そこで、P(y|x'β)を解釈したい場合は、以下のような計算をする必要がある。
# 隣県にオオクチバスが生息するときのパラメーター
> bev <- pr$coefficients[1] + pr$coefficients[2]*1 + pr$coefficients[3]*1 # 隣県にオオクチバスが生息しないときのパラメーター > bnv <- pr$coefficients[1] + pr$coefficients[2]*0 + pr$coefficients[3]*1 # 累積標準正規分布の下側確率を求める # 隣県にオオクチバスが生息するときの下側確率 > pnorm(bev) [1] 0.4351411 # 隣県にオオクチバスが生息しないときの下側確率 > pnorm(bnv) [1] 0.2046290 1969年の時点で隣県にオオクチバスがいるときは、1974年の時点で43.5%で新たにオオクチバス生息域になっていた。これは隣県にいなかったときの確率20.0%よりも高い。 2. 二項ロジット・モデル正規分布は、確率密度関数が積分できないので、モデル拡張には向かない。そこで、正規分布の代わりに、正規分布より計算が容易なロジスティック分布を用いた離散モデルが、ロジット・モデルである。ロジスティック分布は正規分布に良く似ているが、裾野が厚いので、やや分散が大きめにでる(参考ページ)。 R言語で実行する場合は、プロビット・モデルと同様にglm()関数でGLS推計を行う。
lr <- glm(frml, family=binomial(link="logit"))
引数familyが変化するだけである。実用上は二値ロジット・モデルはプロビット・モデルがあれば有用とは言えないが、ロジット・モデルは複雑なモデルに拡張が可能なので利用されることが多い。 3. 多項プロビット・モデルMNPパッケージを用いる事で、マルコフ連鎖モンテカルロ法によるベイジアン多項プロビット・モデルを利用することができる。詳細は、Imai, Kosuke and David A. van Dyk. (2005). "MNP: R Package for Fitting the Multinomial Probit Model," Journal of Statistical Software, Vol. 14, No. 3 (May), pp. 1-32.を参照。また、MCMCpackパッケージでも、マルコフ連鎖モンテカルロ法による多項プロビット・モデルを利用可能。 ベイジアンは良く分からないので、利用例など詳細は省略させてくださいヽ(´д`)ノ 4. 多項ロジット・モデルnnetパッケージ、mlogitパッケージを用いることで、容易に利用可能になっている。mlogitパッケージのほうが、より高度な手法を用いることができるので、mlogitの方が利便性が高いだろう。 5. 条件付ロジット・モデル多項ロジット・モデルと同じく、mlogitパッケージを用いることで容易に利用可能になっている。 6. 複合ロジット・モデルまずはmlogitパッケージをインストールしておく。
install.packages("mlogit")
インストール済みの環境ならば、mlogitパッケージをロードする。
library("mlogit")
mlogitパッケージには、Fishというデータがついているので、ヘルプに習い、それを分析してみる。
data("Fishing")
分析するdataframeの形式データセットを、mlogic()関数で使えるように整理する。被説明変数をchoiseで指定しておくことがポイント。なお、varyingが元dataframe内の条件変数の一括指定([条件変数名].[選択肢名]と.で連結された変数名になっている。pr. boatなど。sepオプションで変更可能)、shapeが元dataframeの形状。
Fish <- mlogit.data(Fishing, varying=c(4:11), shape="wide", choice="mode")
推計をまわしてみる。modeは釣り位置、prは価格(何の価格かは不明)、caは釣った数(平均値?)、incomeは所得の模様。prとcaが質的選択によって変化する条件変数(Conditional Variables)、incomeが個別の属性変数である。
mlr <- mlogit(mode ~ pr + ca | income, data=Fish)
summary(mlr) なおモデル式を、mode ~ 1 | incomeとすると、質的選択によって変化する条件変数が無くなるので、多値ロジット・モデルとして分析できる。逆にmode ~ pr + caとすると、条件付ロジット・モデルとして分析になる。 分析結果は以下のように出力される。
Frequencies of alternatives:
beach boat charter pier 0.11337 0.35364 0.38240 0.15059 Newton-Raphson maximisation gradient close to zero. May be a solution 7 iterations, 0h:0m:1s g'(-H)^-1g = 3.39E-27 Coefficients : Estimate Std. Error t-value Pr(>|t|) altboat 5.2728e-01 2.2279e-01 2.3667 0.0179485 * altcharter 1.6944e+00 2.2405e-01 7.5624 3.952e-14 *** altpier 7.7796e-01 2.2049e-01 3.5283 0.0004183 *** pr -2.5117e-02 1.7317e-03 -14.5042 < 2.2e-16 *** ca 3.5778e-01 1.0977e-01 3.2593 0.0011170 ** altboat:income 8.9440e-05 5.0067e-05 1.7864 0.0740345 . altcharter:income -3.3292e-05 5.0341e-05 -0.6613 0.5084032 altpier:income -1.2758e-04 5.0640e-05 -2.5193 0.0117582 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Log-Likelihood: -1215.1 McFadden R^2: 0.81461 Likelihood ratio test : chisq = 565.17 (p.value=< 2.22e-16) 二項関係だと問題にならない条件変数と属性変数の2種類の説明変数だが、条件変数は全ての選択肢で共通の係数を持ち、属性変数は選択肢ごとに係数を持つことが分かる。 所得が高い人は、ボートを好み、桟橋は好まないらしい。また、費用がかかる釣り位置は回避され、よく釣れる釣り位置が好まれる事が分かる。当たり前の結果だが、分析手法の確認としては好例。 7. 順序ロジット・モデルMASSパッケージ(利用例)とDesignパッケージ(利用例)を用いると容易に扱える模様。 おわりにミクロ計量分析でよく使われる質的選択モデルだが、R言語でも手軽に利用することができる事を確認してきた。特にmlogicパッケージを用いると利用は容易であり、質的選択モデルを頻繁に使う消費者理論の実証などでもR言語が使いやすいと言えるだろう。多項プロビット・モデルは、商用統計解析パッケージでも実装されていないものが多数ある新しい分析手法であり、最新の手法で研究が行えるのも魅力になる。 ところで例のデータが取扱説明書の付属のものだと見栄えが悪いので、ミクロの家計データなどで公表されているものをご存知の方は、ぜひ教えてくださいヽ(´д`)ノ 参考ページ
注意
|
過去ログ |