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

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


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

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


lda関数の判別係数の再利用について
  投稿者:初学者 2019/04/14(Sun) 20:16 No. 22713
青木先生

こんにちは。
いつも掲示板で勉強させていただいております。

このたびは正準判別分析の判別係数を別の研究で再利用できないかと思い、
ご質問させていただきました。
Rのlda関数で分析しております。
print()で出力される判別係数を使って、別の新しいデータの分類に利用できれば、
と考えているのですが、そのようなことは可能なのでしょうか。
R上ではpredict()等で簡単に実現できると思うのですが、
Rが使いこなせない人のために、例えばEXCEL上で判別係数から新データの分類が
簡便にできないかと考えました。

ご教示のほど、どうぞよろしくお願いいたします。

Re: lda関数の判別係数の再利用について
  投稿者:青木繁伸 2019/04/15(Mon) 11:58 No. 22714
library(MASS)
MASS:::predict.lda
の 2 行を入力すると 89 行のソースプログラムが表示されます。
まるまる全部が必須というわけではないですが,その程度の処理が必要ということです。
それを Excel のワークシート機能だけを使って R を知らない人のために用意するというのはあまりお勧めできません。たとえ用意できたとして,使い方もある程度は複雑になるでしょう。データを用意してスタートボタンを押すだけというようにすることもできるでしょうが,かなりたいへんでしょう。それくらいなら predict(obj, newdata) と入力するのが遙かに簡単では?

Re: lda関数の判別係数の再利用について
  投稿者:初学者 2019/04/15(Mon) 12:28 No. 22715
青木先生

早速のご回答、誠にありがとうございました。
たしかに、おっしゃる通りだと感じました。

他の人に使っていただくときには、
・判別分析から学習したobjectをRDSファイル等にセーブ
・RDSファイル、実行するRのscript、簡単な解説、などを渡す
・predict(obj,newdata)してもらう
等の方法で検討してみたいと思いました。

大変勉強になりました。
ありがとうございました。


マンホイットニーのU検定〜同順位の補正〜
  投稿者:確率は難しい 2019/04/03(Wed) 14:29 No. 22708
青木先生 どうかよろしくお願いします。

マンホイットニーのU検定について、先生の解説
http://aoki2.si.gunma-u.ac.jp/lecture/Average/U-test.html
を見ると、同順位の場合に補正する数式が掲載されています。
この数式通りやってみたのですが、同順位が全体の中で多いような分布の場合に確率がうまく出ない場合があることに気がつきました。

例えば、満月の日の翌日は雨が降りやすいという仮説に対して、帰無仮説を立てるとして、
16日間観測して、満月の翌日だけに100ミリの雨が降り、あとは雨が降らなかったとします。
第1群 満月の翌日ではない日→雨量0のデータが15個
第2群 満月の翌日→雨量100のデータが1個
これを上記の同順位の数式も含めて計算してみると、確率が0.0001ということになりました。

しかし、満月の日の翌日は平均気温が低くなるということで考えると、小数点以下3桁ぐらいまで平均気温を計算して同順位がないとすると、
第1群 気温9.253、10.865度・・・・など、9度〜11度の日が15個
第2群 気温8.025度のデータが1個(最も低い)
これで計算してみると、確率は0.1ということになりました。

感覚的には、満月の日だけに雨が降ったとしても、確率が0.0001というのはおかしいと思いまして、後者の確率(0.1)の方が感覚に合っているように思います。

この、同順位の式を数学的に理解することは私には難しすぎると思いますが、それはさておき、

これを苦肉の策で解決するために、データに適当に小さい乱数(他のデータと逆転しないような十分に小さい乱数)を足して無理に順位をつけるような数字に置き換えて計算しました。例えば、0の日の雨量を、 0.143、0.054、0.658・・・・というふうに置き換えた計算です。そうすると、同順位がなくなり、気温の例と同じ確率になります。

このように、乱数をプラスして同順位をなくして解決する方が、同点の場合の数式を使うよりも実態に合っているように思うのですが、致命的な問題はありますでしょうか。
また、他に良い解決策はないでしょうか?
よろしくお願いします。

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:青木繁伸 2019/04/03(Wed) 18:03 No. 22709
stats の wilcox.test で,同順位があるために「タイがあるため、正確な p 値を計算することができません」となる場合は,coin パッケージの wilcox_test を使ってみましょう。
library(coin)
x <- c(rep(0, 15), 100)
g <- factor(c(rep("a", 15), "b"))
wilcox_test(x ~ g, distribution="exact")
以下が表示されます。
	Exact Wilcoxon-Mann-Whitney Test

data: x by g (a, b)
Z = -3.873, p-value = 0.0625
alternative hypothesis: true mu is not equal to 0
Z は,同順位補正をしたときの正規近似の Z 値で,表示されている p 値の算出には使われていません。

同時に表示されている p-value = 0.0625 は正確な p 値です。

coin パッケージの出現前には exactRankTests パッケージの wilcox.exact が使われていました。
library(exactRankTests)
wilcox.exact(x ~ g)
以下が表示されます。
	Exact Wilcoxon rank sum test

data: x by g
W = 0, p-value = 0.0625
alternative hypothesis: true mu is not equal to 0
W は wilcoxon 統計量です(どちらを1群にするかで,もう一方の極端な側の W = 15 が返されることもあります)。表示された p 値は,coin:::wilcox_test と同じ 0.0625 です。

もう少し他の道具でも試してみるには,
http://aoki2.si.gunma-u.ac.jp/R/index.html
の中の,
XV. exact 検定とモンテカルロ法による近似検定
3. マン・ホイットニーの U 検定(exact test)
に記載した R プログラム exact.mw を使ってみると以下のようになります。
exact.mw(x=rep(0, 15), y=1)
U = 0, Z = 3.87298, P 値 = 0.000107511
正確な P 値 = 0.0625
査察した分割表の個数は 2 個
示されている U = 0, Z = 3.87298 は,coin:::wilcox_test の Z と同じですね。

次の行に正確な p 値が示されています。0.0625 で,coin:::wilcox_test,exactRankTests:::wilcox.exact の p 値と同じです。

exact.mw は exact.mw(matrix(c(15, 0, 0, 1), 2)) のようにも書けますが,
http://aoki2.si.gunma-u.ac.jp/exact/utest/getpar.html
にオンライン版を用意してあります。
列数を 2 として,次の画面で
15 0
0 1
の2×2分割表を入力すれば, 0.0625 が表示されます。

ということで,正確な p 値は 0.0625 でよいようです。

====================

最後に,「乱数をプラスして同順位をなくして解決する方が、同点の場合の数式を使うよりも実態に合っているように思う」ということですが,乱数を加えて検定を行うという流儀は存在します。
しかし,今回の場合は
stats:::wilcox.test(1:15, 16)

Wilcoxon rank sum test

data: 1:15 and 16
W = 0, p-value = 0.125
alternative hypothesis: true location shift is not equal to 0
とすることと同じなので,わざわざ乱数を加えて同順位をなくすということがどの程度適切か疑問ですね。

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:確率は難しい 2019/04/04(Thu) 10:16 No. 22710
青木先生

アドバイスありがとうございました。
結局のところ、z値を用いて正規分布を用いるのは近似に過ぎないということに尽きるということだと理解しましたが、そういう理解で合ってますでしょうか?

実際に適用して統計分析をしており、同順位が少なく、データ数が多い場合には比較的正確な値に近く、有効に適用できていると思っています。たくさんのケースをエクセルで分析して、機械的にZ値を算出し、正規分布だと仮定してp値を出しているため、
ひとつひとつ統計ソフトに入力するというのは、出来れば避けたいと考えています。

この正確なp値の計算は、エクセルの表計算で単純に行えるようなものではないのでしょうか?

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:青木繁伸 2019/04/04(Thu) 12:19 No. 22711
それぞれのデータを wilcox_test(x ~ g, distribution="exact") で分析する数行のプログラムを書けばよいだけでしょう。

エクセルの表計算のほうが手間が掛かるのでは?

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:確率は難しい 2019/04/04(Thu) 15:28 No. 22712
先生 ありがとうございます。

タイムリーな回答をいただきまして、非常に助かりました。
本当にありがとうございました。

おそらく、「同順位が多くて、近似的なP値と正確なP値が大きく乖離するケース」はそれほど多くないと思いますので、その部分だけプログラムを書くように勉強してみます。


二項ロジスティック回帰分析
  投稿者: 2019/03/30(Sat) 02:00 No. 22701
青木先生

初めて質問をさせていただきます。
spssを用いて二項ロジスティック分析を行っていますが、どうにもわからなくて困っております。ご教示いただけないでしょうか。よろしくお願いいたします。

仮説に基づいて、独立変数群を5モデル作り、それぞれの影響をみるため順番に強制投入しています。以下の2つの方法を試みましたが、投入したモデルは同じなのに結果が違っています(B,標準偏差、wald,有意差、EXP(B))。なぜ違っているのか、また、どちらのやり方が適しているのかを教えていただけないでしょうか。

‘販変数投入の際にブロックごとに次々と指定していった場合
▲屮蹈奪ごとの指定せずに「総当たり法」を繰り返し、,里箸と同じブロックごとに変数を増やした場合

以上です。

お忙しいところ、申し訳ありませんが、どうか、よろしくお願いいたします。

Re: 二項ロジスティック回帰分析
  投稿者:青木繁伸 2019/03/30(Sat) 09:16 No. 22702
(2)の実際の分析過程があいまいなのではっきりしたことが言えませんが,一般的には,「分析に用いる変数セットにより,結果的に選択される変数セットが異なることはよくあること」

(1)の「独立変数投入の際にブロックごとに次々と指定」する場合は,投入順序により結果が異なることがある。増減法であろうと同じことが起きる。
増加法と減少法,増減法と減増法で結果が違うのもよくあること。

「ブロックごとの指定せずに「総当たり法」を繰り返し、(1)のときと同じブロックごとに変数を増やした場合」も同じ。ブロックごとに総当たり法で増やすべき変数を決めても,結局は(1)と同じように,投入順序で結果が違うこともある。

「分析に用いる変数セット」に依存しない方法は,全ての変数を対象に総当たり法でモデルを選ぶ。もっとも,理論的に解釈可能なモデルであるかはチェックする必要がある。

=======

以下のような思考実験をするとよく分かる。

以上のようにして「最良のモデル」が得られた後,ある独立変数を分析に加えたらどうなるだろうか?

その変数が橋にも棒にもかからないヘボ変数で,結果は変わらない。
その変数が,新たにモデルに加えられる。
その変数が,モデルに加えられる一方,今までにモデルに含まれていた変数が不適切とわかり除去される(変数1個とは限らない)
その変数が,モデルに加えられる一方,今までにモデルに含まれていた変数が不適切とわかり除去される(1個とは限らない),そして今までモデルに含まれていなかった変数がモデルに加えられる(1個とは限らない)。
その際,同じブロック内の変数に入れ替えが生じる。
また,ブロック内に適当な変数がなくなる。
さらに,ブロック内の複数の変数がモデルに含まれる。
などなど,全ての状況が考えられる。

だって,それが多変量解析だもの

Re: 二項ロジスティック回帰分析
  投稿者: 2019/03/30(Sat) 12:09 No. 22703
青木先生

早々に明快なご説明をいただきまして、誠に、ありがとうございます。

,皚△睚竸瑤療蠧順序は同じで、結果的にも変数は除去されなかったのですが、値(B,標準偏差、wald,有意差、EXP(B))が異なるのは、なぜだろうかと不思議に思っていました。
しかし、先生が明示してくださった「全ての変数を対象に総当たり法でモデルを選ぶ」方法で解析を進めようと思います。ありがとうございました。

今後も、どうぞよろしくお願い致します。

Re: 二項ロジスティック回帰分析
  投稿者:青木繁伸 2019/03/30(Sat) 16:58 No. 22704
文意がとりにくかったのですが,「最終的に採用された変数は同じ」ということですか?

データに欠損値があったりしますか?
2つの方法で,実際に使用されたサンプルサイズがいくつかわかりますか?
変数が同じでも,実際に分析に使われるサンプルサイズが異なれば,B,標準偏差、wald,有意差、EXP(B))が異なるのは当たり前です。

Re: 二項ロジスティック回帰分析
  投稿者: 2019/03/31(Sun) 23:27 No. 22705
