104089
統計学関連なんでもあり

統計学に関する質問とか意見などなんでも書き込んでください。回答があるか,回答が正しいかは保証できません。
レポートや宿題の答えを求める記事や,その他不適切と判断した記事は,管理者の独断で即刻削除します。
ローマ数字,丸付き数字などのフォントセット依存文字および半角カタカナは使わないでください。
記事を引用する際には,適切に引用してください(全文引用はしないでください)。
問題が解決した(結局解決しなかった)場合は,その旨報告していただければ,コメントした人や他の読者に参考になるでしょう。


[トップに戻る] [利用上の注意] [ワード検索] [過去ログ] [統計学自習ノート] [管理用]
latest article: No. 22238, 2016/12/09(Fri) 17:36
おなまえ
タイトル
コメント
URL
添付ファイル
暗証キー (英数字で8文字以内)

自分の投稿記事を修正・削除する場合 --> 処理 記事No 暗証キー


因子分析の信頼性と妥当性について
  投稿者:さくら 2016/12/07(Wed) 19:01 No. 22237
先生の掲示板でいつも色々と学ばせて頂いております。

これまでは因子分析で抽出された因子の信頼性評価をクロンバックαに
より求めてきました。

先行研究に、合成信頼性を使った信頼性評価と、平均分散抽出に
よる妥当性が必要と記載の上、検証されています。現在、アンケートに
より回収したデータを用いて因子分析を行っている最中です。

平均分散抽出とはどのような意味で、どのように値を求めたらよいので
しょうか。また合成信頼性についても意味するところと値の求め方が
分からずに困っております。

今後、このあたりについてどのような文献を用いて学んでいったら
良いのか、目星がつけられません。

ご教示頂ければ幸いです。よろしくお願い申し上げます。

Re: 因子分析の信頼性と妥当性について
  投稿者:鈴木康弘 2016/12/09(Fri) 17:36 No. 22238
 日本語のサイトではこれというのが見つかりませんね。
"average variance extracted" と "composite reliability" で探すと

https://www.researchgate.net/post/Can_anyone_tell_me_how_to_calculate_average_variance_extracted_AVE_and_composite_reliability_CR_of_a_single_latent_variable_with_5_indicators

 というのが当たります。英語ですけど。


R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:明石 2016/12/01(Thu) 22:10 No. 22227
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。

何卒どうぞ、よろしくお願いいたします。

ーーー

データフレームの行方向での結合 rbind()関数を適用できる前提条件は、
カラム数・カラム名の並びが同一であることは理解しております。

カラム数・カラム名が異なる複数ファイルの、行方向での結合について
ご教示をいただければ大変に助かります。

結合対象のファイル名は、fllenames[n]のベクトルに格納されているとします。

素人考えですが、以下のようにやればよいのかと、思いを巡らせている段階です。

・まず、fllenames[1]のファイル名で読み込み、データフレームdat1に格納します。

・次に、fllenames[2]のファイル名で読み込み、データフレームdat2に格納します。

・dat1のカラム名と、dat2のカラム名とを比較して、
 rbindができるように、dat1とdat2のカラム名を整形してあわせる(※)。

・これを、逐次的に、fllenames[n]まで、すべてのファイル名について行う

言葉ではさらりと書きましたが、※ の箇所の具体的な手順を書き出すことができません。

もっと良い方法があるように思います。

もし可能ならば、具体的なRプログラムをお示ししていただき、
ご教示をいただければ大変に助かります。

何卒どうぞ、よろしくお願い申し上げます。

失礼いたします。

Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:青木繁伸 2016/12/02(Fri) 09:05 No. 22228
具体的によくわからないけど,merge でできそうに思いますが?

Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:明石 2016/12/02(Fri) 09:56 No. 22229
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

早々にご教示をいただきまして、誠にありがとうございました。

「転置させて、カラム名をキーにしてmergeをする」という意味と理解しました。

これだけでも貴重な手がかりです。
これを手がかりにして、検討をしてみます。

分かりにくい、つまらない質問で、大変に申し訳ございませんでした。
以降、気をつけます。
失礼いたします。

Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:青木繁伸 2016/12/02(Fri) 11:15 No. 22230
具体例を示してくれればよいのだけど?
merge は(大部分が)同じ個体についてのデータフレームをマージするのが目的。
「カラム数・カラム名が異なる複数ファイル」というのは,全く別の個体のデータですか(普通はそんなのは考えられないけど)。
もしカラム名の重複がないデータなら,両者に同じカラム名(以下では id)で値の重複もないカラムを追加して,all=TRUE で merge すればよい。
a = data.frame(id=1:5, x=1:5*3, y=(1:5)^2)
a
b = data.frame(id=11:20, z=runif(10))
b
merge(a, b, all=TRUE)
実行例
> a = data.frame(id=1:5, x=1:5*3, y=(1:5)^2)
> a
id x y
1 1 3 1
2 2 6 4
3 3 9 9
4 4 12 16
5 5 15 25
> b = data.frame(id=11:20, z=runif(10))
> b
id z
1 11 0.45929778
2 12 0.03272818
3 13 0.51413244
4 14 0.07282785
5 15 0.78663924
6 16 0.12525008
7 17 0.87019710
8 18 0.37075708
9 19 0.13116714
10 20 0.44714621
> merge(a, b, all=TRUE)
id x y z
1 1 3 1 NA
2 2 6 4 NA
3 3 9 9 NA
4 4 12 16 NA
5 5 15 25 NA
6 11 NA NA 0.45929778
7 12 NA NA 0.03272818
8 13 NA NA 0.51413244
9 14 NA NA 0.07282785
10 15 NA NA 0.78663924
11 16 NA NA 0.12525008
12 17 NA NA 0.87019710
13 18 NA NA 0.37075708
14 19 NA NA 0.13116714
15 20 NA NA 0.44714621


Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:明石 2016/12/03(Sat) 09:45 No. 22231
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

ご親身にご相談にのってくださり、誠にありがとうございます。
改めて御礼を申し上げます。

やりたいことをどのようにお伝えすればよいのか、色々と考えており、
ご返信が遅くなり、申し訳ございませんでした。

青木先生がお書きになられたことは、よく理解できました。
ありがとうございます。

私がやりたいと思うことは、
個体〜観測変数 というデータの結合ではなく、
データベースから出力した帳票類のデータファイルの結合です。

入力ファイルが3つからなる簡単な例示を、添付画像でお示しをします。
ファイル1の項目は、v1、v2、v3
ファイル2の項目は、v1、v2、v4
ファイル3の項目は、v2、v3、v5

これら3つのファイルを行方向に結合したいのですが、
カラム名が一致しないために、単純にはrbindが適用できません。これが悩みです。

共通する項目については、それを利用し、
共通しない項目については、右側の列に追加しながら、
行方向に結合します。

データがない箇所は、NAで埋めます。

結果として、
行方向の数は、入力ファイルの行数の和
列方向の数は、入力ファイルのカラム名の和集合
となります。

これがやりたいことです。

先の投稿では、説明不足で申し訳ございませんでした。
どうぞ、よろしくお願い申し上げます。


Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:青木繁伸 2016/12/03(Sat) 13:18 No. 22232
なので,merge を使えばできるでしょう?
> file1 = data.frame(v1=rep(1,5), v2=rep(1,5), v3=rep(1,5))
> file2 = data.frame(v1=rep(2,5), v2=rep(2,5), v4=rep(2,5))
> file3 = data.frame(v2=rep(3,5), v3=rep(3,5), v5=rep(3,5))

> file1$id = 1:5 # ユニークな値を取る変数を付加する
> file2$id = 1:5 + 10 # 不要な場合もあるが
> file3$id = 1:5 + 20 # 間違いを防ぐためには必須

> file1
v1 v2 v3 id
1 1 1 1 1
2 1 1 1 2
3 1 1 1 3
4 1 1 1 4
5 1 1 1 5
> file2
v1 v2 v4 id
1 2 2 2 11
2 2 2 2 12
3 2 2 2 13
4 2 2 2 14
5 2 2 2 15
> file3
v2 v3 v5 id
1 3 3 3 21
2 3 3 3 22
3 3 3 3 23
4 3 3 3 24
5 3 3 3 25

> a = merge(merge(file1, file2, all=TRUE), file3, all=TRUE)
> a = a[,c("v1", "v2", "v3", "v4", "v5")] # 列順を決める必要があるならこれを追加
> a
v1 v2 v3 v4 v5
1 1 1 1 NA NA
2 1 1 1 NA NA
3 1 1 1 NA NA
4 1 1 1 NA NA
5 1 1 1 NA NA
6 2 2 NA 2 NA
7 2 2 NA 2 NA
8 2 2 NA 2 NA
9 2 2 NA 2 NA
10 2 2 NA 2 NA
11 NA 3 3 NA 3
12 NA 3 3 NA 3
13 NA 3 3 NA 3
14 NA 3 3 NA 3
15 NA 3 3 NA 3
複数のデータフレームのマージは,
a = read.table(filenames[1])
for (i in 2:n) {
a = merge(a, read.table(filenames[i]), all=TRUE)
}
のようにやる。

【御礼】 Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:明石 2016/12/03(Sat) 17:02 No. 22233
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

青木先生が「mergeでできる」とおっしゃった意味を、納得し、理解できました。

今回も助けていただきました。
心から御礼を申し上げます。
ありがとうほざいました。

失礼いたします。

Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:明石 2016/12/04(Sun) 08:44 No. 22234
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

おさらいをしていて、追加質問です。
どうぞ、よろしくお願いいたします。

私は、データベース操作言語のSQLを利用していますので、
Rのmergeも、その理解でおりました。

結合の種類(内部、左外部、右外部、完全外部)に応じた指定方法がありますが、
結合キーは必須だと思っていました。

青木先生からご教示いただきました方法では、結合キーの指定がないことに、
驚きました。

Rのmergeは、SQLの結合だと思っていた私には、
このような使い方があることを、まったく知りませんでした。

追加質問です。

> file1$id = 1:5 # ユニークな値を取る変数を付加する
> file2$id = 1:5 + 10 # 不要な場合もあるが
> file3$id = 1:5 + 20 # 間違いを防ぐためには必須

idを付加する場合、しない場合について、
私なりに色々と実験をしましたが、この意味合いが理解できません。

ご説明をいただけますと、大変に助かります。
どうぞ、よろしくお願いいたします。

色々とありがとうございます。

Re: R カラム数・カラム名が異なる複数ファイルの、行方向での結合
  投稿者:青木繁伸 2016/12/04(Sun) 11:18 No. 22235
結合キーは,R が自動で指定します
両者に共通の列です
今のような場合には複数のファイルを指定することもあり,id 列を付加することで,面倒な処理を避けることができるでしょう
id 列を付加しないと,思いがけない結果になりびっくりすることになります
以下の例を吟味するとよいでしょう(file2 が前の例と異なります)
file1 = data.frame(v1=rep(1,5), v2=rep(1,5), v3=rep(1,5))
# file2 = data.frame(v1=rep(2,5), v2=rep(2,5), v4=rep(2,5))
file2 = data.frame(v1=rep(1,5), v2=rep(1,5), v4=rep(2,5))
file3 = data.frame(v2=rep(3,5), v3=rep(3,5), v5=rep(3,5))
file1
file2
file3
a = merge(merge(file1, file2, all=TRUE), file3, all=TRUE)
a = a[,c("v1", "v2", "v3", "v4", "v5")] # 列順を決める必要があるならこれを追加
a

Special Thanks(R カラム数・カラム名が異なる複数ファイルの、行方向での結合)
  投稿者:明石 2016/12/04(Sun) 12:13 No. 22236
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

何度も何度もご丁寧にご教示をいただきまして、誠にありがとうございました。

おかげさまで理解できました。

新しい追加拡張パッケージが沢山ありますが、
baseの基本機能を正しく理解することで、対応できることが多いなと、
改めて、Rはよく設計されていると感心しました。

心から御礼を申し上げます。
失礼いたします。


カイ2乗検定について
  投稿者:心理系院生 2016/11/30(Wed) 12:56 No. 22224
2群の割合に違いがあるかを検討するために、カイ2乗検定を行うのがよいと考えました。
しかし、2群の標本数が大きく異なっており、このような場合にカイ2乗検定を行ってよいものかわかりませんでした。
一方は10000対7000ほどで、もう一方は60対10程度です。
このような違いがあってもカイ2乗検定を行ってもよいのでしょうか。
もし他に適切な手法があれば、ご教授いただけましたら幸いです。
お忙しい所恐れ入りますが、何卒よろしくお願いいたします。

Re: カイ2乗検定について
  投稿者:かい 2016/12/01(Thu) 13:14 No. 22225
サンプルサイズが異なっても検定に問題ありませんが,10000対7000の方は十分精度が高い(区間推定の範囲が非常に狭い)ので,実質的には10000対7000の方の点推定値と60対10の方の比率を比較しているのとあまり変わらない結果になります.

Re: カイ2乗検定について
  投稿者:心理系院生 2016/12/01(Thu) 14:54 No. 22226
かい様

お返事ありがとうございます。
よく理解できました。
問題はないということで安心しました。
このまま研究を進めていきたいと思います。


尺度開発
  投稿者:kashio 2016/11/15(Tue) 14:19 No. 22211
