No.21240　mvpartによるポアソン回帰木のエラー　　【dive】　2014/08/03(Sun) 13:20

はじめて質問いたします。
魚の生息場所選択性を調べるため，河川の約100地点について魚の個体数と環境変量（水深，流速など7変量）のデータを取り，個体数を応答変数，環境変量を説明変数としてmvpartパッケージのmvpart()を用いて回帰木分析を行っています。
応答変数である魚の個体数はカウント・データであるため，mvpartの引数をmethod="poisson"として分析していますが，以下のようにエラーが出ます。

> t<-mvpart(fish~., data=x, method="poisson", cp=0.001, minsplit=3, xv="min", xvmult=100)
> X-Val rep : 1 2 ... 100
> Minimum tree sizes
> tabmins
> 2 3 4 5 6 7 8 10 13 15 17 18 20 21 23 31
> 1 20 23 13 16 6 4 6 1 3 1 2 1 1 1 1
> 以下にエラー yval[, 1] : 次元数が正しくありません

一方，method="anova"として通常の回帰木分析を行うとエラーは出ません。
また，rpartパッケージのrpart()でmethod="poisson"とすると，エラーは出ません。

後者のパッケージrpartの説明書を読むと，method="poisson"の場合は，反応変数が2列（1列目：観察時間，2列目：事象の回数）もしくは単列（事象の計数値のみ）のベクトルとなっています（http://cran.r-project.org/web/packages/rpart/vignettes/longintro.pdf）。
当方の反応変数（fish）は，魚の個体数のみの単列ベクトルです。
mvpartとrpartでは，method="poisson"とした場合に，反応変数の列数に対する扱い方が異なるようです。

本サイトには，類似の質問が2007年10月（http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc040/04521.html）にありますが，未解決のようです。
また，海外のサイトにも類似の質問を見つけましたが（https://biostat-lists.wustl.edu/sympa/arc/s-news/2001-05/msg00300.html），解答部分がありません

どなたか，mvpartでmethod="poisson"とした場合にエラーが出ないようにする方法をお教えいただけないでしょうか。
よろしくお願いいたします。
No.21244　Re: mvpartによるポアソン回帰木のエラー　　【青木繁伸】　2014/08/03(Sun) 19:06

debug で，実行経緯をトレースすれば何らかの情報が得られるかも知れません。
自分で出来ないということならば，data=x の実際の値を示す。データが大きすぎるとか，秘密にしたいということなら，エラーを再現できる最小限のデータを示す必要があるでしょう。同時に，あなたが実行したプログラムの一部始終。
要するに，第三者がエラーを再現できる状況を明示することが必要ということ。
No.21245　Re: mvpartによるポアソン回帰木のエラー　　【dive】　2014/08/03(Sun) 19:59

青木先生

拙文をお読みくださりありがとうございます。
以下に私が実行したプログラムのすべてを示します。
> library(mvpart)
> x<-read.table("clipboard", header=TRUE)　#エクセルからデータを読み込む
> x #実際には112行*8列（以下に10行*8列を示す）

reach	CU	width	maxd	depth	vero	cover	fish
A	R	17.0 	0.29 	0.06 	58.4 	0.22 	275 
A	GR	19.0 	0.61 	0.08 	29.0 	0.01 	294 
A	RP	13.8 	0.78 	0.10 	2.0 	0.01 	33 
A	R	10.6 	0.61 	0.12 	12.3 	0.00 	50 
A	RP	12.1 	0.61 	0.17 	25.5 	0.20 	10 
B	PG	15.1 	1.57 	0.14 	16.6 	0.68 	14 
B	P	9.5 	2.89 	1.56 	11.2 	0.00 	5 
B	RP	3.6 	2.27 	2.12 	39.3 	0.01 	149 
B	R	4.9 	0.87 	0.76 	96.1 	0.00 	9 
B	R	5.5 	0.56 	0.35 	134.2 	0.00 	0 

> set.seed(0)
> t<-mvpart(fish~., data=x, method="poisson", cp=0, minsplit=3, xv="min", xvmult=100)
X-Val rep : 1  2 ...100
Minimum tree sizes
tabmins
  6 