青木先生

ご連絡ありがとうございます。
本当に、このように親身になって頂いてありがたいです。

ご質問して頂いてように「最終的に採用された変数は同じ」です。

そして、「実際に使用されたサンプルサイズ」」を確認すべく、もう一度、spssの結果のlog,鉢△鮓比べてみました。すると、青木先生の仰るようにサンプルサイズが違っていました!
,諒法だと投入した全ての変数の欠損値を含むので、△離汽鵐廛襯汽ぅ困茲蠅眈さくなっていました。
これが原因だったのですね。よくわかりました。

ちなみに、,任呂匹離皀妊襪Nの値が同じになりますが、△世肇皀妊襪瓦箸Nの値が違ってくるということになりますね。論文にモデル1、モデル2、モデル3、モデル4、モデル5と階層的に示したいときには、Nを揃えた,諒法が適切なのでしょうか?
何度も質問させていただき、恐縮ですが、ご教示いただけないでしょうか。

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

Re: 二項ロジスティック回帰分析
  投稿者:青木繁伸 2019/04/01(Mon) 08:29 No. 22706
モデルの評価を AIC などで行う場合は,サンプルサイズが同じでないと評価不能ですね。

ただし,対象者の属性により欠損値が生じやすくなるような状況がある場合は,結果の評価に注意が必要になるかもしれません。

なお,変数選択の結果を見て,モデルに残った変数だけを指定して再度分析を行うべきです。
以下のような状況を考察すれば分かるでしょう。
変数が a ~ e,対象者が 9 人。o は測定値がある(欠損値ではない),x は欠損値。
       a b c d e
case1 o o o o x
case2 o o o x x
case3 o o x o o
case4 o o o o o
case5 o o o o o
case6 o o o o o
case7 o o o o o
case8 o o o o o
case9 o o o o o
5変数全部を使って分析を始める場合,c, d, e が欠損値の case1, case2, case3 は分析対象から外されます。分析対象は 6 人です。
その後の変数選択で a, b, d が残ったとします(c, e は採用されなかった)。
これで終わってはいけません。
a, b, d を使って分析することにすると,欠損値を持つとして分析から外されるのは case2 のみとなります。分析対象は 8 人です。 case1 も case3 も分析対象となり,2 人増えます。

Re: 二項ロジスティック回帰分析
  投稿者: 2019/04/01(Mon) 10:56 No. 22707
青木先生

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

ようやく分かってまいりました。
Nを揃えるときの注意点に気をつけながら、モデルに残った変数でやり直しをしてみます。

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


Split plot + RCBD with repeated measures
  投稿者: 2019/03/25(Mon) 18:32 No. 22700
青木先生

自分で調べてみましたが、解決できなかったので質問させていただきます。
本当に申し訳ありませんが、よろしくお願いいたします。

10つの森林を対象に、それぞれの森林に種Aの樹木1個体と樹種Bの樹木1個体を植栽し、それらの個体の高さ(樹高)を2ヶ月おきに観察しています。
私が明らかにしたいのは、樹高に対する種の影響と経時的な変化(測定時期の影響)、それらの交互作用です。
森林をブロックとして捉え、その中に、処理区(樹種)を設置し、さらに、その経時変化を調べていますが、それぞれのブロックの処理区に個体が1つずつしか存在しません。
経時変化の測定はそれぞれの個体を繰り返し測定します(repeated measures)。
これを統計ソフトRで解析したいのですが、一般化線形混合モデルを用いる場合、変量効果をどのようにモデルに組み込めばいいのかが分かりません。
glmer(樹高~樹種+(1|森林/樹種/測定時期)+測定時期+樹種×測定時期)
などの形を考えていますが、これでいいのかも分かりません。
よろしければご教示いただけないでしょうか。

また、2つ目の質問になってしまいますが(不躾で本当に申し訳ありません)、上記の実験デザインの中で、それぞれのブロックの処理区に2つの個体を配置した場合はどのようなモデルになりますでしょうか。
こちらにつきましても、
glmer(樹高~樹種+(1|森林/樹種/個体/測定時期)+測定時期+樹種×測定時期)
などの形を考えていますが、よくわかっていません。
よろしければお教えいただけないでしょうか。

私の勉強不足なところをお尋ねしてしまい、大変申し訳ありません。
お手数をおかけいたしますが、どうぞよろしくお願いいたします。


最小選択肢にほとんど丸のつかない可能性のある尺度開発
  投稿者: 2019/03/11(Mon) 08:13 No. 22695
青木先生
初めて質問させていただきます。
初学者ですよろしくお願い致します。

現在、自己評価尺度の開発を行っています。尺度原案の信頼性妥当性検証を行っている
ところです。
尺度原案のリッカートを6件法「1.全く当てはまらない」、「2.当てはまらない」、「3.どちらかというと当てはまらない」、「4.どちらかというと当てはまる」、「5.当てはまる」、「6.非常に当てはまる」として調査を行いました。
返送されてきた調査票は、「4.どちらかというと当てはまる」「5.当てはまる」に回答が集まる傾向がありました。

このように最小選択肢にほとんど〇のつかない可能性のある尺度はどうなのでしょうか…?

(リッカートを変更する予定ですが…)

Re: 最小選択肢にほとんど丸のつかない可能性のある尺度開発
  投稿者:青木繁伸 2019/03/11(Mon) 09:20 No. 22696
もっと極端な場合,つまり一つの選択肢にほとんど全てが回答するという場合は,そのような質問は不適切ということになりますね。
そのような場合と比べてどうすべきかということでしょう。明確な基準というものはないと思います。他の研究者とのコンセンサスが得られるかどうかも判断基準になるでしょう。

Re: 最小選択肢にほとんど丸のつかない可能性のある尺度開発
  投稿者: 2019/03/11(Mon) 12:20 No. 22697
青木先生

早速お返事をいただきありがとうございました。
極端な場合は、質問項目として除外されることになると考えます。

他の研究者の意見も伺って検討していきたいと思います。ありがとうございました。


重回帰分析の基準変数
  投稿者:北海道 2019/03/08(Fri) 08:59 No. 22692
久しぶりにメールさせていただきます。よろしくお願いします。
重回帰分析の基準変数と言っていいのかわからないのですが、「不安を持ちやすい」という変数をリッカートで大変当てはまる、当てはまる、どちらかというと当てはまる、どちらかというと当てはまらない、当てはまらない、全く当てはまらないの6つで聞きました。最近は6件だとダミー変数を作らない方法でも良いようなので、他の変数と一緒にそのまま投入しspssで解析しました。全く当てはまらないを基準としてそれ以外のベーターが出てきたので、これを表にしました。
論文を投稿したところ、査読者が全く当てはまらないのベーターがない、といってきました。査読者にどのように説明したらいいかご教授願えますでしょうか

Re: 重回帰分析の基準変数
  投稿者:青木繁伸 2019/03/08(Fri) 09:35 No. 22693
「全くあてはまらない」のβは 0 なのです。
それが不都合だと思う人はその人のかって(統計知らずなだけ)なのですが。
あなた自身も,「カテゴリー変数を重回帰分析に使う」ということ,その際に「ある変数を基準とすること」について十分理解されていないようです(あなたはダミー変数を使ったつもりはなくても,解析プログラムはダミー変数を使っています)。それらについて書いてある本とか WEB ページで勉強すれば,あなたの言葉で説明できるようになるでしょう。

あまりよくないページが多いのですが,
http://jojoshin.hatenablog.com/entry/2016/05/08/005747
などは,結果の解釈法についても簡潔に書かれていると思います。

Re: 重回帰分析の基準変数
  投稿者:北海道 2019/03/08(Fri) 09:42 No. 22694
青木先生
早速のお返事ありがとうございました。
よく分かりました。
まだまだ未熟ですが、これからも研究を続けていきます。


Python
  投稿者:波音 2019/03/01(Fri) 17:31 No. 22691
質問でも意見でもなく恐縮ですが、、

何気に今気づいたのですが、青木先生はPythonのスクリプトを大量にまとめられていたのですね。

http://aoki2.si.gunma-u.ac.jp/Python/index.html

少し勉強しようと意気込みましたが、Rでできるからといって学ぶのを呆けておりました。
(先生には感服いたします)


簡単なサンプルサイズの求め方について
  投稿者:清水せつ子 2019/02/26(Tue) 17:59 No. 22686
ある特定の母集団の血液中の化学物質の濃度を求めるために実験を委託したいと考えており,
血液中の化学物質の濃度の90%タイル値を信頼水準95%で予測するためのサンプル数を計算しています.
計算はできる限りシンプルにしたいと思っています.

私は,過去のデータを整理して,血液中の化学物質の濃度を正規分布に従い,母分散が既知と仮定することができるから,母平均の検定に使われる検出力の式で求められるのではないかと考えています.

検出力=P(|z|>Zalpha/2)

ここで,Zは,標本平均を期待値と分散で標準化した値.alphaは有意水準.

そんな中,某研究仲間から,二項分布に従うと仮定すれば,一つのサンプルをとり,それが90%タイル以上となる確率は0.9,n個のサンプルのすべてが90%タイル以上となる確率は0.9^nだから,

1-0.9^n= (1-信頼水準)    (1)

という数式でサンプルサイズnが求まると言われました.

・その研究者に「血液中濃度を二項分布に仮定することはできないのではないか」と問いましたところ,不良品の計算と同じだから,0.9%タイル値を超えるか超えないかといった二値でみればいいのだと言われました.また,血液濃度の分布は,ノンパラメットリックで見ているという回答を得ました.

母集団の分布が変われば,90%タイル値を95%の信頼水準で予測するためのサンプルサイズは異なるはずですが,(1)式では母集団がどのような分布であっても,サンプルサイズは同じとなるのでおかしいと思います.

すなわち,この回答に様々な理論的な間違いが含まれていると思いますが,どこに間違いがあるか今のところ答えが見つかりません.

・さらに,二項分布で仮定して,90%タイル値を超えるか超えないかといった二値でみれるならば,
(1)式は 1-p(n)=(1-信頼水準),ここでp(n)はBin(n,p)ではないかと思います.

以上,不適切な質問かもしれませんが,以上の2点についてご助言いただけませんでしょうか.

なお,私がこのFAQで問いたいのは,上記の間違いを理論的に説明したく,そのご助言をいただきたいということです.参考のため,私が以前に作成した血液濃度のサンプルサイズに関するファイルを添付しようと思いましたが,ファイルサイズが大きく添付できませんでした.この件でまた悩みがでましたら,改めて投稿させていただきます.

Re: 簡単なサンプルサイズの求め方について
  投稿者:青木繁伸 2019/02/26(Tue) 18:50 No. 22687
利用上の注意に書いてあるように,

画像などのバイナリーファイルをアップロードすることが可能です。

添付可能ファイル : PNG, ZIP
最大投稿データ量 : 150 KB

ということです。


> それが90%タイル以上となる確率は0.9
  ? 0.1 では?

> 0.9%タイル値を超えるか超えないか
  ? 90%タイルの書き間違い?

検出力と信頼率(信頼度)を混同している?

Re: 簡単なサンプルサイズの求め方について
  投稿者:清水せつ子 2019/02/27(Wed) 11:17 No. 22688
青木先生

早速の回答ありがとうございます.
添付ファイルについては承知いたしました.
ご丁寧な情報をありがとうございます.

> それが90%タイル以上となる確率は0.9
  ? 0.1 では?

> 0.9%タイル値を超えるか超えないか
  ? 90%タイルの書き間違い?

はい,1-0.9,90%タイルが正しいです.

また,次の文章にも訂正があります.

>n個のサンプルのすべてが90%タイル以上となる確率は0.9^nだから,

0.9^n→(1-0.9^n)です.

どれも私の書き間違いです.
申し訳ありません.

・(1)式は検出力と信頼度を混合しているということですね.
そのとおりだと思いました.
まとめます.

検出力というのは,ある検定方法を用いて,ある値thetaに対する帰無仮説を棄却する確率です.
つまり,標本平均が母平均と同等という帰無仮説を立てた場合,標本平均が母平均と異なる対立仮説が示せる確率,ある標本の不良率が,母不良率と同等であるという帰無仮説を立てた場合は,それが異なる対立仮説が示せる確率を表すものです.