コミュニケーションに関する尺度を開発中です。
いくつかの書籍を参考に因子分析をしました。累積寄与率を示したかったし、スクリュープロットで4因子にかくていし、下位尺度や信頼性の結果も納得できたので投稿したところ査読者から以下のコメントをもらいました。今時 が引っかかり、納得も行かず、相関関係が認められない4因子に分けられる確証のところの意味もわからないのです。
素直に、解析をし直すべきでしょうか?

主因子法とバリマックス回転を用いた因子分析をされていますが、理由が見当たりません。言えなかった状況9カテゴリーと、言わなければよかった状況4カテゴリーから作成された尺度が、必ず相関関係が認められない4因子に分けられるという確証が無い限りは、今時は最尤法とプロマックス回転などの斜行回転を用いるべきであると考えます。

Re: 尺度開発
  投稿者:ポン太 2016/11/15(Tue) 15:05 No. 22212
どのような書籍を参考にされたのかわかりませんが、査読者の指摘の通りだと思います。

主因子法は別にして、実質科学からして、相関関係をゼロと仮定するパリマクスが理由もなく
使われていたら、因子分析を理解されている査読者なら、??となりますよ。

納得いかないのであれば、編集委員会に伝えたらいいでしょう。しかし今回はおそらく却下されてしまいますよ。

良書を選んで読まれてはいかがでしょうか。

Re: 尺度開発
  投稿者:kashio 2016/11/15(Tue) 15:28 No. 22213
ありがとうございます。良書をご紹介いただければと思います。
プロマックス回転に挑戦します。この査読者は、これだけですが、
もう一人の査読者は以下のように言っています。共分散は次の論文でと思っていたのですが、プロマックスに挑戦してもダメでしょうか?

α係数も、0.6台が3因子と、正直高いとは言えません。そのことも考えると、確証的因子分析で本尺度項目での因子構造の確認を行うことまでしたほうがよいのではないでしょうか。

Re: 尺度開発
  投稿者:ポン太 2016/11/15(Tue) 15:49 No. 22214
はい、確認的因子分析まで持ち込む必要があると思います。

「共分散はつぎの、、、、 ダメでしょうか」の意味が理解できません。分析には、分散共分散行列、もしくは相関行列が使われているはずなのですが。

いずれにせよ、あなたがやるべきことは反論ではなく、2人の査読者の指摘されたことだと思います。

Re: 尺度開発
  投稿者:ポン太 2016/11/15(Tue) 15:51 No. 22215
書き忘れましたが、因子間相関がゼロと言えるのであれば、反論の余地はあるでしょう。

Re: 尺度開発
  投稿者:kashio 2016/11/15(Tue) 17:27 No. 22216
ありがとうございました。確証的因子分析してあります。今回は質的研究で質問項目を起こしました。質的研究で分けたカテゴリーから項目を作っていった過程の記述を丁寧に書いてあるので紙面が足りず探索的因子分析までにしようと思いました。頑張ってみます。

お願いです。良書を紹介してください。

Re: 尺度開発
  投稿者:kashio 2016/11/15(Tue) 22:24 No. 22217
ポン太さま

やっぱり、わかったふりはやめます。今時 に引っかかっていたのは、何年か前に還暦をすぎていて、主因子法バリマックス回転で何本か論文を書いているから!でも、今時 だからではなく、因子間相関がゼロでなければ、バリマックス回転はしない。ということがわかりました。バリマックスで一番納得がいく下位尺度ができたのですがー
残念ですが、最尤法プロマックスでやってみます。
確証的因子分析はやってあるし、いわゆる当てはまりもいいようです。
でも、なぜ、確証的因子分析をしなくてはならないのかわかっていません。
反論をしたいのではなく、理解したいのです。査読者は、アルファー係数が0.6と低いので、確証的因子分析を進めていますが、アルファ係数は信頼係数だと理解しています。ご助言ください。

Re: 尺度開発
  投稿者:ポン太 2016/11/16(Wed) 08:43 No. 22218
確認的因子分析をされているので。その意味は理解されていると思ったのですが…。因子構造の単純化です。探索的および確認的因子構造を作図してみれば一目瞭然かと。作図しなくても理解はできますが。

当てはまりはいいかもしれませんが、因子間相関をゼロと仮定するか否かのところで、近年の統計学の理論が反映されていないのが問題です。ただわたしは、純粋な統計学者ではありませんので、数式などを使っての説明はできません。

アルファ係数は高いに越したことはないでしょうが、大学の成績の秀優良可くらいに考えておけばいいのでは?0.6なので可ですか。まあ、この基準を論文に書けるかは別問題になりますが。そもそもこうあるべきというルールはないと思います。アルファ係数は内的信頼性係数のことでいいと思います。青木先生のRの関数を使って、係数を高めることは可能だと思います。もちろん数値だけをみてバンバン削除するのは違いますが。

豊田秀樹氏、狩野裕氏の著書が参考になります。もちろん他にもあります。

Re: 尺度開発
  投稿者:kashio 2016/11/16(Wed) 10:38 No. 22219
ポン太様

ありがとうございます。「理解した」とは自分の言葉で語れること,と思っています。
まだ語れません。きちんと理解できるように、まず、ご紹介いただきました書籍を読みます。

Re: 尺度開発
  投稿者:青木繁伸 2016/11/16(Wed) 23:29 No. 22220
因子間に相関があることを仮定しているのがプロマックス回転であるからして,バリマックス回転はそれに含まれる(相関が0である特殊な場合)。
よって,プロマックス回転をして,ある因子間の相関はほとんど0であるといえばよいだけの話。すべての因子間相関が0に近いならば,「バリマックス回転でよいのだ」といえるわけだけど,因子間相関が0であることをことさらい言う意味がどれほどあるかということ。

まあ,昔は(コンピュータが使えなかったり(そんなばかな),斜交回転までは検討できなかったりで)主因子解・バリマックス回転を行うのが主流でもあったのだけど,現在に至っては「いまさら」なんでしょうね。査読者の「今の時点じゃ,そんなの当たり前でしょ」という,上から目線がや〜〜な感じなのは仕方ない(あなた方も,少し前までは主因子解でバリマックス回転で論文書いてたでしょ...と今更いっても仕方ない)。

確証的因子分析の必要性も,初期の頃には重視されていたとは言いがたい。

尺度が使用に耐えるためにはアルファ係数は 0.8 以上というのは教科書的ではあります。0.7 以上ならまあいいのではないかとかいう人もいたりするわけではありますが,0.6 で十分というのは,私は見聞きしたことはないですね。

統計学を用いた論文の論文たる水準も時代によって変化する(上昇する)ので,そのときどきの水準を満たすことは学術論文として必須のことではあるのでしょう。

いずれにしろ,学会誌に論文を掲載してもらうにはまずは査読者の判定を経なければならないわけで,正当な反論がない限り泣く子と地頭には勝てないので,査読者の指示に従うしかないということになるでしょう。

査読者と喧嘩をしてみるのも面白いかも知れませんね。

Re: 尺度開発
  投稿者:kashio 2016/11/18(Fri) 16:12 No. 22221
ポンタ様、青木先生ありがとうございました。
最尤法、プロマックス回転をしました。
主因子法、バリマックス回転と同じ因子、項目数でした。因子固有値が高くなりかえって良い結果でした。下位尺度ごとに相関はありました。
もう一つ査読者が 項目分析の結果を書けと言ってきました。
これが良くわかりません。分布の偏りや、標準偏差を確かめたり
しました。下位尺度ごとにクロンバックを上げるように項目を出し入れしました。
その都度テスト再テストや基準関連妥当性も見ました。

クロンバックですが、査読者は0.6しかないといっていますが、
第1因子0.68、第2因子0.75、第3因子0.68、第4因子0.68 全体0.81
でした許されませんか?隣の部屋の先生に勧められて、KMOとBartlettなるものもやって、標本妥当性なるものもみました。隣の先生は有名な先生ですけど項目分析はしらないといってます。項目分析を調べたところクロンバックがなんとかなったのでよいのではないかと思うのですが、この項目分析も今時のものですか?、

Re: 尺度開発
  投稿者:青木繁伸 2016/11/18(Fri) 21:04 No. 22222
クロンバックのアルファですが,教科書的には 0.8 というのが閾値ですが,「実用では 0.7 でもよい」とか書いてあるサイトも見かけましたが,その根拠については確認していません。

わたしも,項目分析についてはよくわかりません。

Re: 尺度開発
  投稿者:kashio 2016/11/19(Sat) 21:15 No. 22223
ありがとうございます。参考書を読んでもなぜかがわからず苦しむことがおおいです。もう少しまなんで、探索的因子分析だけでなく、なぜ確証的因子分析をしなくてはならないか理解したいです。


二つの級内相関係数の差の検定
  投稿者:NIIGATA 2016/11/09(Wed) 22:59 No. 22210
お世話になります。
唐突で大変恐縮ですが、2つの級内相関係数の差の検定なるものは存在するのでしょうか?

そして、存在するならば、それは相関係数の差の検定とはどのように異なるのでしょうか?

可能でしたらRのスクリプトも含めご教示いただけますと幸いです。


工程能力指数Cp
  投稿者:かずみ 2016/11/04(Fri) 21:32 No. 22207
「工程能力指数Cpは、規格の中心と標本平均値Xberがほぼ同じ時に求める。」となっていますが、このほぼ同じ位とはどれ位までなのですか?
Cpkは、こうしたことを考えなくて良いので、実務ではCpkを用いていますが、
Cpの場合は、この辺りが曖昧でわかりません。検定で同じとは言えないので、分散分析などのプーリングの考え方などから考えてαを25%と考えてこれを両側として、、、、すると±0.3シグマ位をほぼ同じと考えてよいのでしょうか?
このほぼ同じくらいといった事を統計的に解説した書物、論文はあるのでしょうか?

Re: 工程能力指数Cp
  投稿者:まさき君 2016/11/09(Wed) 08:19 No. 22208
品質管理上の経験則であって、理論的に説明できないとは思いますが、工程ばらつきにおいては、ロット毎に平均値は変動するが、ロット内の分散は殆ど変動しない、と考えるのが一般的で、これに基づけば、Cpにも一定の意味合いはあります。
例えば、Cpkが全然駄目な時、Cpは十分ならXberを補正する事で歩留り改善の余地があるが、Cpも駄目なら手のつけようがない、等。

Re: 工程能力指数Cp
  投稿者:かずみ 2016/11/09(Wed) 18:10 No. 22209
ありがとうございます。
やはり工程能力指数はCpkとしてXberが補正、調整可能な場合の工程能力指数がCpとかんがえるべきなのでしょうか?
であれば色々な書物、文献はそうすべきですね。


集団と個体との比較
  投稿者:さかな 2016/10/27(Thu) 16:11 No. 22194
ある個体が,集団と同じレベルといえるかどうかを調べようとしています。
既存の集団があって,その集団から1個体を抜き出してある処理をした上で,元の集団と差があるかどうかを調べるというものです。
抽象的で分かりにくいので,架空の具体例でお伝えします。

トリュフを探すための豚を80頭飼っており,その中から1頭(Aとする)を抜き出して嗅覚が良くなる薬を投与します。その後,50日間,全頭で同時に同じところでトリュフを探させる実験を毎日します。毎日異なる土地へ連れて行って,探しやすさなどの条件を変えます。日ごとに,各豚が見つけたトリュフの数を数えます。全ての豚は常に個体識別がしっかりとできる状態です。

薬の効果を確認するために,Aが有意にトリュフ探しの能力が向上したか否かを統計的に調べるにはどのようにしたら良いでしょうか。また,もし能力が向上したとしたら,どの程度向上したのか(1.4倍とか)をどのように判断すれば良いのでしょうか。

調べようとする試験区個体が1個体だけでは無理があるのかもしれませんが,この場合非常に高価な薬であるために,1頭分しか用意することができませんでした。また,Aは50回全ての日に参加していますが,その他の豚は参加しなかった日があるため,対照区の豚の頭数は日によって異なります。

Re: 集団と個体との比較
  投稿者:青木繁伸 2016/10/27(Thu) 18:17 No. 22195
実験計画に失敗しましたね
統計学的に意味のないデータです
高い薬と時間と人手が無駄になりました

Re: 集団と個体との比較
  投稿者:鈴木康弘 2016/10/28(Fri) 07:41 No. 22196
その80頭を、1番鼻のよい豚、2番目に鼻のよい豚、と並べることができるなら、薬を与えた豚が2番目までに入っていれば、5%有意水準の両側検定で有効だと、言えるかも?

Re: 集団と個体との比較
  投稿者:青木繁伸 2016/10/28(Fri) 08:45 No. 22197
豚が1頭だと,その測定値が何個得られようと,測定値が表すのは個体内変動であり,その測定値の平均値が大きくても,その豚に限ってはそうかも知れないけど,他の豚にも同じ効果があるとはいえないのではないか?に対しては反論できない。測定回数を増やしても,その豚の特性の測定精度が高くなるだけあり,効果が一般的であるかどうかについては回答できないままである。
統計学は,複数の個体について測定し,その測定値が表すのは個体間変動であり,測定値の平均が誤差変動よりも大きければ効果があるとする。効果が顕著な豚もそうでない豚もあるけど平均すれば効果がある。この場合には,サンプルサイズが大きくなれば効果の測定値の精度が高くなる(より確かな推定値が得られる)。

