「統計学関連なんでもあり」の過去ログ--- 042

No.08919　glmmML　　【nao】　2009/01/15(Thu) 13:07

glmmML()についてお伺いします。

例データ
y = 1,2,2,3,1
total = 3,3,3,3,3
x = 7,8,9,10,11
id = 1,2,3,4,5

上記のようなデータを，

res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id)
で解析できました。

しかし，上記の例データと同じタイプですが，サンプル数が3000ほどあるデータを下記のようにファイルから読み込むと，

res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id, data=data)

##以下にエラー glmmML.res(X, Y, weights, cluster.weights, start.coef, start.sigma, :
Y values must be 0 <= Y <= 1 ##参照ここまで

とエラーが出ます。

数値を直接代入したときは，yの値が2以上でも計算してくれたのに，なぜファイルを読み込んで計算しようとすると"yが1以下でなければいけない"というエラーが出るのでしょうか？

サンプル数が多いため？なのか，何か別の問題があるのか？
どなたかご存知の方がいれば教えてください。
よろしくお願いします。

No.08920　Re: glmmML　　【青木繁伸】　2009/01/15(Thu) 14:17
range(data$y)
range(data$total-data$y)
を計算したら，どのような結果が得られるか教えてください。

それと，以下のようなデータのときの結果と比べてみるといかがでしょうか？
> data <- data.frame(
+ y = c(1,2,10,3,1),
+ total = c(3,3,3,3,3),
+ x = c(7,8,9,10,11),
+ id = c(1,2,3,4,5)
+ )
> res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id, data=data)
 以下にエラー eval(expr, envir, enclos) : y values must be 0 <= y <= 1
 追加情報:  Warning message:
In log(p.hatt/(1 - p.hatt)) :  計算結果が NaN になりました 
あなたの場合は，追加情報以降は出ていないので，呈示した場合と異なるのかもしれません。
ついでに，RとglmmMLのバージョン，OSも明示して質問するのがよいでしょう。

No.08921　Re: glmmML　　【nao】　2009/01/15(Thu) 15:02

青木さま

使用環境は下記になります。
R 2.8.1
Package glmmML version 0.81-3
Windows XP

計算させたところ，それぞれこのようになりました（=以降に記入）
range(data$y) = 1 3
range(data$total-data$y) = 0 2

> res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id, data=data)
以下にエラー eval(expr, envir, enclos) : y values must be 0 <= y <= 1
追加情報: Warning message:In log(p.hatt/(1 - p.hatt)) : 計算結果が NaN になりました

上記と比べると，ファイルを読み込んだ場合も同じメッセージが出ていますが，追加情報以降は出ませんでした。
glmmMLで二項分布を使う場合は，データフレームにするとダメということなのでしょうか？

よろしくお願いします。

No.08922　Re: glmmML　　【青木繁伸】　2009/01/15(Thu) 15:15

> glmmMLで二項分布を使う場合は，データフレームにするとダメということなのでしょうか？

やってみると判りますが，そんなことはありません。
> data <- data.frame(
+ y = c(1,2,2,3,1),
+ total = c(3,3,3,3,3),
+ x = c(7,8,9,10,11),
+ id = c(1,2,3,4,5)
+ )
> res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id, data=data)
では，何のエラーも出ません。

データを減らしていって，同じ問題が生じる最小のデータセットを特定してはいかがでしょう。データの行数が50以下なら，そのデータを掲示してみると良いかもしれません。

No.08923　Re: glmmML　　【nao】　2009/01/15(Thu) 15:39

アドバイスありがとうございます。
まずは，同じ問題が生じるデータセットを特定してみたいと思います。
ひとつ確認させていただきたいのですが，データの行数が3000ほどあってもRで解析する上で問題はないのでしょうか？
大きなサンプルを扱ったことがなく，かつ小さなサンプルではこれまでにエラーが出ていませんでしたので気になりました。

ネット検索をしてみると，10000以上でも計算しているようなので問題ないかとは思うのですが・・・。

No.08924　Re: glmmML　　【青木繁伸】　2009/01/15(Thu) 16:04

あなたを疑うようで申し訳ないのですが，No. 8921 の結果，

range(data$y) = 1 3
range(data$total-data$y) = 0 2

は，「サンプル数が3000ほどあるデータ」についてのものでしょうか？

> データの行数が3000ほどあってもRで解析する上で問題はないのでしょうか？

普通は問題が出るようなことはありません。もしあるとしても，そのような場合には「データ行数がXXX行以上のデータは処理できません」とかのエラーメッセージを出すのが常識でしょう。

Windows XP, R 2.7.2 で，3000行のデータを作って，やってみましたよ。
> library(glmmML)
> data <- data.frame(y=sample(3, 3000, replace=TRUE), total=rep(3, 3000), x <- rnorm(3000), id=1:3000)
> res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id, data=data)

ちゃんとできました

ここで，データを弄って，data$y[1] を data$total[1] より大きくしてしまいましょう。

> data$y[1] <- 10
> res <- glmmML (cbind (y, total-y)~ x, binomial, cluster = id, data=data)
 以下にエラー eval(expr, envir, enclos) : y values must be 0 <= y <= 1

このエラーメッセージも，あなたが最初に書いたものとは違うのですけどね。。。glmmML のバージョンは，0.81-3　かな。
で，range をとってもらったのは，以下のようなことがないのだろうかと思ったんですが。
> range(data$y)
[1]  1 10
> range(data$total-data$y)
[1] -7 2        マイナスになってしまっている（データ数なのに）

No.08925　Re: glmmML　　【nao】　2009/01/15(Thu) 16:15

青木さま

ありがとうございます。
先生の指摘の通り，データレンジがマイナスになっていました。
お伝えしたのと違ったのは決して意図的ではなかったのですが，同名のファイルで1つはy値が3までで，もうひとつはy値の4以上を削除し忘れていたものです。
中身が同じだと思い，特に考えずに両方で作業していましたが，いただいたアドバイスを見て気づきました。

おかげさまで，エラーメッセージも出なくなりました。
ありがとうございます。