一方,信頼度は,帰無仮説が正しいとされる確率です.ここでは,95%信頼水準が当てはまります.

検定を行うことによる間違いが2種類あって,帰無仮説が正しいのに間違いだとする過誤を示すのが信頼度,対立仮説が正しいのにそれを間違いだとする確率を示すのが(1-検出力)です.

(1)式では,明らかに,検出力と信頼度が区別されていません.

・また,信頼度か検出力を求めようとしている式が,1-0.9^nというのも間違い.

もし混合していて,(1)式の左辺が検出力,右辺が1-信頼水準としているのであれば,

左辺は不良率と同じように検出力を求めることができ,1-0.9^nではなく,

1-P(n),

ここで,p(n)は二項分布(n,p),nは標本数,pは不良確率

となります.

以上,これが私の理解です.
間違いがあればご指摘ください.よろしくお願いいたします.

Re: 簡単なサンプルサイズの求め方について
  投稿者:青木繁伸 2019/02/27(Wed) 12:10 No. 22689
> 帰無仮説が正しいのに間違いだとする過誤を示すのが信頼度

違います。「帰無仮説が正しいのに間違いだとする過誤」は第一種の過誤(その過誤を犯す確率がα)。信頼度は 1-α で表される。同じようにしてデータを収集し,信頼限界を計算したとき,その信頼限界が母数(帰無仮説で仮定した統計量)を含む確率。信頼度をを大きくするにはサンプルサイズを大きくする必要がある。

> 検出力というのは,ある検定方法を用いて,ある値thetaに対する帰無仮説を棄却する確率です.

微妙だけど正確ではない。検出力とは,帰無仮説が間違っている場合に,帰無仮説を棄却できる確率(1-β)。βは,帰無仮説が間違っているのに帰無仮説を棄却できない確率。サンプルサイズを大きくすれば,検出力は大きくなる。

つまり,あなたの場合は,「血液中の化学物質の濃度の90%タイル値」を「誤差範囲内」で推定する際に,推定値「90%タイル値±誤差」が真値である「90%タイル値」を信頼度95%で含むようにするためにはサンプルサイズが幾つ必要かということでしょう。しかも,調査は1回しか行われないので,その調査が成功する為には検出力も高くする必要があるので,サンプルサイズはもっと必要であると。

母平均を求めるのならば http://aoki2.si.gunma-u.ac.jp/lecture/SampleSize/muconf.html ,特にその後半の方法が使えるのだが,90%タイル値を求めるのはこれではだめでしょうね。

解を求める直接的な方法はないかもしれない。モンテカルロ法で求めることはできるかも。

Re: 簡単なサンプルサイズの求め方について
  投稿者:清水せつ子 2019/02/27(Wed) 12:33 No. 22690
青木先生

早速のお返事ありがとうございます.

訂正ありがとうございます.
また,モンテカルロ法に対する助言についてもありがとうございます.
大変参考になります.
平均ではなく,90%タイル値の求め方については,少しやってみて,後日,躓きましたら再度投稿いたします.

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


独立変数の水準サイズについて
  投稿者:畠山 2019/02/20(Wed) 17:28 No. 22683
はじめて質問させて頂きます。

ロジスティック回帰分析で各独立変数の水準サイズはどの程度あれば良いのでしょうか?

「作業者」という独立変数が3つの水準(スタッフA、スタッフB、スタッフC)を持ち、
それぞれの数がスタッフA:250、スタッフB:500、スタッフC:1,000であった場合、得られた結果は各水準のサイズの差が大きく影響しますでしょうか?

あるreviewerから、この数の差が大きいので結果が信用ならないとのコメントを頂きました。
全体のサンプルサイズは独立変数×10以上であり、十分かと思うのですが、水準ごとの数についてわからず質問させて頂きました。

お忙しいところ恐縮ですが、宜しくお願い申し上げます。

Re: 独立変数の水準サイズについて
  投稿者:青木繁伸 2019/02/20(Wed) 19:52 No. 22684
査読者は,水準ごとのサンプル数が 2,4 倍とかなり異なることを言っているのでしょう

この違いが恣意的なものなら問題になるでしょうが,各水準での抽出割合を同じにして標本を抽出すると,結果として母集団の構成を反映することになります。このような標本抽出は,むしろ好ましいことでしょう。逆に母集団で 1:2:4 の割合で存在するものを,抽出割合を無視して 1:1:1 でサンプルを取って分析したらその方がよほど問題だと思いますが。

このようなことは事前に抽出数を決定しない変数には常に存在するものです。サンプルを無作為抽出した後,いろいろな変数ごとに集計すると水準ごとのサンプルサイズには相当な違いが出てくるでしょう。全ての変数で,ある要素を持っているものと持たないものの割合が 1:1 になることなんてないでしょう。変数によっては数倍から数十倍になるものもあるでしょう。

どうも,査読者は古典的な実験で,対照群と各処理群のサンプルサイズを全て同じにすべしという例数設定と同じだと誤解しているのではないか?

例えば,社会調査のような非実験系の統計で,平均値の差の検定のような単純な場合でも,データを無作為抽出した後にいろいろな属性で再分類して検定を行うような場合,例数がかなり異なることはよくあることです。「そのような場合に検定は行えないし,行ったとしても信用ならない」なんて,誰もいいませんよね。

平均値の差の検定だろうがロジスティック回帰分析だろうが,「水準ごとのサンプルサイズが違うからだめだ」ということにならないのは明らかでしょう。

Re: 独立変数の水準サイズについて
  投稿者:畠山 2019/02/21(Thu) 10:08 No. 22685
青木先生

ご返答ありがとうございます。
今回は臨床データを後方視的に検討したものだったので、恣意的なものではありませんでした。

各水準のサンプル数を意図的に変えたわけではない旨を上手く査読者の先生に伝えたいと思います。
お忙しいところご返答下さり誠にありがとうございました。


Rのreadlineによる入力待ちについて
  投稿者:桜とお城 2019/01/18(Fri) 13:02 No. 22677
 いつも参考にさせて頂いております。
  
 ある処理の結果を見て,変数に所定の数値をキーボードから入力してから,次の処理をするようにしたいので,以下の例をやってみました。

x=readline("input>")
x=as.numeric(x)
x*2

一行ずつ実行するならうまく行くけど,スクリプトで3行一括で実行すると,入力待ちになりません。

 基本的なことかもしれませんが,解決する方法がありませんか?どうぞよろしくお願い致します。

Re: Rのreadlineによる入力待ちについて
  投稿者:青木繁伸 2019/01/25(Fri) 19:32 No. 22681
長い記事の投稿があって,気づきませんでした。

解決法は以下のようにするとよいかと。

つまり,前後を { } で囲んで,複文にする。

{
x=readline("input>")
x=as.numeric(x)
x*2
}

ただ,これをそのままコンソールにコピーペーストしては何の改善もありません。

R のエディタにこれを書いて,この部分を選択して「実行」(Mac なら command+return)

> {
+ x=readline("input>")
+ x=as.numeric(x)
+ x*2
+ }
input>4
[1] 8

RStudio でも同様

Re: Rのreadlineによる入力待ちについて
  投稿者:桜とお城 2019/01/30(Wed) 10:58 No. 22682
青木先生,ありがとうございました。問題なくできました。


3 群以上の比率の差を対比較で検定する方法
  投稿者:Bob 2018/05/31(Thu) 21:56 No. 22559
青木先生
お世話になっております。

統計の初歩を学習している者でございます。
以下、お教え頂ければ幸甚に存じます。

Q1:二値尺度において、3群以上を比較する場合、その手法は、2群の比較と同様、フィッシャー法でよいと理解していますが、その理解であっていますでしょうか?

Q2:この検定で有意さがついたと前提です。その後個々の検定をするには何法を使用したらよろしいのでしょうか?
数値型でしたら、ダン法やチューキー法かと存じますが比率の差については存じません。

大変初歩的なことで恐縮でございます。種々調べましたがよく分かりませんでした。
お教え頂ければ幸いでございます。

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:青木繁伸 2018/06/02(Sat) 18:42 No. 22560
フィッシャーの正確検定を使うなら,ボンフェローニ法を適用するかな?

カイ二乗検定を使うならば
K 群の比率の差の検定・多重比較
ライアンの方法
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Ryan.html
テューキーの方法
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Tukey.html

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:Bob 2018/06/02(Sat) 19:59 No. 22561
青木先生

ご多忙の折、お教え頂き有難うございました。
とても助かりました。
今後もしっかり統計の勉強をしていきたいと思います。

取り急ぎお礼まで

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:渡辺 2019/01/18(Fri) 18:10 No. 22678
青木先生

古い投稿への相乗りという形で失礼いたします。

「K群 x 2カテゴリー」でなく、

「K x N」のときのやり方を教えていただけますと幸いです。

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:青木繁伸 2019/01/21(Mon) 21:37 No. 22679
そのようなものは見たことがないですが,「K の全ての分割法」×「N の全ての分割法」で検定して,その全ての可能性を n として,個々の検定を α/n でやるということになるのではないかと思います。まあ,ちょっと,現実離れしているかなあとは思います。

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:渡辺 2019/01/23(Wed) 11:22 No. 22680
遺伝学解析で比較したい遺伝子型が複数あり、表現型のクライテリアも4カテゴリーあるデータを得ています。
やはり個別にやって補正するしかないのですね。
2x4のカイ2乗検定を遺伝子型の各組み合わせでやって、組み合わせ数で割ることにします。

お答えいただきありがとうございます。


GLMによるモデルと信頼区間の図
  投稿者:ぐれいぷ 2019/01/17(Thu) 11:31 No. 22674
 一般化線形モデルによる解析による予測値と、その信頼区間の求め方についてお教え願います。

 7本の樹(A〜G)について、ある測定値(measure)と比率(30枚の葉のうち症状の出た葉(A)の割合)に関して、一般化線形モデル(ロジスティック回帰)で解析しています。元のデータとモデルによる曲線、さらに95%信頼区間を図示したいと考えています。実際のデータによるRでの解析は以下の通りです(信頼区間は下側のみ示します)。

 図示した場合、信頼区間の破線がx軸の値の細かさで変動し、4から5で0.001刻みにすると階段状になってしまうのですが、モデルによる曲線と信頼区間を示す適切な方法をご存知でしたらお教えください。よろしくお願いします。

> a
tree measure A N
1 A 4.50 11 30
2 B 4.33 11 30
3 C 4.33 13 30
4 D 4.50 17 30
5 E 4.67 23 30
6 F 4.33 23 30
7 G 5.00 23 30

>
> glm <- glm(cbind(a$A,a$N-a$A)~a$measure,family=binomial)

> glm$coefficients
(Intercept) a$measure
-8.270453 1.900956

> ##############x軸の値を0.1ごと
> x <- seq(4,5,by=0.1)
> pred.y <- 1/(1+exp(-(-8.270453+1.900956*x)))
> plot(a$measure,a$A/a$N,xlim=c(4,5),ylim=c(0,1))
> lines(x,pred.y)
>
> ##############信頼区間を求める
> mean <- length(a$N)/sum(1/a$N)
> mean
[1] 30
> lower <- qbinom(0.025,mean,pred.y)/mean
> lines(x,lower,lty=2)
>
>
> ##############x軸の値を0.001ごと
> x <- seq(4,5,by=0.001)
> pred.y <- 1/(1+exp(-(-8.270453+1.900956*x)))
> plot(a$measure,a$A/a$N,xlim=c(4,5),ylim=c(0,1))
> lines(x,pred.y)
>
> lower <- qbinom(0.025,mean,pred.y)/mean
> lines(x,lower,lty=2)

Re: GLMによるモデルと信頼区間の図
  投稿者:青木繁伸 2019/01/17(Thu) 14:05 No. 22675
qbinom(0.025,mean,pred.y)/mean は,qbinom(0.025,mean,pred.y) 自体が離散値なので,それを定数で割ったところで離散値には変わりない。したがって,「階段状に表示される」のが正確で正常な状態です。

Re: GLMによるモデルと信頼区間の図
  投稿者:ぐれいぷ 2019/01/17(Thu) 14:39 No. 22676
青木先生
 お忙しい中早速のご回答ありがとうございました。モデルと信頼区間について勉強しなおします。今後ともよろしくお願いします。


オフセット項を含む単回帰モデルのパラメタ推定
  投稿者:新春 2019/01/10(Thu) 22:28 No. 22670