「薬を与えた豚が2番目までに入ってい」ても,「その豚はそうだったでしょうけど」ということに変わりはない。統計解析手法の選択の問題ではなく実験計画時点での問題。データを得たあとで分析不能があきらかになるという不都合。

Re: 集団と個体との比較
  投稿者:さかな 2016/10/31(Mon) 10:12 No. 22200
青木先生,鈴木先生

コメントいただきまして,誠にありがとうございます。
実験する豚が1頭だけではその個体に対する効果だけしかわからないので回数を重ねても同じで意味がないということ,理解いたしました。ありがとうございます。

では,この豚に対するこの薬の効果を知りたい場合には,何らかの方法で効果の有無を判断することができるでしょうか。以下のようなことを考えてみましたが,いかがなものでしょうか。

ここで,すべての豚の50回(50日)の成績(トリュフの個数)がわかっています。
例えば,50回のそれぞれにおける全体の中でのAの位置を知るために,偏差値を計算し,全ての回で偏差値66.4以上(上位5%以内)であったとしたら,有意にAの成績が良かったと考えてよろしいでしょうか。
あるいは全ての回で偏差値33.6(下位5%)〜66.4(上位5%)であった場合に,有意に成績が良かったり悪かったりすることがないといえるでしょうか。
さらに,50回のうち何回か(例えば9回)偏差値70以上のケースがあった場合には,どのように判断すればよいでしょうか。

最初と少し議論がずれてしまいますが,薬の豚に対する効果ではなく,あくまでこの個体に対する効果を見たい場合です。

Re: 集団と個体との比較
  投稿者:青木繁伸 2016/11/01(Tue) 08:24 No. 22201
その個体の一連の測定値の中の最小値が他の普通の個体のデータと一緒に外れ値の検定を行って,外れ値だといえればよいかな。

最小値を使うのは,最も控えめに言うため。平均値でもよいかも知れないけど,批判する人は批判するだろうから。

Re: 集団と個体との比較
  投稿者:さかな 2016/11/01(Tue) 09:09 No. 22202
青木先生

ありがとうございます。
少し説明不足のところがあったようです。すみません。

この実験は,毎日条件の異なる土地で実施ています。したがって,それぞれの豚が見つけたトリュフの数は1桁の時もあれば,3桁の時もあります。全体的にほかの豚が多く見つけるときはAも多く見つける傾向がみられます。したがって,各実験日のデータは独立(という用語が適切かわかりませんが)であるため,Aの測定値の最小値や平均を使うのは適切ではないのではと考えます。
つまり,各実験日のデータごとに比較をして,50日分のデータを合わせて差があるといえるのかどうかを判断したほうが良いのではないかと思います。
いかがでしょうか。

その場合には,実験日ごとに外れ値の検定をして,もちろん全実験日で外れ値またはそうでないとの結果であれば,それぞれ差がある,または無いと判断できますが,どちらの場合も含まれるときに,どのようにして判断すればよいのか,ご教授いただけたら幸甚に存じます。

Re: 集団と個体との比較
  投稿者:青木繁伸 2016/11/01(Tue) 09:43 No. 22203
> この実験は,毎日条件の異なる土地で実施ています。したがって,それぞれの豚が見つけたトリュフの数は1桁の時もあれば,3桁の時もあります。全体的にほかの豚が多く見つけるときはAも多く見つける傾向がみられます。したがって,各実験日のデータは独立(という用語が適切かわかりませんが)であるため

その一頭の豚のデータは独立ではないですね。多い日も少ない日もあるでしょうが,潜在能力は一定なのでしょうから。

しかし,まあ,ここまで「測定値」が特徴的だと,どうしようもないのではないですか。

一頭ではあるが,その一頭(他の豚も)なまじ複数回容易に測定できるから複雑になる。
複数頭からの1回のデータというのが大原則。

Re: 集団と個体との比較
  投稿者:さかな 2016/11/01(Tue) 11:46 No. 22204
すみません,例を変えてお聞きいたします。

湖にアヒルが80羽いて,そのうち1羽(Bとします)が野犬に襲われてけがをしてしまいました。嘴にけがをして,餌を食べにくい状況です。

このBを保護していますが,Bが他のアヒルと同じくらい餌を食べるのであれば,湖に返そうと思います。

そこで,1か月間毎日,Bの餌の摂取量とほかのアヒルの摂取量をそれぞれ確認しました。ただし,ここでは1カ月間での回復はないものと仮定します。

このデータによって,Bの餌の摂取量が他のアヒルと比べて有意に少ないか否かを確認することはやはり無理なのでしょうか。

Re: 集団と個体との比較
  投稿者:鈴木康弘 2016/11/03(Thu) 07:15 No. 22205
青木先生の言われるように、外れ値の検定では、いかがでしょう?

Re: 集団と個体との比較
  投稿者:さかな 2016/11/04(Fri) 13:28 No. 22206
青木先生,鈴木先生

いろいろとご教授いただきまして,誠にありがとうございました。


データ間の相関
  投稿者:初心者 2016/10/29(Sat) 13:33 No. 22198
青木先生、いつもお世話になっております。

Aという病気をもつ患者50名とBという病気をもつ患者20名に、それぞれ5つの心理検査を行いました。5つの心理検査のデータのうち、2つの心理検査間での総得点の相関や各心理検査の下位項目間での得点の相関を調べるときに用いる統計学的手法は因子分析になるのでしょうか?ちなみに結果についての仮説はありません。

初歩的な質問で大変申し訳ありませんが、どのような統計学的手法を用いたらいいのか分からず、困っています。

恐れ入りますが、どうぞご教示の程、よろしくお願い致します。

Re: データ間の相関
  投稿者:青木繁伸 2016/10/29(Sat) 18:26 No. 22199
いくつかの変数間の相関関係を見るのなら相関係数と言うことになるのでしょうけど
因子分析をするにしては,サンプルサイズが小さいし
A,B の二疾患をもつ患者をひとまとめにして分析してもよいかということもありますね
初心者とはいえ,いろいろなことを知っていないと適切な分析はできませんねえ
(そもそも,心理検査って,疾患ごとに特徴的な所見が得られたりするんでしょうかね)


回帰直線の傾きの意味
  投稿者:さかな 2016/10/26(Wed) 16:40 No. 22189
初歩的な質問で恐縮ですが,よろしくお願いいたします。

A, B の 2 群の data を比較し,B が A と同じなのか違うのか,違うのであればどの程度違うのかを知りたいのですが,散布図から回帰直線を求めました。

このときの直線が B=0.95A だったとします。
この場合に,B は A に対して 0.95,即ち,「B の値は A の値の 95%である」といえるのでしょうか。

Re: 回帰直線の傾きの意味
  投稿者:青木繁伸 2016/10/26(Wed) 17:22 No. 22190
A 群,B 群のデータを散布図に描いて回帰直線を求めるって?具体的にどうやったのですか?

同じ対象に対して A 法と B 法により二つの測定値を得て...,ということですか?

Re: 回帰直線の傾きの意味
  投稿者:さかな 2016/10/26(Wed) 17:54 No. 22191
早速コメントいただき,誠にありがとうございます。

説明が不十分でした。誠に失礼いたしました。

A群とB群とで対応のあるデータです。図のようになります。


Re: 回帰直線の傾きの意味
  投稿者:青木繁伸 2016/10/26(Wed) 19:27 No. 22192
「A群とB群」という言い方は避けるべきですね。

さて,Y = βX + ε という場合ですが,Y と X の比を推定するには,3 通りの方法があります(スネデカー,コクラン 統計的方法 原書第6版,岩波書店 P.161-165)

(1) εの分散が X に関わらず一定であると仮定する場合には b = ΣXY / ΣX^2
(2) εの分散が X に比例すると仮定する場合には b = ΣY / ΣX = Yの平均 / X の平均
(3) εの分散が X の 2 乗に比例すると仮定する場合には b = Σ(Y / X) / n
(1) は原点を通る回帰直線,(2) は平均値の比,(3) は個々の比の平均値であるが,これらはすべてβの不偏推定量である。どれが最善かは,εの分散の挙動次第である。

同書の最後に,以下のように書かれている。
「これら3つの推定値のどれかを用いようとするときは,データを作図して Y が X に比例するかどうかを常に検べ(ママ),必要であれば,直線が原点を通るという帰無仮説を検定すべきである。どれかの形の比推定値を軽率に採用すると,Y/X が X につれて一定ではないという情報をうしなうことになるかもしれない。」

Re: 回帰直線の傾きの意味
  投稿者:さかな 2016/10/27(Thu) 09:19 No. 22193
青木先生

用語の使い方が誤っておりましたこと,お詫び申し上げます。

丁寧にご説明いただき感謝いたします。

「スネデカー,コクラン 統計的方法」も確認させていただきました。 

誠にありがとうございました。


ロジスティック回帰、Tukey法
  投稿者:ぐれいぷ 2016/10/19(Wed) 09:35 No. 22179
初めて質問します。よろしくお願いします。
生物系の実験をしていて、以下のようなデータの解析をしたいと考えています。
> a <- read.table("clipboard",header=T)
> a
med tree a N
1 A P 2 16
2 A Q 1 15
3 A R 1 16
4 A S 0 17
5 B P 10 16
6 B Q 9 15
7 B R 5 16
8 B S 8 17
9 C P 5 16
10 C Q 4 15
11 C R 4 16
12 C S 0 17
(4本の樹(P,Q,R,S)の葉に3種の薬剤(A,B,C)を処理して、処理した葉N枚のうち、a枚で障害が出た、というデータです。)

 薬剤の効果を比較したいと思い、ロジスティック回帰後にRのmultcompのパッケージを用い、Tukey法で検定しました。
> med <- a$med
> r <- glm(cbind(a$a,a$N-a$a)~med,family=binomial)
> mul <- glht(r,linfct=mcp(med="Tukey"))
> summary(mul)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glm(formula = cbind(a$a, a$N - a$a) ~ med, family = binomial)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B - A == 0 2.7081 0.5737 4.720 < 0.001 ***
C - A == 0 1.3412 0.6026 2.225 0.06470 .
C - B == 0 -1.3669 0.3988 -3.428 0.00169 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

 薬剤の効果は比較できるのですが、「1つの薬剤に関しては単純な4反復になっていて、上の解析では同じ樹を使っている効果が除外されているのでは?」と質問を受けました。

 薬剤の効果と樹の効果については、パッケージcarを使って、
> r2 <- glm(cbind(a$a,a$N-a$a)~a$med*a$tree,family=binomial)
> Anova(r2)
Analysis of Deviance Table (Type II tests)

Response: cbind(a$a, a$N - a$a)
LR Chisq Df Pr(>Chisq)
a$med 36.264 2 1.335e-08 ***
a$tree 7.878 3 0.04861 *
a$med:a$tree 8.477 6 0.20523
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
と解析してみたのですが、上記のTukeyで、樹の効果が分かるような解析方法はありますでしょうか?

Re: ロジスティック回帰、Tukey法
  投稿者:太郎 2016/10/19(Wed) 14:38 No. 22180
 単純に a/N(障害発生葉の比率)をデータとして、薬剤3水準と樹4水準の乱塊法(繰り返しのない2元配置法)を適用するというのはいかがですか?
 
 比率型データなので角変換してから分散分析して、薬剤、樹の多重比較(Tukeyなど)をしてもよいと思います。

Re: ロジスティック回帰、Tukey法
  投稿者:ぐれいぷ 2016/10/19(Wed) 23:43 No. 22181
太郎様
 回答していただきありがとうございます。私も少し前まで、比率の検定→変換という方法を用いておりましたが、変換して分散分析するよりも一般化線形モデルを用いたほうがいいと聞きましたので、上記の方法を用いました。
 上記で樹の効果が分かる方法がなければ、角変換してみます。

Re: ロジスティック回帰、Tukey法
  投稿者:鈴木康弘 2016/10/20(Thu) 16:43 No. 22183
 すみません、multcomp というパッケージのことは初めて聞いたし、ダウンロードもしてないのですが、

> r <- glm(cbind(a$a,a$N-a$a)~a$med+a$tree,family=binomial)

 とやって glht コマンドを使っても、ご希望のようにいかないのでしょうか?

Re: ロジスティック回帰、Tukey法
  投稿者:ぐれいぷ 2016/10/21(Fri) 10:28 No. 22184
鈴木康弘様
 回答ありがとうございます。はい、多重比較の結果は出るのですが、その中にtreeの効果は出ません。

Re: ロジスティック回帰、Tukey法
  投稿者: 2016/10/21(Fri) 20:15 No. 22185
同じ樹から繰り返しデータを取っているため一般化線形混合モデルの適用ではないでしょうか。
ただしこの種の解析をしたことがないため、望みの結果が得られているかどうかは自信がありません。
> library(lme4)
> r <- glmer(cbind(a, N - a) ~ med + (1|tree), data = a, family = binomial)
> summary(r)

Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: cbind(a, N - a) ~ med + (1 | tree)
Data: a

AIC BIC logLik deviance df.resid
51.4 53.3 -21.7 43.4 8