100 
 以下にエラー yval[, 1] :  次元数が正しくありません
set.seed(0)以降は，上記の10行*8列のデータによる結果です。
試みに，カテゴリカル変数であるreachとCU，また比率であるcoverを除いて分析しなおしてみましたが，結果は同じでした。
また，widthだけを説明変数にしても同じです。
以上，ご指導くださいますようお願いいたします。
No.21246　Re: mvpartによるポアソン回帰木のエラー　　【青木繁伸】　2014/08/04(Mon) 16:10

エラーだと思うので，以下の内容を作者に伝えてあげると，はっきりすると思います。

確信はないですが，たぶんということで取りあえずの対策は以下の通り。

エラーの原因は，mvpart:::text.rpart の中で yval と n を書くときに，method="poisson" の場合には yval というベクトルを渡すにもかかわらず，x$funcitons$text で yval[,1] や yval[,2] のように使おうとするためのエラーです。method="poisson" の場合には，x$functions$text の yval 引数は，yval = data.frame(frame$yval[leaves], frame$n[leaves]) のようにしないとだめだと思います。また，ついでに n 引数は n = frame$n[1] とするのでしょう（たぶん）。

これらのことをちゃんとやるためには，プログラムを書き直す必要があります。
手順は以下の通り。
sink("mvpart2.R")
mvpart
mvpart:::text.rpart
sink()
により保存される mvpart2.R を開き，
1 行目を以下のように変更
mvpart2 <- function (form, data, minauto = TRUE, size, xv = c("1se", "min",
177行目

を削除
新たに177行目となった行を以下のように変更
text.rpart2 <- function (x, splits = TRUE, which = 4, label = "yval", FUN = text,
229行目と231行目を
stat <- x$functions$text(yval = data.frame(frame$yval[leaves], frame$n[leaves]),
ylevel = ylevels, digits = digits, n = frame$n[1],
で置き換え。
128行目の plot.rpart を mvpart:::plot.rpart
201行目の rpartco を mvpart:::rpartco
変更したプログラムを保存。

今後，mvpart を使いたいときは，1度 library(mvpart); source("mvpart2.R") を行った後，mvpart2 関数を呼ぶ。

mvpart2(fish~., data=x, method="poisson", cp=0, minsplit=3, xv="min", xvmult=100)

これにより，以下のような図が描けるはず。

mvpart2.R （エンコーディングは UTF-8）

No.21247　Re: mvpartによるポアソン回帰木のエラー　　【青木繁伸】　2014/08/05(Tue) 01:08

質問に上書きしてしまいました...失礼しました

> 2）先生がご指摘のとおり，プログラムの書き直しについて作者に伝えるべきと思いますが，書き直したプログラムから得られた結果を用いて論文を書く場合，論文中ではどのように書くべきでしょうか。とくに，プログラムの書き直しについて，作者の了承が必要になるでしょうか。

free soft ですから，プログラムを書き直して使用することに作者の了承は不要でしょう。ただ，どの部分をどのように書き直したかの詳細を明示する必要はないでしょうが，そのように書き換えたことが適切かについては，私は責任を持てませんので，「このような理由から，このように直しましたがこれで正しいでしょうか」というような確認はしておく方が無難でしょう。

なお，書き直したプログラムは，図の上にリンクポイントがあります。
> mvpart2.R （エンコーディングは UTF-8）

> 1）Error, CV Error, SE が異なる

確かに，今日やってみると，あなたの結果（この下の図）と同じになりますね。
日付に基づく set.seed をどこかでやっているのかな？

No.21248　Re: mvpartによるポアソン回帰木のエラー　　【dive】　2014/08/05(Tue) 12:25

青木先生

ご回答いただきありがとうございました。

1）Error, CV Error, SE が異なる件

再度，分析いただき重ねてお礼申し上げます。
同じ結果が出たとのことで，安心いたしました。

2）プログラム書き直しの作者への了承

英語でこの件を作者に伝えることができるか，はなはだ不安ですが，連絡を取ってみます。
その際には，作成いただいたリンクを使用させていただきます。
このたびは，大変ありがとうございました。

● 「統計学関連なんでもあり」の過去ログ--- 047 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る