オフセット項を含む単回帰モデルの切片の標準誤差の推定方法について,ご教示いただければ幸いです.
以下のようなRのスクリプトを実行したとき,切片の推定値と標準誤差をsummary関数で確認できますが,標準誤差(以下の例では0.1114)をどのように算出しているのかわからず悩んでいます.適切な数式や考え方などが掲載されている資料があれば,ご紹介いただきたくお願い申し上げます.

> n <- 300
> x <- rnorm(n, mean=10, sd=1)
> e <- rnorm(n, mean=0, sd=2)
> a <- 7
> y <- a + x + e
> d <- data.frame(x,y)
> lm.out <- lm(y~1+offset(x), data=d)
> summary(lm.out)

Call:
lm(formula = y ~ 1 + offset(x), data = d)

Residuals:
Min 1Q Median 3Q Max
-6.2107 -1.2351 -0.0889 1.3389 5.1611

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.9713 0.1114 62.59 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.929 on 299 degrees of freedom

> (a.hat <- mean(y) - mean(x))
[1] 6.971308
> y.hat <- a.hat + x
> (sigma.e <- sqrt(sum((y - y.hat)^2)/(n-1)))
[1] 1.929243

Re: オフセット項を含む単回帰モデルのパラメタ推定
  投稿者:新春 2019/01/13(Sun) 16:46 No. 22671
質問させていただいた件,誤差伝播の法則に基づく次の計算で,summary関数で表示される結果と一致する値が得られました.

> n <- 300
> x <- rnorm(n, mean=10, sd=1)
> e <- rnorm(n, mean=0, sd=2)
> a <- 7
> y <- a + x + e
> d <- data.frame(x,y)
> lm.out <- lm(y~1+offset(x), data=d)
> summary(lm.out)

Call:
lm(formula = y ~ 1 + offset(x), data = d)

Residuals:
Min 1Q Median 3Q Max
-5.4845 -1.3463 0.0459 1.3090 4.8227

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.0122 0.1123 62.42 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.946 on 299 degrees of freedom

> var.xbar <- sum((x - mean(x))^2)/(n*(n-1))
> var.ybar <- sum((y - mean(y))^2)/(n*(n-1))
> cov.xybar <- cov(x,y)/n
> (sigma.a <- sqrt(var.xbar + var.ybar - 2*cov.xybar))
[1] 0.1123454

以上,思いつきましたのでご報告いたします.
より適切な方法があれば,ご指摘いただければ幸いです.

Re: オフセット項を含む単回帰モデルのパラメタ推定
  投稿者:青木繁伸 2019/01/13(Sun) 21:30 No. 22672
基本的には,lm や lm.fit など,関連するソースを読んでいけば良いと言うことになります。

ご存じとは思いますが,ソースを見るのは,コンソールに lm などと入力するだけです。

また,実際にどのように解析が進むのかをみるのは debug(lm) などとした後,lm を実行するということでよいでしょう。debug についての事前の知識は必要でしょうが。

Re: オフセット項を含む単回帰モデルのパラメタ推定
  投稿者:新春 2019/01/14(Mon) 16:12 No. 22673
青木先生,

ソースコードの確認やdebugの活用など,ご指摘ありがとうございます.
丁寧にソースコードを追って,どのような計算が行われているのかを確認したいと思います.
ご助言いただけましたこと,また,このような質問の場を維持していただいていることに感謝申し上げます.


【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/02(Wed) 15:18 No. 22662
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
新年明けましておめでとうございます。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

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

文字列の部分的マッチング charmatch
http://www.okadajp.org/RWiki/?R%E3%81%AE%E6%96%87%E5%AD%97%E5%88%97%E5%87%A6%E7%90%86%E9%96%A2%E6%95%B0#b71cbd6d

以下の例題が示されています。
複数の場合は「0」を返すということも確認しました。

> charmatch("me", c("mean", "median", "mode"))
[1] 0

> charmatch("med", c("mean", "median", "mode"))
[1] 2

この2つを組み合わせて、以下も確認できました。
> charmatch(c("me","med"), c("mean", "median", "mode"))
[1] 0 2