Scaled residuals:
Min 1Q Median 3Q Max
-1.7409 -0.2627 0.3620 0.4520 0.7766

Random effects:
Groups Name Variance Std.Dev.
tree (Intercept) 0.1319 0.3632
Number of obs: 12, groups: tree, 4

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.7576 0.5521 -4.994 5.90e-07 ***
medB 2.7654 0.5808 4.761 1.92e-06 ***
medC 1.3595 0.6053 2.246 0.0247 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
(Intr) medB
medB -0.849
medC -0.805 0.766

> library(multcomp)
> mul <- glht(r, linfct = mcp(med = "Tukey"))
> summary(mul)

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts

Fit: glmer(formula = cbind(a, N - a) ~ med + (1 | tree), data = a,
family = binomial)

Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
B - A == 0 2.7654 0.5808 4.761 < 1e-04 ***
C - A == 0 1.3595 0.6053 2.246 0.06171 .
C - B == 0 -1.4059 0.4065 -3.458 0.00151 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)


Re: ロジスティック回帰、Tukey法
  投稿者:ぐれいぷ 2016/10/21(Fri) 22:22 No. 22186
荒様
 回答ありがとうございます。Rを使って解析していただいて助かります。一般化線形混合モデルは今まで扱ったことがないので知りませんでした。勉強してみます。
 今後ともよろしくお願いいたします。

Re: ロジスティック回帰、Tukey法
  投稿者: 2016/10/22(Sat) 20:02 No. 22187
後で気が付いたのですが、このような対応のあるデータにTukey法を適用できないのではと思いました。

ですので、この場合であれば
A-B, A-C, B-C の3通りの組み合せで解析を行い、その後P値をHolm法で調整するのがいいのではないでしょうか。

combn関数で各組み合わせが配列として得られるため、
以下のプログラムで検定結果が得られるはずです。
このままでは結果の表示が貧弱なため、もう少し改変が必要ですが。

comb <- combn(unique(a$tree), 2)
n <- length(comb)
p <- numeric(n) # P値を入れる変数

for (i in seq(n)) {
tmp <- subset(a, tree %in% comb[, i])
r <- glmer(cbind(a, N - a) ~ med + (1|tree), data = tmp, family = binomial)
p[i] <- summary(r)$coefficients[2, 4] # P値
}

p.adjust(p) # デフォルトではHolm法で補正

Re: ロジスティック回帰、Tukey法
  投稿者:ぐれいぷ 2016/10/24(Mon) 14:40 No. 22188
荒様
 たびたびお答え下さりありがとうございます。compなど、今まで使ったことが
ないため、勉強して使ってみます。


2つの関連を求める
  投稿者:たつ 2016/10/18(Tue) 05:50 No. 22178
統計素人です。
Aという物質の乾燥させた後(水分を飛ばした後)の質量を測定しています。本来は、70度で30分間乾燥させた時の質量なのですが、これでは時間がかかるということで、105度で5分間の乾燥をさせた時の質量を測定に変更することを考えています。その際に、統計的な手法を用いて、両方の相関はより70度で30分間乾燥から105度で5分間に変更しても良いと結論を導き変更したいと考えています。
70度で30分間乾燥した後の質量と105度で5分間の乾燥した後の質量の相関は、どのような方法で求めればよいのでしょうか?またその例数はどれくらいが適切でしょうか?

Re: 2つの関連を求める
  投稿者:鈴木康弘 2016/10/20(Thu) 07:01 No. 22182
 乾燥させたら、1つのサンプルを2つの方法で測定できない?というのは置いておき、
Rなら、相関係数は cor() や cor.test() で求めます。

 必要なサンプルサイズは下記が参考になりますか?

http://aoki2.si.gunma-u.ac.jp/R/power-cor-test.html


分位数で区切られた集合の名称は?
  投稿者:渡辺比登志 2016/10/07(Fri) 17:06 No. 22166
統計学の素人ですが、英語の医学論文を読んでいて常々気になっていることがあります。
例えば、「四分位数」(または四分位点)は英語で「quartile」です。「第1四分位数」は「first quartile」(quartile 1: Q1)です。ところが、英語の論文ではしばしば、「最低値から第1四分位数までの範囲に含まれるデータの集合」の意味で「quartile 1」「Q1」が使われています。つまり、全データを4等分した各部分集合をQ1、Q2、Q3、Q4としている論文が多いのです。
これは正しい用語法なのでしょうか。いろいろ調べても「quartile」の語にそのような意味は見つからず、また上記のような部分集合を表す英語表現も見つかりません。この習慣的用語法が誤りだとすると、4等分された各部分集合を英語で正確にはなんといえばよいのでしょうか。
一方、日本語ではこの部分集合を表すのに「四分位」の語が使われることがあります。しかし、「四分位」に相当する英語はやはり見つかりません。日本語の「四分位」は正しい用語なのでしょうか。

Re: 分位数で区切られた集合の名称は?
  投稿者:鈴木康弘 2016/10/10(Mon) 08:07 No. 22170
 >全データを4等分した各部分集合

 これはlower 25%, second 25%, third 25%, upper 25% じゃないでしょうか?
 日本語なら第一四分位範囲、第二四分位範囲...では。

Re: 分位数で区切られた集合の名称は?
  投稿者:渡辺比登志 2016/10/11(Tue) 17:07 No. 22171
鈴木康弘さま、ご回答ありがとうございました。
英語ならたしかに「lower(first) 25%, second 25%, third 25%, upper(fourth) 25%」といえますね。これ以外の(もっと形式的な、専門用語らしい?)英語表現はないのでしょうか。私は医学論文しか知りませんが、他の分野では「lower 25%, second 25%, third 25%, upper 25%」を使っているのでしょうか。その場合、(quartileに対する「Q1」のような)略号はあるのでしょうか。
私が「英語表現が見つからない」と書いたのは、これら4つの部分集合の“総称”のことですが、それはあるのでしょうか。なぜ医学論文では、「quartile」が(誤って?)使われるのでしょうか。
「四分位範囲」(inter-quartile range)は“第3四分位数と第1四分位数の差”で、つまりsecond 25%とthird 25%が含まれる範囲の大きさのことですから、「lower 25%, second 25%, third 25%, upper 25%」の意味ではありませんね。「第1四分位範囲」という表現はできないと思います。
引き続きご教示をお願いします。

Re: 分位数で区切られた集合の名称は?
  投稿者:鈴木康弘 2016/10/13(Thu) 06:56 No. 22172
>これら4つの部分集合の“総称”のことですが、それは

 inter-をつけない "quartile ranges" ではいかが?

Re: 分位数で区切られた集合の名称は?
  投稿者:渡辺比登志 2016/10/14(Fri) 12:47 No. 22173
"quartile ranges" という表現が統計学用語にあるのでしょうか。ネットでは見つけることができません。
"range" は「範囲」を意味しますから、その範囲に含まれるデータの集合を表すことはできないと思いますが。

Re: 分位数で区切られた集合の名称は?
  投稿者:鈴木康弘 2016/10/16(Sun) 06:48 No. 22174
>ネットでは見つけることができません。

 そうなんですか?

>その範囲に含まれるデータの集合

 という質問でしたら、わかりません。"data in 〜" とか?

Re: 分位数で区切られた集合の名称は?
  投稿者:青木繁伸 2016/10/16(Sun) 21:45 No. 22176
たぶん,一言で表現できる用語はないのでしょうね。
日本語の学術用語は,漢字を使って2〜6文字で端的に表現しているものの,ホントにそれでよいのか?というのもあるでしょう。なんたって,それまでになかった概念を表す用語を作るわけだから,statistics を統計学,philosophy を哲学といわれても,「哲」って何?というようなものでしょう。
英語などでは,日本語(学術用語)で数文字で表現されるものも,文章(の一部)として表すことも多いですよね。sum of squares (残差平方和)とか Spearman's rank corerration coefficient(スピアマンの順位相関係数)とか。よくあるのは root mean squares -- RMS と略称されますがこれ,「二乗平均平方根」ですけど sqrt(sum((x-mean)^2)) なら,標準偏差(とか標準誤差)です 。しかし,分野によってはそれ以外のものを指します(計算式は同じだし,数学的な意味的には同じなのだけど)。
よって,「もしそのような用語がない,あるいは定義が不正確・不十分ならば,言語にかかわらず,文章で説明するしかない」というのが解答のように思いますが。

つまり,今の場合ならば,「第一1四分位数より小さいデータの集合」,「第一四分位数より大きく,第二四分位数より小さいデータの集合」...とか,「データを小さい順に並べて同じ数ずつに分割した場合の一番小さいグループ,二番目に小さいグループ...」ということで,英語での記述もそれを英語で述べると...

科学的な表現というのは正確性が必要なのだから,ある概念を数語(数文字)の一般的に了解されていない用語を使って表すことは不適切...ということではないでしょうか?
関係ないけど,「"avoidant purchasing" なんて変だろう "purchase avoidance" とか "consumer avoidance" じゃないか」とか指摘している人がいるけど,「これこれのことを以下ではこう言うよ」というように宣言してから使用すれば,問題ないんだよね(明らかにおかしくなければ,なんだけど)。

もし,端的な表現・用語があるならば,査読者が教えてくれると思います。
不適切,あるいはわかりにくい表現だなあと思う場合は,「ネイティブにはかなわんなぁ」とかあきらめて,次回以降にはそれを使うとか?(^_^)

ちなみに,医学論文に限らないけど,論文に使われている(原著論文として採用されているかに限らず)表現・用語が常に正しいという保証は,ないですね...(著者が統計学に詳しくなくて,査読者,編集委員会がヘボだとそのまま通ってしまう。日本の学術論文では特に)。

Re: 分位数で区切られた集合の名称は?
  投稿者:渡辺比登志 2016/10/17(Mon) 17:35 No. 22177
青木繁伸先生、丁寧なご回答ありがとうございました。

欧米人が英語で誤用をして、それらを読んだ日本人学者がそれを正しい英語表現だと思い込んで、また日本語でさらに誤って「四分位範囲」と書いたり「四分位」と適当な語を作ったりしているのではないかと想像しています。英語の「quartile」を「四分位数」「四分位点」と正確に訳してみれば(ある範囲に含まれるデータの意味にはならないので)、その英語が誤用だと気づくはずですが。

「lower 25%」という鈴木康弘様のご指摘にヒントを得て、「lowest quarter」で検索してみたら、下記などいくつか見つかりました。

Weerts MJ, et al. Oncotarget. 2016 May 17;7(20):29166-76.

「lowest one-third」が出てくる論文も見つかりました。おそらくこれらが正しい用法かと思います。なんということはないですね。日本語でどう言えばよいのかわかりませんが。
どうもありがとうございました。


GLMM(ロジスティック回帰)について
  投稿者:植田 2016/10/07(Fri) 13:17 No. 22165
いつもお世話になっております。
R3.3.1のlme4パッケージ内にあるglmer関数を用いて以下のデータの分析を行いました。

ID1 ID2 response fix
12 9 1 0.5
19 12 0 0.5
25 21 0 0.5
22 7 0 0.5
22 41 0 0.5
12 40 0 0.5
12 40 0 0.5
12 36 0 0.5
10 47 0 0.5
3 32 0 0.5
3 5 0 0.5
3 9 0 0.5
40 17 0 0.5
21 50 0 0.5
25 53 1 -0.5
22 30 1 -0.5
22 30 1 -0.5
12 11 1 -0.5
49 57 1 -0.5
42 42 1 -0.5
42 42 1 -0.5
42 54 1 -0.5
42 41 1 -0.5
56 49 1 -0.5

model<-glmer(response~fix+(1|ID1)+(1|ID2), family="binomial")

という式で分析を行うと、「 警告メッセージ:
checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, で:
Model is nearly unidentifiable: large eigenvalue ratio
- Rescale variables?」という警告メッセージが表示されました。
そして、summaryで出力すると、明らかにおかしな結果がでてきます。
ちなみに、少しデータをいじると(responseのデータに0を増やしました)、それらしき出力結果がでてきました。
このデータセットの構造のどのような点が問題であり、どのようにすれば改善することができるのでしょうか。

まことに恐縮ですが、ご教授いただけましたら幸いに存じます。
何卒よろしくお願い申し上げます。

Re: GLMM(ロジスティック回帰)について
  投稿者:青木繁伸 2016/10/07(Fri) 18:35 No. 22167
この関数はよくわかりませんが,
response~fix+(1|ID1)+(1|ID2)
って,response を fixt と ID1, ID2 を使って説明するということですよねきっと。
でも,response と fit は
-0.5 0.5
0 0 13
1 10 1
となっているので,fix だけで「完全に説明できる」とうことで,これはこれでうれしいことなんでしょうけど,glmer を使ってまで分析するのか???という状況になっているのでしょう。

Re: GLMM(ロジスティック回帰)について
  投稿者:植田 2016/10/07(Fri) 21:51 No. 22168
青木先生

質問者の植田です。
早速のお返事ありがとうございます。
完全に説明できてしまっているということで、このような状態になっているのですね。
最初はlmerで分析をしていたのですが、警告でglmerを使えと言われました。
少し調べてみると、海外の掲示板?に、正規分布以外の場合にはlmerではなく、glmer関数を使った方がよいと書かれていました。

もしよろしければ重ねてご質問をお許しいただけたらと思うのですが、
「fixによる違いが明瞭なので、GLMMがかけられず、他の分析を行った」といった記述を論文中(英語)で行ってもよいものなのでしょうか。
類似した表現を見たことがないため、ご意見をお聞かせいただけましたら幸いです。
どうぞよろしくお願いいたします。

Re: GLMM(ロジスティック回帰)について
  投稿者:青木繁伸 2016/10/07(Fri) 22:16 No. 22169
そもそも,複雑な分析をするまでもなく,単純集計で明らかな事実を,複雑なモデルを仮定して分析するのはどうかなということでは?
結局の所はまだまだ,サンプルサイズが小さいのが問題でしょう。
データによってはそれほどたくさんのデータをとれないといこともあるかも知れませんが,だとしても,そのような少数データに多変量解析(この場合はそれほどの多変量でもないが)を使うというのが時期尚早なのでしょう。(平均値の差の検定なんていう基礎的な統計解析であっても,検定は出来るけどその妥当性はどうかな?ということもあるわけですよね)。


R 論理ベクトルの論理和の計算
  投稿者:明石 2016/10/04(Tue) 09:47 No. 22162
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

昨日は、「"を含む文字列の連結」について、色々な方法をご教示いただき、
誠にありがとうございました。

とても良い勉強になりました。
御礼を申し上げます。

Rプログラムについて、ご教示いただきたいことが出てきました。

私も自分で動くプログラムは作成できましたが、
より良い方法を知りたくて、投稿させていただきました。

どうぞよろしくお願いします。

-----------------------------------------------------

以下のような、論理値からなる行列(m行 n列)があります。

論理ベクトルは列ベクトル(縦ベクトル)で、m、nの意味は以下のとおりです。
m; 論理ベクトルの要素の個数
n; 論理ベクトルの本数

お示しをしました例では、m=30,n=10 です。
この行列の名前を mx とします。