さて、以下の例題を考えます。
v <- c("aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
s1 <- c("aa","bb")
s2 <- c("ba","ab")

charmatch(s1, v)
charmatch(s2, v)

> charmatch(s1, v)
[1] 1 2
> charmatch(s2, v)
[1] 4 3

複数の場合は「0」を返すという理解ですので、なぜ、「0」にならないのか、分かりません。
ご教示いただけましたら助かります。

新年早々、お手数をおかけいたします。
何卒どうぞ、よろしくお願いいたします。

Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:青木繁伸 2019/01/03(Thu) 15:08 No. 22663
一番最初の定義があるからでは?

> 正確なマッチが一部だけのマッチより優先される(つまり、ターゲットの先頭部分に 値が正確にマッチし、更にターゲットのほうがより長いケース)。
> charmatch(c("aa","bb"), c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
+ )
[1] 0 3
> charmatch(c("aa","bb"), c("aa", "aas", "bb", "ab", "ba", "abab", "abba", "baba", "baab"))
[1] 1 3
> charmatch(c("aa","bb"), c("aaw", "aas", "bb", "ab", "ba", "abab", "abba", "baba", "baab"))
[1] 0 3


Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/03(Thu) 16:24 No. 22664
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
新年明けましておめでとうございます。

青木先生がご指摘してくださいました箇所
> 正確なマッチが一部だけのマッチより優先される
> (つまり、ターゲットの先頭部分に 値が正確にマッチし、更にターゲットのほうがより長いケース)。

の意味を完全に理解できないでおりました。
大変に失礼いたしました。

---------------------------------------
青木先生、追加のご質問がございます。
何卒どうぞよろしくお願いいたします

> s <- c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
> v <- c("aa", "bb", "ab", "ba")
>
> grep(v[1], s)
[1] 1 2 9
> grep(v[2], s)
[1] 3 7
> grep(v[3], s)
[1] 4 6 7 8 9
> grep(v[4], s)
[1] 5 6 7 8 9

sを対象として、ベクトルvを、grep()関数で部分一致照合させた添字番号を取得したい
と思っていますが、
grep()関数では、第一引数は単体で与える必要があることから、
上記のように、vの要素ごとにgrep()関数を呼ぶことになります。

ベクトル v と s を与えて、上記の結果(4本のベクトル)をリストで得る
Rプログラムをご教示いただければ大変に勉強になります。

新年早々、お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:青木繁伸 2019/01/03(Thu) 18:32 No. 22665
sapply(v, grep, s) または lapply(v, grep, s) でできます
> s <- c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
> v <- c("aa", "bb", "ab", "ba")
> sapply(v, grep, s)
$aa
[1] 1 2 9

$bb
[1] 3 7

$ab
[1] 4 6 7 8 9

$ba
[1] 5 6 7 8 9

> lapply(v, grep, s)
[[1]]
[1] 1 2 9

[[2]]
[1] 3 7

[[3]]
[1] 4 6 7 8 9

[[4]]
[1] 5 6 7 8 9


Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/04(Fri) 08:19 No. 22666
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
新年早々に、大変に良い勉強をさせていただきました。
この例題を通じて、今まで不得手であった sapply、lapplyも理解できました。
心から御礼を申し上げます。

Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/04(Fri) 17:07 No. 22667
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

ベクトル v と s を与えて、部分一致照合の添字を得る方法のご教示をいただき、誠にありがとうございました。
大変に良い勉強をさせていただきました。

s <- c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
v <- c("aa", "bb", "ab", "ba")

lapply(v, grep, s)

[[1]]
[1] 1 2 9

[[2]]
[1] 3 7

[[3]]
[1] 4 6 7 8 9

[[4]]
[1] 5 6 7 8 9

取得した添字を用いて、要素を得ることができます。
> res <- lapply(v, grep, s)
> s[res[[1]]]
[1] "aa" "aa" "baab"
> s[res[[2]]]
[1] "bb" "abba"
> s[res[[3]]]
[1] "ab" "abab" "abba" "baba" "baab"
> s[res[[4]]]
[1] "ba" "abab" "abba" "baba" "baab"

---------------------------------------
青木先生、追加のご質問がございます。
何卒どうぞよろしくお願いいたします

ベクトル v と s を与えて、部分一致照合の要素
"aa" "aa" "baab"
"bb" "abba"
"ab" "abab" "abba" "baba" "baab"
"ba" "abab" "abba" "baba" "baab"
を得る方法をご教示いただければ、大変に勉強になります。

お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:青木繁伸 2019/01/04(Fri) 18:19 No. 22668
まあ,一行で書かなきゃいけないということでもないでしょうし,また,他に良い方法があるかもしれませんが,以下のようにすることでも出来なくはないです
> sapply(sapply(v, grep, s), function(i) s[i])
$aa
[1] "aa" "aa" "baab"

$bb
[1] "bb" "abba"

$ab
[1] "ab" "abab" "abba" "baba" "baab"

$ba
[1] "ba" "abab" "abba" "baba" "baab"

Re: 【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/04(Fri) 18:58 No. 22669
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

sapply(v, grep, s)で得た添字を、s[]に代入する試みをしていましたが、結局できませんでした。

ご教示いただきました関数を定義する方法を理解できました。

ありがたいご教示に、心から、心より感謝いたします。
ありがとうございました。


【R】青木先生のAHPプログラムの使用方法
  投稿者:明石 2018/12/20(Thu) 16:47 No. 22658
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

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

青木先生が作成されたAHP
http://aoki2.si.gunma-u.ac.jp/R/AHP.html

を使用するためのデータx y の与え方について、以下の例題でご確認をさせてください。

・評価基準(価格、品質、外見)
・代替案(商品A、商品B、商品C、商品D)
・一対比較表
 添付の画像ファイルにお示しをします。

# 各評価基準の重要度

x <- c(1/5, 1/7, 1/3)

# 各代替案の重要度
y <- matrix(c(1/3, 1/5, 1/3, 1/3, 3, 5, 7, 7, 3, 1/3, 1/7, 1/5, 1/3, 3, 5, 3, 5, 3), 6, 3)

x と y の与え方は、よろしいでしょうか?

ご確認させてください。
お手数をおかけいたします。


Re: 【R】青木先生のAHPプログラムの使用方法
  投稿者:青木繁伸 2018/12/21(Fri) 10:54 No. 22660
それでよいと思いますが

Re: 【R】青木先生のAHPプログラムの使用方法
  投稿者:明石 2018/12/21(Fri) 18:10 No. 22661
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

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

y <- matrix(c(1/3, 1/5, 1/3, 1/3, 3, 5, 7, 7, 3, 1/3, 1/7, 1/5, 1/3, 3, 5, 3, 5, 3), 6, 3)
で、行と列のサイズの与え方に不安がありましたので、ご確認をさせていただきました。

お手数をおかけいたしました、おかげさまで助かりました。

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

青木先生に、Rプログラムのご要望がございます。

意思決定分野のAHPとDEA(包絡分析法)について、利用できるRプログラムを調べております。

AHPについては、CRANのRパッケージがありますが、データ作成方法が煩雑で、苦慮していました。

AHPについては、青木先生のプログラムのおかげで、とても助かりました

DEAについても、CRANのRパッケージがありますが、マニュアルの理解で苦慮しています。

DEAについても、ぜひ、青木先生のRプログラムの作成、公開をお願いしたいと思います。

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

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


確証的因子分析のモデル適合について
  投稿者:さとう 2018/12/17(Mon) 22:33 No. 22655
はじめて質問させていただきます。
現在、確証的因子分析を行っています。
重みづけのない最小二乗法で実施しましたらGFI、AGFI、NFIが0.9以上でした。
RMSEAの値が出ないため。なんの値を見ていければ良いでしょうか。
教えていただけないでしょうか?
サンプル数は293 正規性の検定では、正規性はみられませんでした。
よろしくお願いいたします。

Re: 確証的因子分析のモデル適合について
  投稿者:波音 2018/12/20(Thu) 12:18 No. 22657
出力されているGFIやNFIを見て判断する、でいいと思いますが何か懸念されていることがあるのでしょうか、、

指標は挙げられているのも以外にもいくつかありますが、(個人の経験的に)どの指標を見ても何か結論が大きく変わることはないと思っています。

参考URL
http://daas.la.coocan.jp/sem/fitting_and_testing.htm

Re: 確証的因子分析のモデル適合について
  投稿者:さとう 2018/12/20(Thu) 22:52 No. 22659
コメントありがとうございます。
現在、共分散構造分析を学んでいる最中でして、学習不足で申し訳ありません。
指標の持つ意味を理解することが大切であり、どの指標を見ても結論は大きく変わらないということが、わかりました。質問して学びになりました。
ありがとうございました。


時間依存性共変量のpropensity score
  投稿者:徒弟 2018/12/18(Tue) 09:28 No. 22656
お世話になっております。

掲題の件で相談させて下さい。

ある患者コホートで、介入Xが、生存時間に影響するか評価したいと思ってます。
ただ、介入のタイミングが患者毎に異なっております。
すなわち介入は時間依存性共変量ということになります。

介入の有無を時間依存性共変量として、アウトカムを生存時間とするCox回帰で評価を行ってきましたが(結果は有意でした)、介入の有無で年齢等の患者背景が異なっており、propensity scoreマッチングのようなことをしたいと思ってます。

通常プロペンシティスコアは介入を従属変数としたロジスティック回帰で計算するのが一般的ですが、今回の場合介入が時間依存性共変量なので、Cox回帰でpropensity scoreを計算してはどうかと考えてます。

propensity scoreをCox回帰で計算するのは妥当でしょうか?
ご意見・ご指導賜れれば幸いです。


名義尺度の最小二乗フィッティングで追加制約がある場合の標準誤差の計算方法
  投稿者:小西 2018/12/05(Wed) 18:18 No. 22654
長文失礼致します。
以下についてお分かりになる方がいましたらご教授頂けますでしょうか?

現在、交互作用のない2因子名義尺度に対する最小二乗フィッティングを行うプログラムを作成しているのですが、特殊な制約を加えた場合の標準誤差の計算方法がわからず困っております。
目標は2種類の加工装置の各個体が製品のばらつきにどの程度寄与しているかを推定するため、平均値(切片)、ばらつきがどの個体に起因しているか(回帰係数)、それら係数の信頼区間を求めることです。
例えば装置Aが3水準(A1,A2,A3)、装置Bが4水準(B1,B2,B3,B4)あるとして、AとBの組み合わせごとに製品の測定値が与えられます。ただし、いくつかの組み合わせが欠損値となることもあり得ます。

基本的にはJMPの計算結果を参考にしているのですが、測定値の偶然性や欠損値などによってデータが線形従属になった(解に自由度がある)場合、JMPでは一部の水準を除外(係数を0に固定)する計算方法を採用しています。
この方法では入力データの水準の並び順が変わっただけで除外される水準が変化してしまいますし、特定の水準を特別扱いしていることにもなるので、作成中のプログラムではこの点を改良したいと考えています。
※JMPにおけるランク落ちの対処法の詳細は
 http://www.jmp.com/japan/support/help/13/flm-sls.shtml
 の下部にある「1次従属性がある場合」のサブトピックで説明されています。

そこで、特定の水準を特別扱いせずに解を一意に決定する方法として、各個体の回帰係数(切片は除く)の2乗和を最小にする制約を設けることを考えました。
(理論的根拠をもとに考えたわけではないので、統計的に不適切だったりより良い方法があるかもしれません)
そして、等式制約付き最小二乗問題を解くことで、そのような回帰係数を求めることもできました。
しかし、信頼区間を求めるのに必要な標準誤差をこの制約下でどのように計算すればよいかがわかりません。※入力データによっては自由度が不足して標準誤差を求められないケースがあることは承知しております。
この場合の標準誤差はどのように計算すればよいでしょうか?(下記詳細もご参照ください)

■参考:JMPと同等な計算方法
・入力データ例(製品測定値と使用装置) ※線形従属になる一例
 11 A1 B1
 13 A1 B2
 12 A2 B1
 15 A2 B2
 14 A3 B3
 10 A3 B4

 1列目の測定値ベクトルをy、データ数をnとします。
 ※この例では各組合せの測定値が高々1個になっていますが、実際には制限はありません。

・各因子水準をコード化(因子内の回帰係数の和は0に制約)
     A13 A23
 A1 →  1  0
 A2 →  0  1
 A3 →  -1  -1
 これを行列Aとします。

     B14 B24 B34
 B1 →  1  0  0
 B2 →  0  1  0
 B3 →  0  0  1
 B4 →  -1  -1  -1
 これを行列Bとします。

・コード化した入力データ
 切片 A13 A23 B14 B24 B34
  1  1  0  1  0  0
  1  1  0  0  1  0
  1  0  1  1  0  0
  1  0  1  0  1  0
  1  -1  -1  0  0  1
  1  -1  -1  0  -1  -1

・線形従属列の除外
 コード化した入力データの1列目、1〜2列目、1〜3列目、1〜4列目、1〜5列目、1〜6列目と順にランクを調べると6列目が増えたときにランク落ちしていることがわかるので、6列目を削除します。
 ※より大きな入力データの場合は後続の列も同様に繰り返します。

 切片 A13 A23 B14 B24
  1  1  0  1  0
  1  1  0  0  1
  1  0  1  1  0
  1  0  1  0  1
  1  -1  -1  0  0
  1  -1  -1  0  -1
 これを入力データ行列X、Xのランクをrとします。

・削減した係数を復元する行列
 ブロック行列
 (1 0 0)
 (0 A 0)
 (0 0 B)
 から上記で除外した線形従属列(この例ではB34列)を同様に除外した行列をZとします。

・最小二乗問題Xa=yを解く
 解ベクトルaから切片,A1,A2,B1,B2が決まり、Zaですべての係数が求まります。
 平均二乗誤差:MSE=(y'y-a'X'Xa)/(n-r)
 標準誤差:SE=sqrt(MSE*diag(Z*inv(X'X)*Z'))
 ※実際にはQR分解を利用してこれらと等価な計算を行います。

■詳細:検討中の2乗和を最小化する計算方法
※以下ではI_kはk次単位行列です。

・各因子水準のコード化(ここでは因子内の回帰係数の和を0に制約しません)
 A=I_3, B=I_4

・コード化した入力データ
 切片 A1 A2 A3 B1 B2 B3 B4
  1  1  0  0  1  0  0  0
  1  1  0  0  0  1  0  0
  1  0  1  0  1  0  0  0
  1  0  1  0  0  1  0  0
  1  0  0  1  0  0  1  0
  1  0  0  1  0  0  0  1
 これを入力データ行列X、Xのランクをrとします。

・等式制約付き最小二乗問題を解く
 Ca=0       ※切片以外の二乗和を最小化
 s.t. Da=0    ※因子内の回帰係数の和を0に制約
    X'Xa=X'y  ※正規方程式で最小二乗誤差に制約

 ここでCは以下のブロック行列
 (0 I_7)

 Dは以下の行列
 (0 1 1 1 0 0 0 0)
 (0 0 0 0 1 1 1 1)

 また、以下のブロック行列をEとします。
 ( D )
 (X'X)

 これにより解ベクトルaは求められますが、標準誤差SEの計算方法がわかりません。

 素人考えでJMPの計算式の形だけまねて、
 平均二乗誤差:MSE=(y'y-a'X'Xa)/(n-r)
 標準誤差:SE=sqrt(MSE*diag(pinv(X'X))) ※pinvは疑似逆行列
 としてみましたが、JMPで列の除外がない入力データの場合にJMPの結果と一致しませんでした。
 JMPの計算ではdiagの中の式(共分散行列?)にDa=0相当の制約条件が反映されているため、同様にDa=0相当およびEの零空間で二乗和を最小にする制約条件をdiagの中に反映する必要がある気がしています。


その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/01(Sat) 10:59 No. 22639
いろいろ枝葉末節にこだわりすぎた感がありました。
サンプルサイズが小さい場合に精度が低いのはやむを得ない(当然)なので,結局は,対数正規分布のパラメータとそれを対数変換した正規分布のパラメータの関係式から解を求めるしかないようですね。
対数正規分布の標本平均,標本分散から,非線形連立方程式で正規分布のパラメータを求め,それに基づき信頼限界を計算し,それを指数変換する。

元データなしで,信頼限界を求める
getCL = function(MEAN, VAR, n) {
func = function(p) {
return(c(exp(p[1]+p[2]/2)-MEAN, exp(2*p[1]+p[2])*(exp(p[2])-1)-VAR))
}
library(nleqslv)
ans = nleqslv(c(2,2), func) # 対数正規分布の平均値と分散から
mean = ans$x[1] # 対数をとったときの正規分布の平均値と
sd = sqrt(ans$x[2]) # 標準偏差を非線形連立方程式を解いて求める
cat(sprintf("mean = %f, sd = %g\n", mean, sd))
cl.normal = mean+c(1, -1)*qt(0.025, n-1)*sd/sqrt(n-1) # 正規分布目盛りでの信頼限界
cl = exp(cl.normal) # 対数目盛りに変換
cat(cl)
}
対数正規分布にしたがうデータから,標本平均・標本分散 MEAN = 3.390011 VAR = 37.060
が求められている(元データはない)

これから対数をとった正規分布の標本平均,標本標準偏差を求める
mean = 0.500339, sd = 1.20041
となる
それに基づいて,信頼区間を求めると
0.667087 4.077615

サンプルサイズが小さいと rlnorm で生成されたデータの平均値と分散は,母数とはかなりかけ離れたものになっている

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:Saito 2018/12/03(Mon) 10:29 No. 22642
コメントありがとうございます。
しかし、一つ前のスレッドでも書きましたが、こちらの手元にあるのは平均値とその標準「誤差」であり、元データの平均値と標準「偏差」ではありませんでした。そのため、微妙に私のやりたいこととずれていることに気づきましたので、再度ご教授頂ければ幸いです。何度も何度も申し訳ありません。

まず、元データの平均値と標準偏差が使える場合に合うことを示します。

> set.seed(1)
> n <- 100000
> u <- 1
> sig <- 0.5
>
> #生物の密度データがあるとする(非負)。平均密度とその信頼区間を知りたい。
> y <- exp(rnorm(n, u, sig))
> hist(y)
>
> #yというデータがあるもとで、正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
>
> res2 <- summary(lm(log(y)~1))
> lo2 <- exp(res2$coef[1]-1.96*res2$coef[2])
> up2 <- exp(res2$coef[1]+1.96*res2$coef[2])
>
> MEAN = mean(y)
> VAR = var(y)
> lo2;up2
[1] 2.706802
[1] 2.723691
> getCL(MEAN, VAR, n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906>

上記のように、MEAN=mean(y), VAR=var(y)として、元データyの平均値と分散が使えるという状況であれば、getCLとlo2, up2がほぼ一致することを確認しました(なお、sigが大きいと合わないようなので小さくし、小サンプルであることの問題も省きたいため、サンプルサイズnは大きくしてあります)

ところが、私の手元にあるのは、yの平均値と標準偏差ではなく、解析結果であるyの平均値とその標準誤差となっています。つまり、res$coef[1]とres$coef[2]が手元にあるということでして、mean(y)とvar(y)の値が手元にあるわけではありませんでした・・・。私の説明が拙かったようで、本当に申し訳ありません。
確認してみましたが、当然ながら、こちらとは合いません。

> MEAN = res$coef[1]
> VAR = res$coef[2]^2
> lo2;up2
[1] 2.706802
[1] 2.723691
> getCL(MEAN, VAR, n)
mean = 1.124605, sd = 0.0016874
3.07897 3.079034>

何度もお尋ねするのは心苦しいのですが、res$coef[1]とres$coef[2]が手元にあり、mean(y)とvar(y)の値(とn)が手元にない条件で、lo2, up2に準じるものは出せるでしょうか。何度も申し訳ありません。

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

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/03(Mon) 11:23 No. 22643
> 私の手元にあるのは、yの平均値と標準偏差ではなく、解析結果であるyの平均値とその標準誤差となっています

標準誤差=標準偏差/sqrt(サンプルサイズ) ですけど?
標準偏差=標準誤差*sqrt(サンプルサイズ)
分散=標準偏差^2

でしょう?

> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
lm する必要なんてなくて,
res$coef[1] = mean(y)
res$coef[2] = sd(y)/sqrt(n)
ですよね。

簡単なところを複雑にやっているので,頭が混乱します

> set.seed(1)
> n <- 100000
> u <- 1
> sig <- 0.5
> #生物の密度データがあるとする(非負)。平均密度とその信頼区間を知りたい。
> y <- exp(rnorm(n, u, sig))
> hist(y)
> #yというデータがあるもとで、正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
> lo1
[1] 3.068823
> up1
[1] 3.089189
> res$coef[1]
[1] 3.079006
> mean(y)
[1] 3.079006
> res$coef[2]
[1] 0.005195524
> sd(y)/sqrt(n)
[1] 0.005195524
> (res$coef[2]*sqrt(n))^2
[1] 2.699347
> var(y)
[1] 2.699347
>

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/03(Mon) 11:39 No. 22644
つづき

> mean(y)
[1] 3.079006
> var(y)
[1] 2.699347
> n
[1] 1e+05

以下で mean, var を使っているけど yは必要ない。平均値と分散を引数として与えているだけ。

> getCL(mean(y), var(y), n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906

以下のように必要なのは数値だけ
ただ,ここで示したように有効数字7桁では正確ではないので,コンピュータ上で持っている正確な値を使っただけ。
あなたは,y は持っていないが,計算された平均値と分散(をサンプルサイズで割ったものの平方根である標準偏差)を持っているので,その二つの数値を指定するだけ。

> getCL(3.079006, 2.699347, 100000)
mean = 0.999331, sd = 0.500551
2.70805 2.724906

標準誤差と分散の関係をちゃんとすれば,res$coef[1], res$coef[2] を使っても,当然同じ結果になる

> res$coef[1]
[1] 3.079006
> (res$coef[2]*sqrt(n))^2
[1] 2.699347


> getCL(res$coef[1], (res$coef[2]*sqrt(n))^2, n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906 当然ながら同じになる

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/03(Mon) 13:49 No. 22647
また枝葉末節なこと(基本的なこと)になりますが,あなたが使っている 1.96 は標準正規分布由来のものだが,本来は t 分布に基づく値でなければならない。

http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean2.html
の「母分散が不明の場合」にあたるので。

いまはテストデータでやっているから n がとても大きく,1.96 で問題ないが,実際に運用するときには問題になる。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:Saito 2018/12/04(Tue) 08:56 No. 22650
コメントありがとうございます。
不勉強で申し訳ありません。理解できていないのですが、標準誤差から標準偏差に戻すときに、nが未知でも戻るのでしょうか?今手元にあるのは、平均値とその標準誤差のみで、サンプルサイズnは不明なのです。getCLの中身を見る限り、何だか無理なような気がしてきています・・・。
もし戻らないとなると、当初の目的を果たすことは不可能でしょうか。

t分布については失念していました。
n=10だと、2.228ですね。ありがとうございます。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/04(Tue) 09:38 No. 22652
> 標準誤差から標準偏差に戻すときに、nが未知でも戻るのでしょうか?

定義式に書かれているとおり。
nが分からないなら無理でしょう。
現代統計学では,サンプルサイズこそ一番大事なものです。

妥当と思われるnを仮定するか,何通りかのnで代替してこの範囲という結果を示すしかないでしょう。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:Saito 2018/12/04(Tue) 11:25 No. 22653
コメントありがとうございます。
nがわからなければ不可能という旨、承知いたしました。
色々とご教授くださりありがとうございました。

以下余談です。
何故こうなるのかわからず、おそらく間違っていると思って相談しませんでしたが、この件について私ではない方(前任者)のメモで

Lower=exp(log(res$coef[1])-1.96*res$coef[2]/res$coef[1])
Upper=exp(log(res$coef[1])+1.96*res$coef[2]/res$coef[1])

として、平均値を対数を取って、なぜかCVを足す(引く)というやり方で正規分布を仮定して得られた平均値と標準誤差(res$coef[1], res$coef[2])を使って、対数正規分布を仮定したときの95%信頼区間(Lower, Upper)になるというメモがありました。やってみても、getCLと全然違う範囲になりますし、CVという単位が消えた値を足すのも奇妙なので、間違っているのでは、と思っています。とはいえ、他の用途でも構いませんので、この式に見覚えがありましたら、ご教授くだされば幸いです。


片側検定の帰無仮説
  投稿者: 2018/12/03(Mon) 09:19 No. 22640
例えば、母比率の差の検定で、
対立仮説が「A≠B」の場合、
帰無仮説は、両側検定の場合、「A=B」と認識しています。
対立仮説が「A>B」の場合、
片側検定の場合の帰無仮説は、
「A=B」でよいのでしょうか?
厳密には「A≦B」と書くべきではないでしょうか?

Re: 片側検定の帰無仮説
  投稿者:青木繁伸 2018/12/03(Mon) 10:25 No. 22641
検定統計量は,帰無仮説が正しいつまり,母数が特定の一つの値であると仮定したときのもの
片側検定でも,帰無仮説は「母数=特定の値」
http://aoki2.si.gunma-u.ac.jp/lecture/Kentei/p-value.html
の図で,「H0の下での統計量の分布」というのは,たとえば A=B のときに計算できるもの。A≦B というようなとき,分布は特定できないでしょう?

片側検定でも,両側検定でも「H0の下での統計量の分布」は同じものです。
どの部分の面積(確率)をP値とするかが違うだけです。

Re: 片側検定の帰無仮説
  投稿者: 2018/12/03(Mon) 12:06 No. 22645
ご回答ありがとうございます。

同様の質問になりますが、
母比率の差の検定のp値の意味についてご教授頂けましたら幸いです。

両側検定も片側検定も、
実際には差がないのに、差があると誤ってしまう確率
ということになるのかと思うのですが、
片側検定の場合、例えば、
実際にはA≧Bなのに、A<Bと誤ってしまう確率
とすることは不適当なのでしょうか?

Re: 片側検定の帰無仮説
  投稿者:青木繁伸 2018/12/03(Mon) 13:15 No. 22646
検定の基本的なところなので,正確に述べる必要があります。
「似たようなこと」をいっているのでよいだろうと言うことにはなりません。
前提と叙述の順序と用語法がとてもだいじ。
======
P 値とは,「帰無仮説のもとで,観察された統計量より極端な値を取る確率」

http://aoki2.si.gunma-u.ac.jp/lecture/Kentei/p-value.html において

両側検定の場合は PL, PU の和

片側検定の場合は,対立仮説の方向によりPL か PU のいずれか

もし,P 値が(たとえば)0.05 より小さい場合は帰無仮説を棄却する。

両側検定の場合は,帰無仮説はたとえば 母比率=0.5 で,対立仮説は 母比率≠0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率≠0.5 と結論する。
しかし,帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P

片側検定の場合は,帰無仮説は同じく 母比率=0.5 で,対立仮説はたとえば 母比率>0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率>0.5 と結論する。
しかし,帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P

最後の行で,「帰無仮説が正しい(母比率≦0.5)」ではないことに留意

Re: 片側検定の帰無仮説
  投稿者: 2018/12/03(Mon) 15:39 No. 22648
ご教授頂きまして誠にありがとうございます。
少しずつわかった気がいてきました。

両側検定の場合は,帰無仮説はたとえば 母比率=0.5 で,対立仮説は 母比率≠0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率≠0.5 と結論する。

帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
これを
帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
としてもよいのでしょうか?

同様に

片側検定の場合は,帰無仮説は同じく 母比率=0.5 で,対立仮説はたとえば 母比率>0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率>0.5 と結論する。
帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
これを
帰無仮説が正しい(母比率=0.5)のに母比率>0.5 とする確率は P
としてもよいのでしょうか?
それとも
帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
でしょうか?

Re: 片側検定の帰無仮説
  投稿者:青木繁伸 2018/12/03(Mon) 17:16 No. 22649
> 帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
> これを
> 帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
> としてもよいのでしょうか?

「帰無仮説を棄却する」は
「対立仮説を採択する」かつ「対立仮説は母比率≠0.5」なので「母比率≠0.5とする」
と同じですから,よいでしょう。

> 帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
> これを
> 帰無仮説が正しい(母比率=0.5)のに母比率>0.5 とする確率は P
> としてもよいのでしょうか?

両側検定の場合と同じです。よいでしょう。

> それとも
> 帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
> でしょうか?

「母比率>0.5 とする確率」と「母比率≠0.5とする確率」は違うので,後者は間違い。
そもそも,方向性のある片側対立仮説との比較なのに,なぜ両側仮説が出てくるのか?
むりやり考えると
「母比率>0.5 とする確率」と,逆側の「母比率<0.5 とする確率」を合計したものが 2P ということか。
つまり,片側対立仮説「母比率>0.5」も,「「母比率<0.5」も 統計学的には帰無仮説「母比率=0」から見れば同じく棄却域にあるので,下側確率と両側確率の和だが,今の場合分布は対称なので確率は 2P

Re: 片側検定の帰無仮説
  投稿者: 2018/12/04(Tue) 09:28 No. 22651
ありがとうございました。
だんだん理解できました。


正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法
  投稿者:Saito 2018/11/29(Thu) 17:11 No. 22635
お世話になっております。
標記の通りなのですが、生データがなく正規分布を仮定した際の平均値と標準偏差しかない状況で、対数正規分布を仮定していた場合の信頼区間を算出する方法について、私のコーディングが間違っているだけかもしれませんが、うまく推定できないので、どなたかご教授頂けると幸いです。

下記に例を示します。

> set.seed(1)
> n <- 10000
> u <- 0.5
> sig <- 0.8
> x <- exp(rnorm(n, u, sig))
>
> #1.平均値
> mean(x)
[1] 2.269916
>
> #2.標準偏差
> sd(x)
[1] 2.123908
>
> #3.信頼区間下限値
> mean(x)-1.96*sd(x)
[1] -1.892943
>
> #4.信頼区間上限値
> mean(x)+1.96*sd(x)
[1] 6.432775

このとき、手元にあるのは、1-2の平均値、標準偏差、そして95%信頼区間(の近似値)である、3-4の値のみで、生データであるxは利用不可能な状況であるとします。

上記の信頼区間の算出方法ですと、xが負の値を取りえないのに、マイナスの値をとってしまっていることが確認できると思います。この部分を正の値の下限値を出すようにしたいと考えています。

一応、対数正規分布の信頼区間から正規分布の信頼区間を得る方法が使えるのかと思い、下記のURLにある計算式を打ち込んでみますと、実際の95%点とずれが出るように思われます。
http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf

#推定値
> mean(x)/exp(1.96*sqrt(log(1+(sd(x)/mean(x))^2)))
[1] 0.4797148
> mean(x)*exp(1.96*sqrt(log(1+(sd(x)/mean(x))^2)))
[1] 10.7408

#実際の95%点
> quantile(x, c(0.025, 0.975))
2.5% 97.5%
0.3267348 8.0988751

たぶんこういう風に逆算でやるのはダメ、ということなのだと思うのですが・・・。
試しにlogの世界でやってみるとおおよそ合います。
> quantile(log(x), c(0.025, 0.975))
2.5% 97.5%
-1.523258 2.489656
> mean(log(x)-1.96*sd(log(x)))
[1] -1.490756
> mean(log(x)+1.96*sd(log(x)))
[1] 2.477682

どこかで間違っているのだと思いますが、どなたかmean(x)とsd(x)だけを使って、quantile(x, c(0.025, 0.975))で出力される95%点(負の値を取らない)を出すやり方をご存知の方がおられましたら、ご教授ください。

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

Re: 正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法
  投稿者:青木繁伸 2018/11/29(Thu) 19:46 No. 22636
いろいろあるので,順を追って

set.seed(1)
n <- 10000
u <- 0.5
sig <- 0.8
としたとき,
x0 <- rnorm(n, u, sig)
の x0 の平均値と標準偏差は当然ながら 0.5 と 0.8 ではないということ

> x0 <- rnorm(n, u, sig)
> mean(x0)
[1] 0.4947704
> sd(x0)
[1] 0.8098852

正確に指定された平均値と標準差を持つ正規乱数を生成するのは,一度正規化した後に,標準偏差を掛けて平均値を加える(標準化の逆をやる)。

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8

この指数を取って対数正規分布にする。

> x2 <- exp(x1)
> mean(x2)
[1] 2.268486
> sd(x2)
[1] 2.130959

しかし,平均値と標準偏差を計算するのはそもそもおかしい。対数正規分布において計算される平均値と標準偏差のもつ意味は,正規分布のそれらが持つ意味とは全く異なる。たとえば,対数正規分布の平均値±1.96*標準偏差の範囲内にあるデータの割合は95%などではない。

また,mean(x)-1.96*sd(x) を「信頼区間下限値」としているが,それは上で書いたように,あるいは,あなたが最後の方で quantile といっていることからも,データの下から 2.5% 目の推定値でしょう。あなたは「母平均の信頼区間」を得たいのかな?ちがいますよね。

>> 試しにlogの世界でやってみるとおおよそ合います。
>> > quantile(log(x), c(0.025, 0.975))
とありますが,

あなたの対数正規乱数 x は,正規乱数の指数をとったもので,私が上で作った正確な例では x2 です。
> quantile(log(x2), c(0.025, 0.975))
2.5% 97.5%
-1.083910 2.063209
はあたりまえですが,以下と同じですね。 log(x2) == x1, x2 ==exp(x1) ですから。
> quantile(x1, c(0.025, 0.975))
2.5% 97.5%
-1.083910 2.063209

解はあるか?

その前に
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/wrong-sd.html
を読んでみてください。

そこにある図で,正規分布での平均値±標準偏差が対数正規分布では exp(平均値)*/exp(標準偏差) に対応すると書いてあります。

しかし,逆はできないのです。

あなたは図の上側の対数正規分布を対象にしているのですが,対数正規分布で計算した標準偏差は,図の下側の正規分布には変換するすべがないのが分かるかと思います。

解はあると思いますが,私は知らない。

http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf
は,正規分布のパラメータμとσと対数正規分布の平均値と分散の関係式であるが,解析的には成り立つが,実際のシミュレーションデータでは正確には成り立たない。
set.seed を外して何回かやってみると分かるが,x1 の 平均値,標準偏差はいつも指定された値(0.5, 0.8)になるが,それを exp したとたん,x2 の平均値,標準偏差は毎回変わる(つまり,解析的な値とは異なる.)。そのため,解析的な関係式にしたがわない。

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
>
> x2 <- exp(x1)
> mean(x2)
[1] 2.274771 ★これと
> sd(x2)
[1] 2.183789 ●これと

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
>
> x2 <- exp(x1)
> mean(x2)
[1] 2.266348 ★これ
> sd(x2)
[1] 2.08723 ●これ

Re: 正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法
  投稿者:Saito 2018/11/30(Fri) 11:32 No. 22637
コメントありがとうございます。
色々と間違っていたようで申し訳ありません。
質問のコメントについて考えているうちに、私の質問自体がよくなかったことに気づきましたので、最後のほうで新しく例を作りました。二度手間にしてしまって申し訳ありません。私の質問としては、データではなく信頼区間について聞きたかった(ことに気づきました)ということです。申し訳ありません。大変恐れ入りますが、ご教授くだされば幸いです。

先に頂いたコメントについてご返信させていただきます。

> x0 <- rnorm(n, u, sig)
> の x0 の平均値と標準偏差は当然ながら 0.5 と 0.8 ではないということ

についてですが、おっしゃる通りです。ただ、中心極限定理で、0.5と0.8の周りに正規分布してくれるような方法でもOKです。

> しかし,平均値と標準偏差を計算するのはそもそもおかしい。対数正規分布において計算される平均値と標準偏差のもつ意味は,正規分布のそれらが持つ意味とは全く異なる。

仰る通りです。ただ、そういう出力(正規分布を仮定した平均値と標準偏差)がすでにされてしまって(紙になってしまって)いるので、生データを抜きに、それをなんとかしようということでした。

>また,mean(x)-1.96*sd(x) を「信頼区間下限値」としているが,それは上で書いたように,あるいは,あなたが最後の方で quantile といっていることからも,データの下から 2.5% 目の推定値でしょう。あなたは「母平均の信頼区間」を得たいのかな?ちがいますよね。

すみません、mean(x)-1.96*sd(x)はご指摘の通り信頼区間ではありませんでした。また、この指摘を受けて、私も聞きたいことが聞けていないことに気づきましたので、下で改めてシミュレーションデータを作らせていただきました。二度手間になってしまい、申し訳ありません。

>その前に
>http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/wrong-sd.html
>を読んでみてください。

上に同じです。読んでいくうちに、私の聞きたいことが聞けない質問内容になっていることに気づきました。本当に申し訳ありません。

>http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf
>は,正規分布のパラメータμとσと対数正規分布の平均値と分散の関係式であるが,解析的には成り立つが,実際のシミュレーションデータでは正確には成り立たない。

すみません、新しい質問を下記でさせていただく前に、ここについてお聞かせください。解析的に成り立たなくとも、中心極限定理で近い値にはなるのかと思っていたのですが、うまく推定できません。

set.seed(1)
n <- 10000
u <- 0.5
sig <- 0.8

B <- 1000
x <- matrix(exp(rnorm(n*B, u, sig)), n, B)

xlo1 <- apply(x, 2, mean)/exp(1.96*sqrt(log(1+(apply(x, 2, sd)/apply(x, 2, mean))^2)))
xup1 <- apply(x, 2, mean)*exp(1.96*sqrt(log(1+(apply(x, 2, sd)/apply(x, 2, mean))^2)))

xlo2 <- apply(x, 2, quantile, 0.025)
xup2 <- apply(x, 2, quantile, 0.975)

hist(xlo1, xlim=c(0.3, 0.55))
hist(xlo2, add=T, col=2)

hist(xup1, xlim=c(7, 13.5))
hist(xup2, add=T, col=2)

xlo1とxlo2、xup1とxup2はほぼ一致するのかと思ったのですが、どこか誤解しているようです。

以下似たような話ですが、データとパラメータで知りたいことが違っていましたので、質問を修正致します。大変申し訳ありません。

set.seed(1)
n <- 10
u <- 1
sig <- 2

#生物の密度データがあるとする(非負)。平均密度とその信頼区間を知りたい。
y <- exp(rnorm(n, u, sig))
hist(y)

#yというデータがあるもとで、正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
res <- summary(lm(y~1))
lo1 <- res$coef[1]-1.96*res$coef[2]
up1 <- res$coef[1]+1.96*res$coef[2]
> lo1;up1
[1] -1.684872
[1] 22.9317

この信頼区間だと、下限がマイナスの値になってしまっていて、都合が悪いです。
ですので、誤差に対数正規分布を仮定して解析をやり直したものが下記になります。

#yというデータがあるもとで、対数正規分布を誤差分布に仮定したときの平均密度(元のスケール)の95%信頼区間
res2 <- summary(lm(log(y)~1))
lo2 <- exp(res2$coef[1]-1.96*res2$coef[2])
up2 <- exp(res2$coef[1]+1.96*res2$coef[2])
> lo2;up2
[1] 1.345521
[1] 9.318763

私がやりたいこととしては、yを使わず、res$coef[1]とres$coef[2]を使って、lo2とup2を出せないか、ということを考えています。
厳密にlo2とup2に一致しなくとも、正規分布を仮定したときの平均値と標準誤差しかないが、それらから対数正規分布を仮定したときの信頼区間を出す方法、を探しています。
乱数種を替えればyが変わるのでlo1とup1も毎回変わりますが、マイナスのlo1が出る仕組み自体がまずいのでに、それを正の値しか出ない方法に変えたいのです。

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


Rによる行列の処理
  投稿者:KUMON 2018/10/24(Wed) 20:13 No. 22626
Rのプログラムに関する質問です。

各成分が0,1 の2値からなるn×m行列があります。
任意の列に対して、i行目が1ならi+1行目も1となるように変換するコードがうまく書けません。ご教授よろしくお願いいたします。

Re: Rによる行列の処理
  投稿者:青木繁伸 2018/10/24(Wed) 20:52 No. 22627
以下のような列があったとき,あなたはどのような結果をお望みですか?

0
1
0
0
1
1
0
0
0
1
1
0
1
0
1
1
0
1
0

Re: Rによる行列の処理
  投稿者:KUMON 2018/10/24(Wed) 21:45 No. 22629
青木先生

ご返信ありがとうございます。

0
1
1
0
1
1
1
0
0
1
1
1
1
1
1
1
1
1
1

という結果です。

Re: Rによる行列の処理
  投稿者:青木繁伸 2018/10/25(Thu) 21:30 No. 22632
以下のようにすることも出来ますが,
a = c(1, 0, 0, 1, 1, 0, 0, 0, 1, 1,  0, 1, 0, 1, 1, 0, 1, 0)
cat(a)
z = (a | c(0, a[-length(a)]))+0
cat(z)


Re: Rによる行列の処理
  投稿者:KUMON 2018/10/30(Tue) 21:07 No. 22633
青木先生

ご回答ありがとうございます。

列ベクトルa の要素を1つずらしたb を使用して(頭は0で最後尾は削除したもの)
ifelse(a+b>0,1,0)と自解していましたが、

先生のコード (a|b)+0 がよくわかりません。解説していただけないでしょうか?


Re: Rによる行列の処理
  投稿者:KUMON 2018/10/30(Tue) 21:25 No. 22634
分かりました!

論理和を使用しているのですね。+0を加えることでTRUE→1、FALSE→0と変換するわけですね。

参考になります。処理スピードが格段に上がりました。ありがとうございました。


相関係数のp値
  投稿者:コロン 2018/10/25(Thu) 17:16 No. 22630
いつもお世話になっております。

毎回稚拙な質問で恐縮なのですが、以下のデータの相関係数行列はもとまるのですが、それぞれのp値を行列表として求めてくれる関数はないのでしょうか。頑張って探したのですが...。

申し訳ございませんが、よろしくお願いいたします。

q1 q2 q3 q4 q5
3 4 2 5 5
1 2 2 4 5
5 4 3 4 1
4 3 5 2 2
1 2 2 1 3
3 3 3 3 3
4 5 5 4 5

Re: 相関係数のp値
  投稿者:コロン 2018/10/25(Thu) 17:57 No. 22631
青木先生

申し訳ございません。自己解決いたしました。先生のサイトの「相関係数行列の計算と検定」を発見いたしました。


Eviewsにおけるイベントスタディ分析について
  投稿者:ゆう 2018/10/23(Tue) 16:24 No. 22625
・題名の通り、Eviewsを使ってどのような手法でイベントスタディ分析を行うのか手順を教えて頂きたいです。

補足
・回帰分析までは出来ます。そこからAR.CARの出し方がわかりません。
どなたかご存知の方よろしくお願いいたします。苦しんでいます。。。

(エクセルのやり方でも大丈夫です。)


複合アウトカム
  投稿者:あおいひげ 2018/10/17(Wed) 17:43 No. 22623
青木先生

ご回答有り難うございます。
正準相関分析・・・参考になりました。
Rを使用するのはかなりハードルが高いと思っています。
再度、研究デザインから検討し直します。

このスピードで、知りたいことに端的にお返事頂いていることに大変感銘を受けました。今後も、お世話になるかと思います。宜しくお願いします。

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/17(Wed) 22:39 No. 22624
今度は,スレッドを維持するように


複合アウトカム
  投稿者:あおいひげ 2018/10/16(Tue) 16:30 No. 22621
青木先生
回答有り難うございます。

複合アウトカムについては、死亡と再入院など、意味合いが異なるイベントを組み合わせることの意義について見直されてきているようです。
 現実的には、列挙したアウトカムの中でも1番重要視しているものをプライマリーアウトカムとし、それ以外はセカンダリーアウトカムとしてオマケぐらいの位置づけで解析をやるしかないのかもしれませんね。
 因みに、正準相関分析は私が使用しているEZRでは処理できそうにないのですが、他のソフトなら可能なものなのでしょうか。
 

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/16(Tue) 17:55 No. 22622
EZR は,R のシュガーコートなので,R の全てに対応しているわけではない
(EZR にどっぷりというのは避けた方がよい)

R には,cancor という関数がある(stats に含まれているので,library とかしないで,すぐ使える)
使用例(example)も見ることができる
> example(cancor)

cancor> ## No test:
cancor> ## signs of results are random
cancor> pop <- LifeCycleSavings[, 2:3]

cancor> oec <- LifeCycleSavings[, -(2:3)]

cancor> cancor(pop, oec)
$cor
[1] 0.8247966 0.3652762

$xcoef
[,1] [,2]
pop15 -0.009110856 -0.03622206
pop75 0.048647514 -0.26031158

$ycoef
[,1] [,2] [,3]
sr 0.0084710221 3.337936e-02 -5.157130e-03
dpi 0.0001307398 -7.588232e-05 4.543705e-06

以下略


複合アウトカム
  投稿者:あおいひげ 2018/10/15(Mon) 18:15 No. 22619
初めて投稿致します。

医学系臨床研究のデザインを検討しています。統計初心者で、いろいろ試行錯誤で苦しんでおりますのでご教授お願い致します。統計ソフトはEZRを使用しております。

集団をある値が高い群(H群)と低い群(L群)の2群に分けました。100〜200例 VS 300〜500例 (集計を開始していないので見積もり概算です) 。
 アウトカムとして、連続変数と二値変数の2種類の変数を全部で10個くらい想定しています。交絡因子は年齢、性別、既往歴、など10種類くらいを検討しています。
 多変量解析になるかと思いますが、複合アウトカムの処理方法が分かりません。アウトカムを一つずつ、重回帰やロジスティックなどを多重性に気を付けながらやるしかないのでしょうか。

曖昧な質問ですみません。
お教え下さいませ。

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/15(Mon) 19:21 No. 22620
結果が複数の変数ということならば,正準相関分析(重回帰分析の一般形)ということになるかな?実際に使った研究例は多くはないかも。


【R】 リスト構造からの読み取り
  投稿者:明石 2018/10/02(Tue) 21:46 No. 22616
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

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

ベイジアンネット bnlearn パッケージの出力結果のリスト構造からグラフ構造を取得したいと思います。
以下、グラフ描画までのコード例をお示しします。

library(bnlearn) # ベイジアンネットのパッケージ

df <- learning.test # 組み込みのデータセット

res <- hc(df)   # 構造学習 hill-climbing法を適用
plot(res, main = "hill-climbing法による構造学習") # 添付の画像ファイル

str(res)

取得したリスト構造に注釈 銑┐鯢嬪燭靴燭發里髻∨尾にお示しをします。
注釈 銑┐硫媾蠅ら、グラフの親子構造(親は→の元、子は→の先)を出力する
Rプログラムのご教示をいただければ大変に助かります。

誠にお恥ずかしいのですが、まったく手も足も出ません。

所望する出力です。
A,B
A,D
B,E
C,D
F,E

--------------------------------
以下、str(res) の出力結果に、 銑┐涼躰瓩鯢嬪燭靴泙靴拭

List of 3
$ learning:List of 8
..$ whitelist: NULL
..$ blacklist: NULL
..$ test : chr "bic"
..$ ntests : num 40
..$ algo : chr "hc"
..$ args :List of 1
.. ..$ k: num 4.26
..$ optimized: logi TRUE
..$ illegal : NULL
$ nodes :List of 6
..$ A:List of 4             ⇒
.. ..$ mb : chr [1:3] "B" "C" "D"
.. ..$ nbr : chr [1:2] "B" "D"
.. ..$ parents : chr(0)
.. ..$ children: chr [1:2] "B" "D"   ⇒
..$ B:List of 4             ⇒
.. ..$ mb : chr [1:3] "A" "E" "F"
.. ..$ nbr : chr [1:2] "A" "E"
.. ..$ parents : chr "A"
.. ..$ children: chr "E"        ⇒
..$ C:List of 4              ⇒ァ      
.. ..$ mb : chr [1:2] "A" "D"
.. ..$ nbr : chr "D"
.. ..$ parents : chr(0)
.. ..$ children: chr "D"        ⇒
..$ D:List of 4
.. ..$ mb : chr [1:2] "A" "C"
.. ..$ nbr : chr [1:2] "A" "C"
.. ..$ parents : chr [1:2] "A" "C"
.. ..$ children: chr(0)
..$ E:List of 4
.. ..$ mb : chr [1:2] "B" "F"
.. ..$ nbr : chr [1:2] "B" "F"
.. ..$ parents : chr [1:2] "B" "F"
.. ..$ children: chr(0)
..$ F:List of 4              ⇒
.. ..$ mb : chr [1:2] "B" "E"
.. ..$ nbr : chr "E"
.. ..$ parents : chr(0)
.. ..$ children: chr "E"        ⇒
$ arcs : chr [1:5, 1:2] "A" "A" "C" "B" ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "from" "to"
- attr(*, "class")= chr "bn"

何卒どうぞ、よろしくお願いいたします。
いつもありがとうございます。


Re: 【R】 リスト構造からの読み取り
  投稿者:青木繁伸 2018/10/02(Tue) 23:14 No. 22617
この件については,何も知りませんが,
以下のようにすれば,簡単に解が得られのでは?
> res$arcs
from to
[1,] "A" "B"
[2,] "A" "D"
[3,] "C" "D"
[4,] "B" "E"
[5,] "F" "E"
See bn-class for details.

arcs: the arcs of the Bayesian network (a two-column matrix, whose columns are labeled from and to). Undirected arcs are stored as two directed arcs with
opposite directions between the corresponding incident nodes.

Re: 【R】 リスト構造からの読み取り
  投稿者:明石 2018/10/03(Wed) 19:59 No. 22618
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

今回も大変に良い勉強をさせていただきました。
心から御礼を申し上げます。


ゼロでの割り算
  投稿者:青木繁伸 2018/09/09(Sun) 11:10 No. 22612
中澤さんが,2x2 より大きい行列で R の mcnemar.test の結果が NaN になる場合があることについて,書いている。

http://minato.sip21c.org/bulbul2/20180908.html
http://minato.sip21c.org/bulbul2/20160913.html
http://minato.sip21c.org/bulbul2/20160911.html

「対角成分の和がゼロの組合せを計算からスキップできる拡張マクネマー検定の関数」として
mct <- function(x) {
L <- NROW(x)
X <- 0
N <- 0
for (i in 1:(L - 1)) {
for (j in (i + 1):L) {
s <- (x[i, j] + x[j, i])
if (s > 0) {
N <- N + 1
X <- X + (x[i, j] - x[j, i]) ^ 2 / s
}
}
}
return(c( X2 = X, df = N, p = (1 - pchisq(X, N))))
}
を挙げている。
しかし,「このやり方が正しいという理論的根拠を見つけられないので悩ましいところだが」とあるのだが。

「N23とN32がゼロ,という状態だと,(N23+N32)が分母になる項で"Division by zero"エラーがでてしまって,カイ二乗値の計算結果自体がNaNになってしまう」そこで「分母がゼロになる組合せをスキップする」ということを書いている。

この記述は正しくない。N23=N32=0のとき,つまり(0-0)/(0+0) = 0/0 の計算結果は,NaN である("Division by zero"ではない。Python だと,ZeroDivisionError: division by zero になる)。何かに NaN を足すと NaN になり,それ以降 NaN に何を足しても NaN のまま。ということで,R の mcnemar.test は NaN を返すと,中澤さんも書いている。

N23=N32=0 のときの計算結果が NaN になった場合は足し算をスキップすればよいわけで,N23=N32=0 の場合は計算をスキップしても同じではある。またその場合には自由度のカウントはしない。

以下のような関数を書いてみた。
mct2 <- function(x) {
i <- col(x)[lower.tri(x)]
j <- row(x)[lower.tri(x)]
a <- mapply(function(i, j) (x[i, j] - x[j, i]) ^ 2 / (x[i, j] + x[j, i]), i, j)
chisq <- sum(a, na.rm = TRUE) # NaN も除かれる
df <- sum(!is.nan(a)) # !is.na(a) でもよい
return(c(X2 = chisq, df = df, p.value = pchisq(chisq, df, lower.tail = FALSE)))
}
「ゼロ割」は避けるのではなく,積極的に正しく対処するべきということか?(同じようなことではあるが)

蛇足ではあるが,自由度を減ずるのは,カイ二乗分布の再生性から明らか(自由度 1 のカイ二乗分布にしたがう確率変数 k 個の和は自由度 k のカイ二乗分布にしたがう)。

R の mcnemar.test を修正するとすれば,
STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)])

STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
PARAMETER <- sum(!is.nan(y[upper.tri(x)]^2/x[upper.tri(x)]))
にする(前の方の PARAMETER <- r * (r - 1)/2 は不要)。
プログラムにおいて,特殊な条件の発生を見逃すと命取りというところ。
しかしさすがに,sum(y[upper.tri(x)]^2/x[upper.tri(x)]) とは,にくいねえ。
> x = matrix(c(c(7, 2, 5, 3, 11, 2, 0, 3, 4, 0, 5, 9, 3, 15, 14, 10)), 4, 4)
> mct(x)
X2 df p
15.428836864 5.000000000 0.008678852
> mct2(x)
X2 df p.value
15.428836864 5.000000000 0.008678852

> mcnemar.test(x) # 元のまま

McNemar's Chi-squared test

data: x
McNemar's chi-squared = NaN, df = 6, p-value = NA

> mcnemar.test2(x) # 修正したもの

McNemar's Chi-squared test

data: x
McNemar's chi-squared = 15.429, df = 5, p-value = 0.008679
なお,「"Division by zero"を避けるために全部のセルに0.01を足す」というのは明らかに間違い(自由度が減らない)。

Re: ゼロでの割り算
  投稿者:青木繁伸 2018/09/09(Sun) 12:49 No. 22613
2 0
0 5

というのと

2 1
1 5

というのは,0-0 = 1-1 = 0 という情報を与えるという意味では同じなので,自由度は減らさなくてもよいか

STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
だけの変更で,これは実質的には NaN の場合に 0 を加えるというのと同じであるから
PARAMETER はそのまま PARAMETER <- r * (r - 1)/2

とすれば,mcnemar.test(x+0.01) もありか。ただし,加える数は他のまともなセルのことも考えるとできるだけ小さい方が良いだろうとは思う。例えば +1e-15 とか

Re: ゼロでの割り算
  投稿者:中澤 2018/09/09(Sun) 21:30 No. 22614
記事ありがとうございます。
確かにNaNが正しいですね。不正確な記述でした。
自由度の扱いは悩ましいですね。スキップでなくて分子がゼロなのでゼロとみなすということですね。勉強になりました。

Re: ゼロでの割り算
  投稿者:中澤 2018/09/19(Wed) 17:55 No. 22615
http://minato.sip21c.org/bulbul2/20180915.html
に書きましたが,カテゴリ数が3以上の場合の拡張マクネマー検定と同様の帰無仮説を検定する方法としてBhapkarの検定というのがあって,irrパッケージのbhapkar()で実行できました。いま書いている論文ではこれを使おうと思います。


パネルデータの欠損値補完について
  投稿者:羽島 2018/09/04(Tue) 09:21 No. 22609
青木先生

羽島と申します、この度はお世話になります。

Rを使って以下のようなデータを分析しており、その中で住所コードの欠損値を補完する必要が出てきました。LOCFのような手法での補完では精度に不安があるため、個人IDを参照することで欠損値を補完したいのですが、どのようにプログラムを書けば良いでしょうか?

お忙しいところ大変恐縮ですが、ご教示頂けましたら幸いです。

個人ID  年   住所コード

100   2015    810
100   2016    NA
100   2017    810
201   2015    NA
201   2016    NA
201   2017    666

Re: パネルデータの欠損値補完について
  投稿者:青木繁伸 2018/09/04(Tue) 13:03 No. 22610
あまり簡単なプログラムは書けないかな
d <- data.frame(ID = c(100, 100, 100, 201, 201, 201, 330, 330, 330),
year = c(2015, 2016, 1017, 2015, 2016, 2017, 2015, 2016, 2017),
code = c(810, NA, 810, NA, NA, 666, NA, 543, NA))
tbl <- unique(na.omit(d[c("ID", "code")]))
for (i in seq_len(nrow(d))) {
if (is.na(d$code[i])) {
d$code[i] = tbl[which(d$ID[i] == tbl[,1]), 2]
}
}
d

Re: パネルデータの欠損値補完について
  投稿者:羽島 2018/09/04(Tue) 18:15 No. 22611
青木先生

お返事の方ありがとうございます。

元のデータは3万件ほどの個人データがあり手作業での修正を覚悟しておりましたが、
先生のお力で無事補完することが出来ました。

自身の勉強不足を恥じるばかりですが、この度は本当にありがとうございました。

[1] [2] [3] [4] [5] [6] [7]

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