> mx
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[4,] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[5,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[6,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[7,] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[8,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[9,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[10,] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
[11,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[12,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13,] FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE
[14,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[15,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[16,] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
[17,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[18,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[19,] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
[20,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[21,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[22,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
[23,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[24,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
[26,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[27,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[28,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
[29,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[30,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

n本の論理ベクトルの論理和を求めます。

お示しをしました例では、m=30,n=10と小さい値ですが、
実際には、m、nの値はとても大きい値ですので、
n本の論理ベクトルの論理和を効率的に計算する方法を知りたいと思います。

私は、以下のように考えました。

・まず、論理値を、0/1の2値に変換して、
 (過去の投稿で、青木先生からご教示いただきました有用な方法です)
・次に、行列の横計を求めて、
・横計が1以上ならば、1に置き換えて、
・最後に、2値化された横計を、論理値に戻す

mx <- mx + 0

rs <- rowSums(mx)

test<- (rs > 1)
rs[test] <- 1
rs <- as.logical(rs)

rs が、n本の論理ベクトルの論理和で、求めるものです。

より良い方法があると思いますので、ご教示をいただければ大変に助かります。
どうぞ、よろしくお願いします。

いつもありがとうございます。

Re: R 論理ベクトルの論理和の計算
  投稿者:青木繁伸 2016/10/04(Tue) 15:07 No. 22163
示されたプログラムで十分です。

ただいくつか無駄がありますので,無駄を省くと以下のようになります

rs <- rowSums(mx)
rs <- rs >= 1

すなわち,logical についての数値演算の際の変換は自動的に行われるので mx <- mx + 0 は不要です。
また, rs >= 1 とすれば結果は logical なので,rs <- rs >= 1 で最終的な結果が得られます。

一行でやるなら rowSums(mx) >= 1 で,記述も簡潔になります。

【御礼】 Re: R 論理ベクトルの論理和の計算
  投稿者:明石 2016/10/04(Tue) 17:14 No. 22164
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

今回も有難いご教示をいただき、大変に勉強になりました。
心から御礼を申し上げます。

ありがとうございました。


R "を含む文字列の連結
  投稿者:明石 2016/10/03(Mon) 18:14 No. 22159
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。

初歩的な質問で大変に恐縮ですが、調べものをしましたが、見つかりませんでした。

何卒どうぞ、よろしくお願いいたします。

−−−

"アヒル" の文字列から、
"を含む、以下の文字列の連結を作成したいと思います。

"動物名 == "アヒル""

ご教示をいただければ大変に助かります。
どうぞ、よろしくお願いします。

Re: R "を含む文字列の連結
  投稿者:青木繁伸 2016/10/03(Mon) 18:44 No. 22160
詳しくは ? "\"" で出てくるオンラインヘルプを見てください
s = "アヒル"
t = sprintf('動物名 == "%s"', s)
print(t)
cat(t)
t2 = sprintf("動物名 == \"%s\"", s)
print(t2)
cat(t2)
u = paste0('動物名 == "', s, '"')
print(u)
cat(u)
u2 = paste0("動物名 == \"", s, "\"")
print(u2)
cat(u2)


【御礼】Re: R "を含む文字列の連結
  投稿者:明石 2016/10/03(Mon) 21:34 No. 22161
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

今回も助けていただきました。
心から御礼を申し上げます。


重回帰での切片有り無しモデルについて
  投稿者:さくらとお城 2016/09/29(Thu) 17:24 No. 22153
いつもお世話になっております。

下記のデータで切片有無の重回帰を行った結果の解釈に悩んでいます。
R G B N
1 99.0 128.0 27.0 1.61
2 91.5 99.5 40.5 2.16
3 90.0 113.0 30.0 2.32
4 72.5 72.0 45.0 3.07
5 80.0 77.5 46.0 2.97
6 110.0 142.0 24.0 1.58
7 96.5 114.5 31.5 1.87
8 98.0 119.5 30.0 2.02
9 112.5 139.5 21.5 1.61
10 88.0 104.0 36.5 2.02
11 81.5 74.0 44.0 2.87
12 90.0 102.5 34.5 1.99
13 84.5 92.5 39.0 2.17
14 99.5 120.5 34.5 2.49
15 81.0 79.0 44.0 2.99
16 91.0 103.0 38.0 2.25
17 95.5 111.5 32.5 2.07
18 85.5 88.0 42.5 2.31
19 83.0 79.5 45.5 2.59
20 81.0 89.0 38.5 2.36

data=read.table("data.txt",header=TRUE)

lm1=lm(N~R+G+B+0,data=data)
lm2=lm(N~R+G+B,data=data)

summary(lm1)
AIC(lm1)

summary(lm2)
AIC(lm2)

以下はRでの結果の一部だけです。

> summary(lm1)
lm(formula = N ~ R + G + B + 0, data = data)

Residual standard error: 0.2448 on 17 degrees of freedom
Multiple R-squared: 0.9904, Adjusted R-squared: 0.9888
F-statistic: 587.5 on 3 and 17 DF, p-value: < 2.2e-16

> AIC(lm1)
[1] 5.21581

> summary(lm2)
lm(formula = N ~ R + G + B, data = data)

Residual standard error: 0.2407 on 16 degrees of freedom
Multiple R-squared: 0.765, Adjusted R-squared: 0.7209
F-statistic: 17.36 on 3 and 16 DF, p-value: 2.764e-05

> AIC(lm2)
[1] 5.327225

切片無しモデルはAdjusted R-squared: 0.9888で,良すぎるようにも見えて,大丈夫でしょうか?

どうぞよろしくお願いします。

Re: 重回帰での切片有り無しモデルについて
  投稿者:青木繁伸 2016/09/30(Fri) 07:06 No. 22154
基本的には,「理論的に考えて」独立変数が 0 のときに従属変数は 0 になるだろうかということに尽きるでしょう。
データだけから見る場合,R, G を独立変数とした単回帰では原点を通らないので,重回帰でも R, G, B が 0 のとき原点を通るというモデルはどうなんでしょうか?
切片なしのモデルと切片ありのモデルでは決定係数の計算方法が異なるので,両者を比較することはできません。

Re: 重回帰での切片有り無しモデルについて
  投稿者:さくらとお城 2016/09/30(Fri) 08:29 No. 22155
青木先生,早速のご返事ありがとうございます。

>基本的には,「理論的に考えて」独立変数が 0 のときに従属変数は 0 になるだろうかということに尽きるでしょう。

今回のデータは,RGBでNを予測する目的で実験で得たものです。理論的にRGBが0なら,Nも0になることについて,確証もないし,確かめることも難しい。

>切片なしのモデルと切片ありのモデルでは決定係数の計算方法が異なるので,両者を比較することはできません。

これに関連した図書や参考論文があれば,ぜひ教えてください。

Re: 重回帰での切片有り無しモデルについて
  投稿者:青木繁伸 2016/09/30(Fri) 17:00 No. 22156
以前はこの件のページがあったのですが,Excel の馬鹿馬鹿しさ,および Excel をひいき・寵愛・擁護するひとからのクレームによりそのページ(と Excel に関する全ページ)は削除しました。

当たり障りのないところを再現すると

======= 引用開始
スネデッカー・コクラン著(畑村,奥野,津村訳)「統計的方法 原書第6版」岩波書店の6章18節「原点をとおる直線のあてはめ」というところに詳しく書いてあります。

SPSS はどんな風にしているかなと思ってやってみて,以下のようなことであるらしいと思います。
  x  y  x^2  x・y  予測値     y^2  予測値^2
1 1 1 1 0.40779221 1 0.16629448
2 2 4 4 0.81558442 4 0.66517794
3 2 9 6 1.22337662 4 1.49665036
4 4 16 16 1.63116883 16 2.66071176
5 3 25 15 2.03896104 9 4.15736212
6 3 36 18 2.44675325 9 5.98660145
7 2 49 14 2.85454545 4 8.14842975
8 2 64 16 3.26233766 4 10.6428470
9 3 81 27 3.67012987 9 13.4698533
10 4 100 40 4.07792208 16 16.6294485
合計 385 157 76 64.0233766

a = 0.4077922 = 157 / 385
R^2 = 0.8424129 = 64.0233766 / 76 =予測値の平方和 / yの平方和
R = 0.9178305

分散分析表
  平方和  自由度 平均平方    F値
回帰 64.023377   1   64.0233766 48.1112557
残差 11.976623   9  1.33073593
全体 76     10
回帰の平方和は「予測値の平方和」,全体の平方和は「yの平方和」

プログラムの出力に下手でさっぱりわけの分からない日本語訳の注が付いています(SPSS Release 6.1.2)。
曰く,

Note # 10572
原点(切片なしのモデル)の場合,R二乗は,回帰によって証明(ママ;「説明」のミスタイプだと思う)される原点について,yにおける変動性の割合を測定します。これは,切片を含むモデルのR二乗値とは,比較できません。

========== 引用終了

まあ,定義が違うのだから比較できるわけはないのですけどね。

R の summary.lm がやっていることを抜粋すると以下のようになります。
> func = function(z) { # R.squared を求める箇所
+ r = z$residuals
+ f = z$fitted.values
+ if (attr(z$terms, "intercept")) { # 切片のあり/なし判定
+ # 切片ありのモデル 予測値と予測値の平均値からの偏差の二乗和
+ mss = sum((f - mean(f))^2)
+ } else {
+ # 切片なしのモデル 予測値の二乗和(予測値と0からの偏差の二乗和)
+ mss = sum(f^2)
+ }
+ rss = sum(r^2) # 残差の二乗和
+ r.squared = mss/(mss + rss) # 決定係数
+ c(r.squared = r.squared, mss = mss, rss = rss)
+ }
>
> lm1 = lm(Sepal.Length ~ ., iris) # 切片あり
> summary(lm1)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.17127 0.27979 7.760 1.43e-12
Sepal.Width 0.49589 0.08607 5.761 4.87e-08
Petal.Length 0.82924 0.06853 12.101 < 2e-16
Petal.Width -0.31516 0.15120 -2.084 0.03889
Speciesversicolor -0.72356 0.24017 -3.013 0.00306
Speciesvirginica -1.02350 0.33373 -3.067 0.00258

Residual standard error: 0.3068 on 144 degrees of freedom
Multiple R-squared: 0.8673, Adjusted R-squared: 0.8627
F-statistic: 188.3 on 5 and 144 DF, p-value: < 2.2e-16

> func(lm1)
r.squared mss rss
0.8673123 88.6118483 13.5564851

> lm0 = lm(Sepal.Length ~ . + 0, iris) # 切片なし
> summary(lm0)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
Sepal.Width 0.49589 0.08607 5.761 4.87e-08
Petal.Length 0.82924 0.06853 12.101 < 2e-16
Petal.Width -0.31516 0.15120 -2.084 0.03889
Speciessetosa 2.17127 0.27979 7.760 1.43e-12
Speciesversicolor 1.44770 0.28149 5.143 8.68e-07
Speciesvirginica 1.14777 0.35356 3.246 0.00145

Residual standard error: 0.3068 on 144 degrees of freedom
Multiple R-squared: 0.9974, Adjusted R-squared: 0.9973
F-statistic: 9224 on 6 and 144 DF, p-value: < 2.2e-16

> func(lm0)
r.squared mss rss
0.9974049 5210.2935149 13.5564851


Re: 重回帰での切片有り無しモデルについて
  投稿者:青木繁伸 2016/10/01(Sat) 20:58 No. 22157
うーーと。解説が必要なんでしょうか ?

切片なしのモデルの場合,R二乗は,回帰によって説明される,原点からのyの変動の割合を測定します

切片ありのモデルの場合,R二乗は,回帰によって説明される,平均値からのの変動の割合を測定します

Re: 重回帰での切片有り無しモデルについて
  投稿者:さくらとお城 2016/10/03(Mon) 10:11 No. 22158
青木先生

 説明など,ありがとうございます。理解できていませんが,定義が異なるということは,よくわかりました。

 ありがとうございました


Rを使った成績などの判定
  投稿者:初学者 2016/09/26(Mon) 20:19 No. 22150
お忙しいところ失礼致します。
Rについて、ご教示いただきたいと思いました。
データを集計していたのですが、例えば図のように数値と成績が対応している場合、各対象の数値から成績判定するには、どのようにコードを組めばよいのでしょうか。
初歩的なことで申し訳ございません。ご教示いただけますと幸いです。

よろしくお願いいたします。


Re: Rを使った成績などの判定
  投稿者:青木繁伸 2016/09/26(Mon) 21:24 No. 22151
やり方はいろいろありますが,区切りが不等間隔なこともあるので,一般的な方法を示しましょう。findInterval 関数を使います。
x <- seq(6.0, 10.5, by=0.1)
y <- 11-findInterval(x, c(0, 6.7, 6.9, 7.1, 7.3, 7.6, 8.0, 8.5, 9.1, 9.8))
data.frame(x, y)
findInerval の第2引数 c(0, 6.7, 6.9, 7.1, 7.3, 7.6, 8.0, 8.5, 9.1, 9.8)
0 以上 6.7未満を 1,6.7 以上 6.9 未満を 2,...,9.8 以上を 10 にするための区切りです。
11 からその値を引けば,欲しい数値になりますね。

Re: Rを使った成績などの判定
  投稿者:初学者 2016/09/26(Mon) 22:45 No. 22152
早速のお返事をいただき、誠にありがとうございます。
ご教示いただいた通り、findinterval関数を使って一度やってみます。

ありがとうございました。
今後とも何卒よろしくお願いいたします。


R 分割表から元データを再現
  投稿者:明石 2016/09/22(Thu) 19:56 No. 22147
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。

何卒どうぞ、よろしくお願いいたします。

−−−

青木先生の以下のサイトについて、ご質問がございます。

 http://aoki2.si.gunma-u.ac.jp/R/tenkai.html
 分割表から元データを再現 

使用例として、以下のように記載されています。

f <- matrix(c(2,3,1,4,3,5,4,3,2,1,2,4), nrow=3, byrow=TRUE)

result <- tenkai(f) # 結果はリストで返される
x <- result$x
y <- result$y

x
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3

y
[1] 1 1 2 2 2 3 4 4 4 4 1 1 1 2 2 2 2 2 3 3 3 3 4 4 4 1 1 2 3 3 4 4 4 4

この例題に、行名と列名を、以下のように与えたとします。
rownames(f) <- c("r1", "r2", "r3")
colnames(f) <- c("c1", "c2", "c3", "c4")

上記のx、yの値を、数字ではなく、行名と列名で返したいと思います。

そのRプログラムをご教示いただければ、大変に助かります。

どうぞ、よろしくお願いいたします。

Re: R 分割表から元データを再現
  投稿者:青木繁伸 2016/09/22(Thu) 21:15 No. 22148
せっかく与えた行名と列名を使うなら,
rownames(f)[x]
colnames(f)[y]
ただ,tenkai が何をやっているのかを見て,単に r, c で始まる名前でよいならば,行名,列名を与えてやる必要はなくて,
paste("r", rep(row(f), f), sep="")
paste("c", rep(col(f), f), sep="")
で十分。これを関数にしてやると面倒くささがなくなる。
paste0 が使えるならば,以下のように。
paste0("r", rep(row(f), f))
paste0("c", rep(col(f), f))

【御礼】 Re: R 分割表から元データを再現
  投稿者:明石 2016/09/23(Fri) 08:45 No. 22149
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

早々にご教示いただき、誠にありがとうございました。
今回も助けていただきました。

私も自分で、以下のような方法を考えてみました(途中で断念)。
・x,y からデータフレームを作成し、
・x,y が返す数字を、因子の水準番号と見なして、
 因子の水準に対応する名前(ここでが、行と列の名前)を与えて、
・最後に、数値を文字列表示に変換する

ご教示いただき、誠にありがとうございました。
心から御礼を申し上げます。


R ベクトル要素の完全一致検索
  投稿者:明石 2016/09/19(Mon) 15:17 No. 22144
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。

初歩的な質問で大変に恐縮ですが、何卒どうぞ、よろしくお願いいたします。

−−−

Rに組み込みの文字列ベクトル state.name(アメリカ50州の名前)があります。

state.name

"Alabama" "Alaska" "Arizona" "Arkansas" "California" "Colorado" "Connecticut" "Delaware" "Florida" "Georgia" "Hawaii" "Idaho" "Illinois" "Indiana" "Iowa" "Kansas" "Kentucky" "Louisiana" "Maine" "Maryland" "Massachusetts" "Michigan" "Minnesota" "Mississippi" "Missouri" "Montana" "Nebraska" "Nevada" "New Hampshire" "New Jersey" "New Mexico" "New York" "North Carolina" "North Dakota" "Ohio" "Oklahoma" "Oregon" "Pennsylvania" "Rhode Island" "South Carolina" "South Dakota" "Tennessee" "Texas" "Utah" "Vermont" "Virginia" "Washington" "West Virginia" "Wisconsin" "Wyoming"

state.name の50要素から、
文字列の長さが最大、最小となる要素(州の名前)を
州の名前、要素の番号を対にして、すべて取り出すことを考えます。

nagasa <- nchar(state.name)

文字列の長さが最大となる要素(州の名前)は、複数あり、
> state.name[ nagasa == max(nagasa) ]
[1] "North Carolina" "South Carolina"

文字列の長さが最小となる要素(州の名前)は、複数あり、
> state.name[ nagasa == min(nagasa) ]
[1] "Iowa" "Ohio" "Utah"

のように、州の名前はすべて取り出すことはできますが、
その要素の番号までは、とりだすことができません。

初歩的な質問で大変に恐縮ですが、何卒どうぞ、よろしくお願いいたします。

Re: R ベクトル要素の完全一致検索
  投稿者:青木繁伸 2016/09/19(Mon) 22:50 No. 22145
which を使えば良いでしょう

> which(nagasa==max(nagasa))
[1] 33 40
> state.name[which(nagasa==max(nagasa))]
[1] "North Carolina" "South Carolina"
> which(nagasa==min(nagasa))
[1] 15 35 44
> state.name[which(nagasa==min(nagasa))]
[1] "Iowa" "Ohio" "Utah"
> which(nagasa==12) # 任意の長さの州名を取り出す
[1] 34 38 39 41
> state.name[which(nagasa==12)]
[1] "North Dakota" "Pennsylvania" "Rhode Island" "South Dakota"

蛇足ではあるが,which.max, which.min では仕様にあわないが,それはそれで使い道がある
> which.max(nagasa)
[1] 33
> state.name[which.max(nagasa)]
[1] "North Carolina"
> which.min(nagasa)
[1] 15
> state.name[which.min(nagasa)]
[1] "Iowa"

【御礼】 Re: R ベクトル要素の完全一致検索
  投稿者:明石 2016/09/20(Tue) 08:32 No. 22146
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

which.max()関数、which.min()関数を使うことを考えましたが、
最初に検出した1つに要素のみを返すことから、
重複の場合も考慮した、which.max()関数、which.min()関数のような関数を
探しましたが、見つけることができませんでした。

青木先生からご教示いただきました内容で、
which.max()関数、which.min()関数の活用の幅が広がりました。

よく分かりました。

ありがとうございました。


母集団数の調整
  投稿者:蒲郡 2016/09/11(Sun) 19:55 No. 22139
お世話になります。

いくつかのイベントの発生率を地域間で比較するときに、母集団数の差を調整して計算しようと思います。

経験ベイズ法で計算した結果、すべての地域の発生率が同値になったイベントがありました。

この場合は「母集団数の差を調整する必要がない」と判断して、このイベントのみを非調整値として、他のイベント発生率とともに論じることは妥当なのでしょうか。
あるいは、すべてのイベントを非調整値で論じることが妥当なのでしょうか。

そもそも、先生方は、母集団数の差を調整する/しないは、どのように決めていますか。
ご教示いただけますと幸いです。

Re: 母集団数の調整
  投稿者:青木繁伸 2016/09/12(Mon) 17:05 No. 22140
母集団数とは,「母集団のサイズ」つまり「母集団の構成要素の個数」ですか?(そんな用語は統計学にはないと思いますが)。
例えば,47都道府県を47個の母集団とし,母集団のサイズとは各都道府県の人口みたいなことですか?

このような場合は,標本調査ではなく全数調査なので,母数(パラメータ)が観察されるので,検定とか推定は無意味でしょう?
つまり,ある県の例えば「○○の悪性新生物による死亡率」が人口10万対 12.345 ならそれは母比率ですということ。別の県で人口10万対 9.876 だとすれば,差は -2.469 で,それ以上でもそれ以下でもない。人口で調整する意味もない(性・年齢別人口で調整する意味はある--年齢調整死亡率)。

Re: 母集団数の調整
  投稿者:蒲郡 2016/09/12(Mon) 19:51 No. 22141
青木先生侍史

お忙しい中、ご教示をいただきありがとうございます。
用語を誤っており、お手数をおかけして失礼しました。
先生が推測して下さった通りの内容をお尋ねしたかった次第です。

特定の年齢別の人口で調整した場合に、年齢調整死亡率がすべての都道府県で同じ値になった場合は、「調整した結果、差はなかった」と考えて良いのでしょうか?

私の最初の質問を読みなおしてみますと「このイベントのみを非調整値として、他のイベント発生率とともに論じる」のは真実と異なっていると考えました。

Re: 母集団数の調整
  投稿者:青木繁伸 2016/09/13(Tue) 00:48 No. 22142
> 特定の年齢別の人口で調整した場合に、年齢調整死亡率がすべての都道府県で同じ値になった場合は、「調整した結果、差はなかった」と考えて良いのでしょうか?

例では都道府県としましたが,47 の年齢調整死亡率がすべて同じになった?そんなことあるのかなあ?

まあ,同じ値になったのなら「全く差はない」ということでいいでしょう。

> 「母集団数の差を調整する必要がない」

調整したから差がないことがわかったのでしょう?だったら,調整する必要はあるのでしょう。違いますか?

Re: 母集団数の調整
  投稿者:蒲郡 2016/09/14(Wed) 01:16 No. 22143
青木先生侍史

ご教示をいただきありがとうございます。
同じになることがあるのか??という疑問点については、ソフトウェアをRに替えて計算し直したところ、わずかですが差がありました。
ありがとうございました。とても勉強になりました。


比の検定
  投稿者:瑯崎 2016/09/07(Wed) 20:28 No. 22136
皆様平素よりお世話になります。下記のことを検定したい場合、どのような検定手法になるかご教示頂けないでしょうか。
ある同じ品種の植物を5個体育てて、地上部(茎や葉等)と地下部(根)の質量を測定した。地上部と地下部の質量比が一定であるを帰無仮説にして危険率5%で検定したいのです。
個体No 地上部質量 地下部質量
1 20.3 5.6
2 30.7 7.8
3 16.1 3.1
4 16.8 6.9
5 25 7.3

よろしくお願いします

Re: 比の検定
  投稿者:鈴木康弘 2016/09/09(Fri) 07:44 No. 22137
 思うに、方法のひとつは、質量比の値5つについて、その平均の信頼区間を示すこと。

http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean2.html

 もう一つの方法は、地上部と地下部の相関係数を出して、それが0でないことの検定。

http://aoki2.si.gunma-u.ac.jp/lecture/Corr/corr.html

 かな?

Re: 比の検定
  投稿者:瑯崎 2016/09/09(Fri) 08:37 No. 22138
ありがとうございます


R 文字列の連結
  投稿者:明石 2016/08/28(Sun) 21:02 No. 22125
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。
何卒どうぞ、よろしくお願いいたします。

−−−

文字列ベクトルがあります。
簡単な例示です。
"1", "12", "12345", "1234", "123"

各要素の文字列の長さが不揃いですが、
最大の長さ(ここでは、要素"12345"の長さ5)にあわせて、
必要に応じて、先頭に"0"を連結して、
すべての要素の文字列の長さをあわせることを考えます。

出来上がりは、以下のようになります。

"00001", "00012", "12345", "01234", "00123"

ベクトルの長さが膨大(数十万件)なことから、効率的な方法を検討しています。

以下の2点について、ご教示をいただければ、大変に助かります。
(1) 先頭に"0"を連結する方法
(2) applyを使う方法

どうぞ、よろしくお願いいたします。

Re: R 文字列の連結
  投稿者:青木繁伸 2016/08/29(Mon) 07:53 No. 22126
やりかたはいろいろあるけど,一番簡単(?)なのは sprintf を使うもの。
a = c("1", "12", "12345", "1234", "123")
n = max(nchar(a))
unname(sapply(a, function(b) sprintf("%0*s", n, b)))
* には n の値がつかわれ,実際には sprintf("%05s", b) と同じになる。
別の方法は
unname(sapply(a, function(b) paste(paste(rep("0", n-nchar(b)), collapse=""), b, sep="")))
のようなものがあるけど,見た目でもわかるが,sprintf を使うより実行時間は若干長い。
一番外側の unname は,文字列要素に名前が付くのをはぎ取るためなので,場合によっては不要。

【御礼】 Re: R 文字列の連結
  投稿者:明石 2016/08/29(Mon) 08:51 No. 22127
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

今回も助けていただきました。
心から御礼を申し上げます。

例えば、"00000"という文字列を作成するには、
paste( rep( "0", 5) , sep="", collapse="" )
ということまでは自分でも分かりましたので、
個々の要素の文字列の長さから、先頭につける文字列の長さの計算はできました。

ただし、数十万件をループで回している方法は何とも回避したいと思っておりました。

大変に良い勉強をさせていただきました。
ありがとうございました。

Re: R 文字列の連結
  投稿者:青木繁伸 2016/08/29(Mon) 11:27 No. 22128
sapply する必要はなかったですね。これだと unname も不要です。
a = c("1", "12", "12345", "1234", "123")
n = max(nchar(a))
sprintf("%0*s", n, a)

Re: R 文字列の連結
  投稿者: 2016/08/29(Mon) 14:57 No. 22129
青木先生

横から失礼します。
R-3.2.5(Linux)では"0"の代わりに半角スペースで埋められてしまいました。

?sprintfで調べると
## Platform-dependent bad example from qdapTools 1.0.0:
## may pad with spaces or zeroes.
sprintf("%09s", month.name)
のような例が載っていました。

unname(sapply(a, function(b) paste(paste(rep("0", n-nchar(b)), collapse=""), b, sep="")))
では0で埋まりました。

参考までに。

Special Thanks(Re: R 文字列の連結)
  投稿者:明石 2016/08/29(Mon) 20:07 No. 22130
青木先生、荒様

いつもお世話になり、ありがとうございます、明石と申します。

今回、大変に良い勉強をさせていただきました。

2つ、ご報告がございます。
・私は、Windows7、R-3.3.1の環境ですが、荒様と同じ結果となりました。
・100万件のデータに適用しましたが、爆走でした。

勇気を振り絞って、投稿させていただいて本当に良かった、と思いました。
ありがとうございました。

Re: R 文字列の連結
  投稿者:青木繁伸 2016/08/29(Mon) 21:34 No. 22131
Platform dependent ですか。

もとの文字列に空白が含まれないならば,
gsub(" ", "0", sprintf("%0*s", n, a))
の方が速いかもしれません。

Re: R 文字列の連結
  投稿者:明石 2016/08/30(Tue) 19:21 No. 22135
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

お示しをしてくださいました
gsub(" ", "0", sprintf("%0*s", n, a))

恐ろしいほど、早いです。

Rの文字列の連結は、paste()関数を使用することしか知りませんでしたので、
今回、sprintf()関数を教えていただきましたことを契機に、調べ物をしています。

大きな武器を得ました。
心から御礼を申し上げます。
ありがとうございました。


対応のある2要因の分散分析における交互作用についてのpost-hoc test
  投稿者:ポスドクR 2016/08/29(Mon) 22:28 No. 22132
青木先生、皆様
お世話になります。

現在、あるリハビリ介入に関する論文を投稿しているのですが、標記のことに関して査読者からの指摘に関して皆様のお知恵をお借りいたしたく投稿いたします。

今回の研究では、健常被験者20名に、Aという介入とBという介入を行い、その前後でのパラメータの変化(即時効果)を調べています。それぞれの被験者は、AまたはBの介入をランダムな順番で、十分な間隔(数週間以上)開けて行っています。
得られたデータとしては、各被験者でのA介入前、A介入後、B介入前、B介入後の4つです。
私たちは、介入の前後 と 介入の種類 の2要因での対応のある分散分析と考えて解析を起こったところ、それぞれの主効果は有意差なく、交互作用のみが有意な結果となりました。
我々はpost-hoc testとして介入前後、介入の種類のそれぞれの要因について、単純主効果を検討し、記載しました。(それぞれ前/後、A/Bしかないので多重比較は行わず、各条件におけるpaired t-testを行いました。)
すると、査読者から、「2要因の分散分析で有意な交互作用が得られた場合、すべてのパターンについての一元配置分散分析(A前、A後、B前、B後の4つについての分散分析)を行い、多重比較補正を用いて有意な組み合わせを検出することが必要である。」
との指摘を受けました。

いろいろな文献や資料を参考に調べてみましたが、査読者の指摘している内容が正しいというう確証が得られず、さりとて、それが誤っているとのはっきりした証拠も得られず、大変困っています。

post-hoc testに関する査読者の主張の妥当性について、ご意見をいただけましたら幸いに存じます。
よろしくおお願い致します。

Re: 対応のある2要因の分散分析における交互作用についてのpost-hoc test
  投稿者:青木繁伸 2016/08/29(Mon) 23:27 No. 22133
「2要因の分散分析で有意な交互作用が得られた場合」という長いキーワードでネット検索するといくつか出てきますが,そのような場合には「単純主効果の検定を行う」らしく,SPSS での分析法について説明してあるページもあるようですね。

「すべてのパターンについての一元配置分散分析(A前、A後、B前、B後の4つについての分散分析)」というのは,サンプルサイズの水増し(4倍,あなたのデータだとサンプルサイズが80!!)になるし対応のあるデータであるということについても考慮していないので,疑問だと思います。

Re: 対応のある2要因の分散分析における交互作用についてのpost-hoc test
  投稿者:ポスドクR 2016/08/30(Tue) 13:05 No. 22134
青木先生

早速のコメントありがとうございます。
やはり対応のあるデータで無理やり一元配置に持っていくのはちょっと問題ということですね。
心強いお返事ありがとうございました。

何とか差読者のご機嫌を損ねないようにやんわり返事を書いてみます。

ありがとうございました。


R merge()関数でのエラーメッセージへの対処方法
  投稿者:明石 2016/08/24(Wed) 22:04 No. 22122
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。
何卒どうぞ、よろしくお願いいたします。

−−−

merge()関数で、現在の私の知識、経験では理解できない事象に遭遇し、
途方にくれています。

2つの csvファイルを読み込み、データフレームに代入します。
変数名は、dat1、dat2 とします。

dat1 を左側、dat2 を右側として、左外部結合でmergeします。

結合キー"注文番号"のデータ型は、integer です。

dat3 <- merge(dat1, dat2, by.x="注文番号", by.y="注文番号", all.x = TRUE, all.y = FALSE)

これら2つの csvファイルを、EXCELのxlsxファイルとして保存し、
xlsx形式で読み込むと、なんと、上記のマージ処理で、エラーが発生します。

Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column

キー"注文番号"はユニークであることは、unique()関数で確認済みです。
重複はございません。

csvファイルならばよいが、xlsxファイルならばエラーになるという事象に途方にくれています。

EXCELのxlsxファイルの読み込みは、library(openxlsx) を使いましたが、
他のEXCELファイル読み込みのパッケージ(XLConnect、readxl)でも同様でした。

このエラーメッセージへの対処として、どこを確認すればよいのでしょうか?

ご助言をいただければ、助かります。
どうぞ、よろしくお願いします。

失礼いたします。

Re: R merge()関数でのエラーメッセージへの対処方法
  投稿者:青木繁伸 2016/08/25(Thu) 07:45 No. 22123
あなたのやったことと環境を再現できそうにないので(私のパソコンには EXcel はない),コメントはできません
CSV で読めるのなら,CSV で保存すればよいと思います
ほかの人とデータを共用するのでなければ,Excel を介在する意味はないと思います

【御礼】 Re: R merge()関数でのエラーメッセージへの対処方法
  投稿者:明石 2016/08/25(Thu) 12:14 No. 22124
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

貴重なコメントをありがとうございました。
御礼を申し上げます。

お手数をお掛けして、大変に申しわけございませんでした。
失礼をいたします。


就職に最も役立つ生物統計学
  投稿者:ひの 2016/08/23(Tue) 09:00 No. 22121
ずいぶんご無沙汰しております。質問というわけではありませんが、この記事によると、
http://forbesjapan.com/articles/detail/13256

2016年のランキングで「就職に最も役立つ修士号」に輝いたのは、生物統計学なのだそうです。日本ではどうなのかな?


分散ゼロの場合の平均値の比較
  投稿者:ヒロ 2016/08/18(Thu) 23:22 No. 22111
下記の2つの事例に対する判断の成否を教えてください。いずれものある生物のバイオマス増加量データに関するもので、A,B群(いずれも反復数は3)の平均値に有意差があるか知りたいというのがもともとの動機です。各々の判断は正しいでしょうか?

事例1
A群:20、20、25; B群:10、10、10
B群の分散がゼロであり、パラメトリックの検定はできない。ノンパラを行うにしては反復が少なすぎて有意差はそもそも出ない。よって、2群の間に有意差があるかの判断を保留し、反復数を増やして再試験を行う。

事例2
A群:20、20、25; B群:0、0、0
B群の平均値、分散がともにゼロであり、パラメトリックの検定はできない。ノンパラを行うにしては反復が少なすぎて有意差は出ない。しかし、検定はできないが、B群は明らかにバイオマス増加量がゼロなので、A群の方がバイオマスの増加が大きい、と判断する。

アドバイスをいただければ幸いです。
よろしくお願いいたします。

Re: 分散ゼロの場合の平均値の比較
  投稿者:青木繁伸 2016/08/19(Fri) 10:03 No. 22113
> B群の分散がゼロであり、パラメトリックの検定はできない。
> B群の平均値、分散がともにゼロであり、パラメトリックの検定はできない。

そんなことないでしょう?実際に計算してみましたか?
検定の計算式を見ても,どちらかの分散が0であっても,計算できる(0での割り算は発生しない)。たとえば,
http://aoki2.si.gunma-u.ac.jp/lecture/Average/t-test.html
にある式をみてみましょう。
R での計算結果も示しておきましょうか。
> t.test(c(20, 20, 25), c(10, 10, 10))

Welch Two Sample t-test # ウェルチの方法

data: c(20, 20, 25) and c(10, 10, 10)
t = 7, df = 2, p-value = 0.0198 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
4.495579 18.837755
sample estimates:
mean of x mean of y
21.66667 10.00000

> t.test(c(20, 20, 25), c(10, 10, 10), var.equal=TRUE)

Two Sample t-test # 等分散を仮定

data: c(20, 20, 25) and c(10, 10, 10)
t = 7, df = 4, p-value = 0.002192 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
7.039258 16.294075
sample estimates:
mean of x mean of y
21.66667 10.00000

> t.test(c(20, 20, 25), c(0, 0, 0))

Welch Two Sample t-test # ウェルチの方法

data: c(20, 20, 25) and c(0, 0, 0)
t = 13, df = 2, p-value = 0.005865 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
14.49558 28.83775
sample estimates:
mean of x mean of y
21.66667 0.00000

> t.test(c(20, 20, 25), c(0, 0, 0), var.equal=TRUE)

Two Sample t-test # 等分散を仮定

data: c(20, 20, 25) and c(0, 0, 0)
t = 13, df = 4, p-value = 0.000202 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
17.03926 26.29408
sample estimates:
mean of x mean of y
21.66667 0.00000
しかし,いずれの場合も,サンプルサイズが少ないことにクレームが付く可能性が大。

ついでに,

> B群は明らかにバイオマス増加量がゼロ

私は,「明らかにゼロ」とは思いませんが...

Re: 分散ゼロの場合の平均値の比較
  投稿者:ヒロ 2016/08/19(Fri) 11:41 No. 22114
青木先生
お返事ありがとうございます。分散ゼロの場合は検定できないと決めてかかっていました。すみません。計算式を見て納得しました。

関連してもう1点お尋ねしたいのですが、3群以上の多重比較を行う場合で、分散が他の群と異なるものがある場合、2群ずつウェルチのt検定を行い、ボンフェローニで補正したP値(3群の場合はP=0.05→P=0.017など)に照らして有意性を判断する、というやり方は適正でしょうか?

再びですが、よろしくお願いいたします。

Re: 分散ゼロの場合の平均値の比較
  投稿者:青木繁伸 2016/08/19(Fri) 21:11 No. 22116
分散が異なろうが同じであろうが,Welch の方法を使えば良いというのは,実例でも裏付けのあることです。

その上で,多重比較を行えばよいというのも,よく知られていることです。

Re: 分散ゼロの場合の平均値の比較
  投稿者: 2016/08/20(Sat) 09:19 No. 22119
3群以上の場合にはGames-Howell法を試してはどうでしょうか。

- http://aoki2.si.gunma-u.ac.jp/lecture/Average/Tukey1.html
- http://aoki2.si.gunma-u.ac.jp/R/tukey.html

Re: 分散ゼロの場合の平均値の比較
  投稿者:ヒロ 2016/08/20(Sat) 23:04 No. 22120
青木先生、荒様、
ご返信、ありがとうございます。

荒様ご紹介の方法は知りませんでした。試してみます。

今後もご指導、よろしくお願いいたします。


交互作用の検定について
  投稿者:立山 2016/08/16(Tue) 23:23 No. 22105
青木先生 始めまして。cox回帰分析のモデルでのinteractionについて質問させてください。

認知症の発症に関して解析したいのですが、要因A〜Dがあるとして、要因B〜Dで調整した要因Aによるハザード比を出すモデルは、説明変数に要因A〜要因Dを投入して行うことができました。

そして要因A*要因Bの交互作用があるのかも知りたいのですが、モデルに「要因A、要因B、要因A*要因B」の項を投入して要因A*要因Bが有意であるか調べればよいのでしょうか?それとも「要因A〜要因D」と「要因A*要因B」を投入して、要因C要因Dも調整したうえでの交互作用をみるのがよいのでしょうか。
手持ちの本など調べたのですが記載がなく困っております。

お忙しいところ恐縮ですがよろしくおねがいいたします。

Re: 交互作用の検定について
  投稿者:鈴木康弘 2016/08/19(Fri) 07:21 No. 22112
 青木先生からのコメントがないところをみると、お好きなように、ということかもしれませんね。

 主効果に「要因A〜要因D」の4つ入れると、「要因A*要因B」以外の交互作用はどうなんだ?と査読者に言われることはあるかもしれない。

Re: 交互作用の検定について
  投稿者:青木繁伸 2016/08/19(Fri) 21:19 No. 22117
いつも私の意見を期待されるというのも困った状況ですが...

これも,多変量解析は分析に関与するすべての変数に基づいて分析すべしという大原則からすれば,A,B,C,D の主効果,AB, AC, AD, BC, BD, CD, ABC, ABD, ACD, BCD, ABCD の交互作用からなるフルモデルをまず検討すべきでは?もっとも,普通は高次の交互作用の影響度は低いので
A,B,C,D,AB, AC, AD, BC, BD, CD
A,B,C,D,AB, AC, AD, BC, BD, CD, ABC, ABD, ACD, BCD
A,B,C,D,AB, AC, AD, BC, BD, CD, ABC, ABD, BCD, BCD, ABCD
の順でよいと思うのではあるけれども。

Re: 交互作用の検定について
  投稿者:立山 2016/08/20(Sat) 08:43 No. 22118
鈴木様 
アドバイスありがとうございます。
特にA×Bについて興味があったので調べてみようと考えたのですが、それ以外の交互作用についても検討をしたいと思います。

青木先生
アドバイスありがとうございます。
アドバイスを参考にさせていただき、優先度の高い順で検討したいと思います。


単変量を積み重ねても多変量にならない件
  投稿者:青木繁伸 2016/08/18(Thu) 18:53 No. 22108
以下で用いるデータセット zip 形式で圧縮,エンコーディングは utf-8

Download:22108.zip 22108.zip


Re: 単変量を積み重ねても多変量にならない件
  投稿者:青木繁伸 2016/08/18(Thu) 19:12 No. 22109
上のデータを使って分析した結果を吟味してみよう。
> d = read.table("glm.dat", header=TRUE, fileEncoding="euc-jp")
相関係数を見ると,y と x1, X2, x4 はかなりの相関があるが,x3 とはほとんど無相関(それぞれ散布図を描いてみるとよい)。
> cor(d)
x1 x2 x3 x4 y
x1 1.00000000 0.6130828 -0.08943241 0.4913103 0.5505327
x2 0.61308282 1.0000000 0.30814005 0.8306170 0.6977745
x3 -0.08943241 0.3081401 1.00000000 0.2534457 -0.0259393
x4 0.49131031 0.8306170 0.25344571 1.0000000 0.6663762
y 0.55053270 0.6977745 -0.02593930 0.6663762 1.0000000
x1 〜 x4 を用いて単変量ロジスティック回帰を行う(出力の一部を省略)
x1, x2, x4 が有意
> summary(glm(y ~ x1, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -11.33739 2.34670 -4.831 1.36e-06
x1 0.19621 0.04275 4.590 4.43e-06
Residual deviance: 80.772 on 98 degrees of freedom
> summary(glm(y ~ x2, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -19.6555 4.0133 -4.898 9.70e-07
x2 0.3542 0.0746 4.748 2.06e-06
Residual deviance: 51.258 on 98 degrees of freedom
> summary(glm(y ~ x3, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.65418 1.13867 -0.575 0.566
x3 -0.00582 0.02245 -0.259 0.795 有意ではない
Residual deviance: 118.52 on 98 degrees of freedom
> summary(glm(y ~ x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -14.13175 2.85497 -4.950 7.43e-07
x4 0.25013 0.05271 4.745 2.09e-06
Residual deviance: 65.143 on 98 degrees of freedom
x1,x2,x4 を使って分析すると
> summary(glm(y ~ x1+x2+x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -36.07419 9.16181 -3.937 8.24e-05
x1 0.22755 0.08223 2.767 0.00566
x2 0.27575 0.09711 2.839 0.00452
x4 0.15913 0.07817 2.036 0.04178
Residual deviance: 35.413 on 96 degrees of freedom AIC: 43.413
そらみろ,ちゃんと全部有意になったじゃないかと,思いきや
総当たり法でやってみて,AIC が小さい順に挙げると,x1, x2, x4 を含むモデルは 5 番目。
上位 4 モデルは「すべて x3 が入っている!!」
x3 って,y と相関なかったよね?
> print(all.logistic(d), sort.by="AIC")

deviance AIC Formula
26.56728 36.56728 y ~ x1 + x2 + x3 + x4
29.30262 37.30262 y ~ x1 + x2 + x3
33.52754 41.52754 y ~ x2 + x3 + x4
35.80734 41.80734 y ~ x2 + x3
35.41286 43.41286 y ~ x1 + x2 + x4
40.34653 46.34653 y ~ x1 + x2
43.94204 51.94204 y ~ x1 + x3 + x4
45.99639 51.99639 y ~ x1 + x4
48.08899 54.08899 y ~ x2 + x4
51.25798 55.25798 y ~ x2
58.07721 64.07721 y ~ x3 + x4
65.14314 69.14314 y ~ x4
80.77157 84.77157 y ~ x1
80.65916 86.65916 y ~ x1 + x3
118.52330 122.52330 y ~ x3
全部使うモデルはいくらなんでも
それに x4 は有意じゃないよ
> summary(glm(y ~ x1+x2+x3+x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -36.24623 10.78017 -3.362 0.000773
x1 0.22241 0.10053 2.212 0.026950
x2 0.49536 0.17206 2.879 0.003989
x3 -0.21480 0.09236 -2.326 0.020037
x4 0.14700 0.09586 1.533 0.125160
Residual deviance: 26.567 on 95 degrees of freedom AIC: 36.567
じゃあ,x4 なんか除いてしまえ。
いいんじゃないですか。
> summary(glm(y ~ x1+x2+x3, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -31.04391 9.05752 -3.427 0.000609
x1 0.20348 0.09389 2.167 0.030219
x2 0.56413 0.15877 3.553 0.000381
x3 -0.21993 0.08598 -2.558 0.010529
Residual deviance: 29.303 on 96 degrees of freedom AIC: 37.303
その他のモデルもやってみよう。
> summary(glm(y ~ x2+x3+x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -22.31070 5.92023 -3.769 0.000164
x2 0.49410 0.14808 3.337 0.000848
x3 -0.21307 0.07859 -2.711 0.006705
x4 0.11349 0.07827 1.450 0.147039
Residual deviance: 33.528 on 96 degrees of freedom AIC: 41.528
> summary(glm(y ~ x2+x3, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -19.33143 4.71465 -4.100 4.13e-05
x2 0.54806 0.13307 4.119 3.81e-05
x3 -0.21356 0.07522 -2.839 0.00452
Residual deviance: 35.807 on 97 degrees of freedom AIC: 41.807
じつは,このデータセット,他に x5 というのがあるのだけどこれが強力で,x5 一つだけで十分予測ができてしまうというシロモノ。分析に使った(使えた)変数だけを弄っているだけでは,そんな強力な変数があるなんて,お釈迦様でも気がつかないわけです。

Re: 単変量を積み重ねても多変量にならない件
  投稿者:青木繁伸 2016/08/18(Thu) 19:21 No. 22110
2群の判別分析の場合に,2群で平均値の全く同じ変数を使うと判別率が上がるというのがわかりやすく図示している
http://aoki2.si.gunma-u.ac.jp/lecture/Discriminant/ans2.html
分析例
http://aoki2.si.gunma-u.ac.jp/BlackBox/how2use/dataset2.html

つまり,2群の平均値の差が有意な変数を使って判別分析をするのは NG ということ

重回帰分析の場合も同様

Re: 単変量を積み重ねても多変量にならない件
  投稿者:さくらとお城 2016/08/19(Fri) 16:07 No. 22115
青木先生,いつもありがとうございます。大変わかりやすい例示で,よくわかりました。


単変量で有意差なし、多変量で有意
  投稿者:内科医8 2016/08/10(Wed) 10:03 No. 22096
青木繁伸先生、みなさま
上記に関して質問させて頂きます。
AはZに対して正の関連をするのではないかとの仮説に基づいた横断研究において、AとZとの関連を見るときに、単変量解析では有意な相関関係はありませんでしたが、多変量解析ではAはZに対して独立した正の関連因子でした。
 Aと同時に因子Bに関しても解析をしており、AとBは正の相関を認めるのに対し、BとZは負の相関を認めておりました。先の多変量解析にはBも含まれおり、BはZに対して独立した負の関連因子でした。Bの影響により、AとZの関連が単変量解析ではマスクされていたと考えました。
このような考えは統計学において、認められる考えなのでしょうか。また、このような現象を示す統計学用語などありますでしょうか。ご教授のほど、よろしくお願い申し上げます。

Re: 単変量で有意差なし、多変量で有意
  投稿者:鈴木康弘 2016/08/15(Mon) 07:02 No. 22100
 これまで青木先生が何度も書かれているように、単変量解析では有意な相関はないが、多変量解析では有意な相関がみられるというのは、何の不思議もありません。

 このような現象を示す統計学用語は?..さあ、わかりません。

Re: 単変量で有意差なし、多変量で有意
  投稿者:さくらとお城 2016/08/17(Wed) 11:22 No. 22106
横からすみません。

>これまで青木先生が何度も書かれているように、単変量解析では有意な相関はないが、多変量解析では有意な相関がみられるというのは、何の不思議もありません。

確かに青木先生が書かれているコメントは,私も何回か読んだことがあります。ですが,”どうして”と突っ込まれると,返事に困っている状況です。何かわかりやすい説明とか,具体的な例示がありませんか?

特に単回帰で相関の高い変数を選んでも,重回帰で必ずしも良くならない例示などです。どうぞよろしくお願いします。

Re: 単変量で有意差なし、多変量で有意
  投稿者:青木繁伸 2016/08/18(Thu) 18:46 No. 22107
別記事で


統計方法について、ご教示ください。
  投稿者:初心者 2016/08/15(Mon) 07:43 No. 22101
青木先生、いつもお世話になります。

初歩的な質問で、申し訳ありません。

私は、認知症発症の危険因子を調べるために、地域における65歳以上の高齢者を対象とした疫学調査を行いました。初回の疫学調査では、65歳以上の高齢者住民1000人のうち、650人を健常高齢者(認知症ではない高齢者)と診断し、高血圧や糖尿病、脂質異常症などの内科的疾患、うつ病などの既往について聴取しました。そして、その7年後に第2回目の疫学調査を行ったところ、初回調査で健常高齢者と診断した650人のうち、35人が認知症を発症していました。

そこで、ご質問があります。

―蕾鵑猟敢困之鮠鏐睥霄圓反巴任気譟2回目の調査では認知症と診断された人のうち、初回の調査時に高血圧や糖尿病、脂質異常症、うつ病などの既往がなくても、2回目の調査時までの7年間の間に、これらの疾患を発症された方が存在します(2回目の調査時に聴取しました)。これらの疾患が、認知症発症の危険因子となったかどうか調べるための統計方法として、初回調査時における疾患の存在の有無と2回目調査時における疾患の存在の有無について、おのおのにロジスティック回帰分析を行うという方法でよろしいのでしょうか?

認知症を発症した人が35人とサンプルサイズが小さいために、本当は危険因子となる疾患であるにもかかわらず、ロジスティック回帰分析で有意差がでない可能性も考えられると思います。その場合、effect sizeを求めたらいいのでしょうか?

初歩的な質問で申し訳ありませんが、何卒よろしくお願い致します。

Re: 統計方法について、ご教示ください。
  投稿者:青木繁伸 2016/08/15(Mon) 22:30 No. 22102
初回調査で健常高齢者と診断した650人を対象としてロジスティック回帰を行うということでしょう。それ以外ありますか?
有意差が出ないから effect size をもとめるということではないでしょう。排反じゃないでしょう。ロジスティック回帰分析をする,effect size も示すということでしょう?

Re: 統計方法について、ご教示ください。
  投稿者:初心者 2016/08/16(Tue) 13:16 No. 22104
青木先生、御教示ありがとうございました。
今一度、自分なりに勉強してみます。
これからもよろしくお願い申し上げます。

[1] [2] [3]

- J o y f u l  N o t e -
Modified by i s s o