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

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


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

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


消去された質問
  投稿者:青木繁伸 2016/06/28(Tue) 19:12 No. 22052
用が済んだらあとは証拠隠滅か...なんだかな

質疑応答の概要

質問:句読点の使用頻度など,文章の特性が文章作成指導法によって変化するかどうか,被検者にたいして,指導前と指導後に14000字程度の文章を書かせ,特徴をデータ化。例えば,
    読点使用有り 読点使用なし
指導前 回数     回数
指導後 回数     回数
のようになっている。どのような検定を使えば良いか。

回答1:フィッシャーの正確確率検定を使えば良いのではないか。

再質問:フィッシャーの正確確率検定は,データが少ないときと書いてあるが,多くても使って良いのか。
査読者には,マクネマー検定を使うべきではないかといわれているが,そうなのか?

回答2:データが多いときでも使える(どんなときにも正確なんだから)。マクネマー検定は対応のあるデータの場合なので,今のデータは対応がないんだから,マクネマー検定なんかつかえない。
ところで,被検者は一人なのか?

それに対する回答:一人だ。

回答3:被検者一人じゃ統計学で扱うには不適切。サンプルサイズを決めることもできない。検定ではなく,何が何パーセント増えた(減った)で十分なんじゃないか?

それに対する回答:一人だけど,何らかの方向性をいうことはできるだろう。
でも,回答ありがとう。これを参考に査読者に回答,再投稿する。

そして,きょうの夕方,発端の投稿が消去されたので,それに端を発する質疑応答が,全部消えた。

このデータにマクネマー検定を提案するようなヘボ査読者は,この実験デザインの問題点には全く気づかないということなんだろうか。なんという学会だろうね?


偏F値によるステップワイズ変数選択の参考資料をご教示願います。
  投稿者:米谷 2016/06/23(Thu) 17:12 No. 22040
お世話になっております。
以下のウェブサイトにある偏F値によるステップワイズ変数選択について質問があり、投稿をさせていただきました。
http://aoki2.si.gunma-u.ac.jp/R/sreg.html

こちらにある関数のうち、fmax、fmin、step.outの各関数内で使われている理論式を説明がある参考資料がございましたらご教示いただきたく存じます。

特にfmax、fminのsweep処理や、step.out関数のF値や偏回帰係数の標準偏差等について、参考資料がございましたら、そちらから内容を理解したく、ご教示いただければ幸いです。どうぞよろしくお願いいたします。

米谷

Re: 偏F値によるステップワイズ変数選択の参考資料をご教示願います。
  投稿者:青木繁伸 2016/06/23(Thu) 18:13 No. 22042
古い本になりますが,奥野忠一他「正・続 多変量解析法」日科技連,に書かれていたと思います。

Re: 偏F値によるステップワイズ変数選択の参考資料をご教示願います。
  投稿者:米谷 2016/06/23(Thu) 19:19 No. 22044
青木先生

お忙しいところご教示いただき誠にありがとうござました。

偏F値を使った変数選択の理論的な流れを紹介する資料がなかなか見つからず困っておりました。その本で確認させていただきます。

米谷


Shirley-Williams検定
  投稿者:赤岳 2016/06/23(Thu) 16:38 No. 22039
いつもご指導ありがとうございます。
さて、Shirley-Williams検定をRで実行するスクリプトを青木先生のサイトで見つけました。
source("http://aoki2.si.gunma-u.ac.jp/R/src/Shirley-Williams.R", encoding="euc-jp")
出力結果をみると
> Shirley.Williams(data, group, method="up")
4 3 2
2.189899 1.091791 1.243468
となっていますが、ここに表示されている2.189899等は対照群を1としたときの大きさだと思われますが、統計量やP値は、出力されないのでしょうか。
どうかP値の出し方を教えてください。よろしくお願いします。

Re: Shirley-Williams検定
  投稿者:青木繁伸 2016/06/23(Thu) 18:11 No. 22041
ページに書いてあるように
永田靖・吉田道弘(1997)「統計的多重比較法の基礎」サイエンティスト社
を参考にプログラムを書いたのですが,P 値の求め方は,書いてなかったと思います。
付表を参照だったかな?
出力される数値が何を表すかを含め,前述の本を見てください。

Re: Shirley-Williams検定
  投稿者:赤岳 2016/06/23(Thu) 18:33 No. 22043
青木先生
早速お返事いただきありがとうございました。
参考本を確認してみます。


非正規分布の検定について
  投稿者:下山 2016/06/18(Sat) 18:43 No. 22030
度々お世話になります。

ある要素Aの有無に対してスコアB(%)、スコアC(%)というN=150程度の表が有ります。
A B C
有 10% 17%
無 15% 18%
有 20% 28%
・・・・

A-B、A-Cについてそれぞれ要素Aの有無とスコアの間に差が有るか、有るとすればどちらが大きいかを検定したいと考えています。B、Cの分布は正規性検定ではどちらも非正規という結果でしたので、U検定をしようかと思いました。U検定ではBもCもP<0.01でした。分散分析を行ったところA-Bは有無の2群間で分散は等しくなく、A-Cは等分散という結果になりました。

A-Cは中央値で評価しようかと思いましたが、A-Bについては有群と無群どちらがスコアが高いかをどう評価しようか、と壁に当たりました。
どのような方法がありますでしょうか?中央値を用いても良いのでしょうか?

もし正規性の検定以降の手順についても、間違っているところがありましたらご指摘いただけましたら大変ありがたく存じます。どうぞよろしくお願いいたします。

Re: 非正規分布の検定について
  投稿者:青木繁伸 2016/06/18(Sat) 20:32 No. 22031
そもそも,

> A B C
> 有 10% 17%
> 無 15% 18%
> 有 20% 28%
> ・・・・

というデータの見方がわかりません。

それと,検定に必要なのは,% ではなく,実数です。
本当のデータ(実数)が明らかになって困るとは思えないので,それぞれが何を表すかは言う必要はないのだから,ちゃんと説明する方がよいと思いますけど?
じゃないと,返答のしようがない。

Re: 非正規分布の検定について
  投稿者:下山 2016/06/19(Sun) 08:38 No. 22032
青木先生、失礼いたしました。
BとCはそれぞれ異なる計算式から算出した予測死亡率で、どちらもいくつかの調査項目から計算しています。ですので%単位になっています。

取り急ぎお返事させていただきました。

Re: 非正規分布の検定について
  投稿者:青木繁伸 2016/06/19(Sun) 17:49 No. 22033
検定できないでしょう

Re: 非正規分布の検定について
  投稿者:下山 2016/06/19(Sun) 20:39 No. 22034
青木先生 ご回答ありがとうございます。
このデータに関しては検定できない、ということで分析から外す方向で考えます。

後学のために以下のような場合にどう評価していくかも教えていただけましたらありがたく存じます。

N=150
A=某塾に通っているかどうか
B=期末試験の点数(通塾の有無で分けた場合等分散ではない)
C=統一模試の点数(通塾の有無で分けた場合等分散である)

として、点数がどちらも正規分布していない場合、通塾の有無で期末試験の点数に差があるか、あるとすればどちらが点数が高い傾向があるのか、および通塾の有無で統一模試の点数に差があるか、あるとすればどちらが点数が高い傾向があるのかを評価したい。この場合であればデータは実数となるかと思いますが、ご教授いただけましたらありがたく存じます。
大変恐縮でございますが、どうぞよろしくお願いいたします。

Re: 非正規分布の検定について
  投稿者:青木繁伸 2016/06/19(Sun) 21:04 No. 22035
この例と,前のと,どのような関係があるのですか?
本当は関係がないのに,別の問題だとして回答したことを(つまり,新たに提示されたデータについてはこれこれの検定が適切でしょうということを),あなたが勝手に(状況は全く違う)この問題にも適用するかもしれないという懸念があるのに,回答できないでしょう?

抽象的なことばかり言っていないで,データの性質も含めて記述しないと,なんの回答もできません(私のせいにされるのは心外で,そんなことこわくて,できません)。

Re: 非正規分布の検定について
  投稿者:下山 2016/06/20(Mon) 14:06 No. 22036
青木先生

検定を根本的に理解しておらず多々考え違いをしていることを痛感いたしました。ご気分を害してしまいましたことお詫びいたします。

理解が乏しい中で自分なりに解釈すると、正規分布が仮定できる場合にしろ分割表にしろ非正規分布の場合にしろ、実際のデータを吟味しないといけなくて、初心者向けの統計本に記載されているような、パラメトリックなら対応の有無を考慮してt検定やANOVAで、分割表ならχ二乗で、ノンパラメトリックなら・・・で、というように馬鹿の一つ覚えで1対1に決めてできるものではないということかと考えました。

必要十分なデータの性質についての情報を記載できる自信はありませんが(データの性質といわれても自分には分布や分散、外れ値くらいしか思いつかない程度の理解で)、添付の画像のようなデータです。Logitはそれぞれ異なる論文に記載されていた回帰式より計算しております。

はじめに先生がご回答されたように、こういうデータは検定できないのでしょうか(あるいは検定を行う意味・意義がない、検定すべきではない)。


Re: 非正規分布の検定について
  投稿者:鈴木康弘 2016/06/22(Wed) 07:15 No. 22037
 せっかく具体的なデータを出していただいたのに、一般論ですが、
 正規分布し分散も等しい2変数を利用して、2群に差があるかを検定するなら、多変量分散分析を使うはずです。

 でもこの条件に合わなければこれは使えない。
 だから2変数を一度に使うことはあきらめざるをえない(と思う)。

 正規分布もしてないし分散も異なる2群を単変量で検定するのだったら、Brunner-Munzel検定というのがあるようです。

Re: 非正規分布の検定について
  投稿者:下山 2016/06/22(Wed) 20:12 No. 22038
鈴木先生

御助言ありがとうございます。Brunner-Munzel検定は初めて聞きました。
調べてみたいと思います。


手法について
  投稿者:コロン 2016/06/12(Sun) 12:02 No. 22026
いつもお世話になっております。

早速ではございますが、統計手法についてご教授ください。

やりたいことは、○○の性格の学習者は△△の単語学習ストラテジーを取るようだ、を見たいと思います。

以下は案です。

性格については心理学の質問紙法を用いて被験者集団を分類することとします。また単語の学習ストラテジーについても単語学習のストラテジー質問項目を作成して、因子分析を行って分類します。

ただこれらをどのように関連付けて分析すべきなのかよくわかりません。因果の分析なのか、単純に、例えば、3×4の表を作成して、偏りの検定なのか?。

ご教授お願い致します。

Re: 手法について
  投稿者:青木繁伸 2016/06/12(Sun) 22:01 No. 22027
まずは,「○○の性格の学習者は△△の単語学習ストラテジーを取るようだ」ということについてある程度の客観性のある事実を確認する必要があるというのは置いておくとして。
「単語学習のストラテジー質問項目を作成」は重要なポイントでしょう。
「性格については心理学の質問紙法を用いて被験者集団を分類」というのも,適切に被検者集団が分類できるのかというのも。
これらについて,適切な測定ができれば,そのデータの分析方法というのは決まるでしょう。適切な測定ではないデータであれば,どんな分析をしても異論が出ることでしょう。

Re: 手法について
  投稿者:コロン 2016/06/13(Mon) 09:58 No. 22028
青木先生

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

性格については,あくまで「例」ではありますが,4つくらいに分けられるようなものを使う予定です。これによって,「この学習者はタイプA,この学習者はタイプB」というように,カテゴリカルな値を被験者に与えます。

語彙学習ストラテジーですが,語彙学習ストラテジーに関するアンケート項目は,これまでにもたくさんされており,それらを効率的に使用し,今回の調査のための質問項目を作成いたします。それを因子分析して,例えば,4因子構造だとなった場合,それぞれの因子得点を利用したいと考えています。

「客観性ある事実を確認する必要がある」というのはその通りだと思います。また性格も4通りだけに分けることが可能なのかといった異論は出てくると思います。あくまで例えばということで,ご理解いただき,上述のデザインへのコメント,そして分析手法についてアドバイスをいただければと思っております。

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

Re: 手法について
  投稿者:青木繁伸 2016/06/15(Wed) 20:34 No. 22029
多変量解析手法をすべてを統合した,共分散構造分析(構造方程式モデリング、SEM)とか?


2X3分割表のFisher検定の結果の記述について
  投稿者:下山 2016/06/04(Sat) 12:51 No. 22018
青木先生、お忙しいところ恐れ入ります。

ある患者層で特定の薬の処方の有無と性別や年齢層などとの関連を調べるような調査を行いました。
性別は2x2の分割表、年齢層は2x3の分割表を作成しそれぞれについてFisher検定を行いました。2x2表の結果についてはP値、Odds比、95%信頼区間を記載し、2x3表についてはOdds比や信頼区間はそもそも計算できないのではないかと考えP値のみ記載し、OR、CIはnot availableとした表を結果として提出したところ、

・単変量解析はフィッシャーの正確確率検定を行ったと記載がありますが、通常この検定は2×2表に適応されると思われます。年齢層の項目は3つに分かれておりますがこれもフィッシャーの正確確率検定が行われたのでしょうか?
・年齢層のORがnot availableであることや、p valueが一つしか無い点が、なぜか解りませんでした。variableのどれかがreference groupとなり、そのreferenceに対するodds ratio、その95%CI,p valueがそれぞれ出るはずです。間違った分析をされているか、方法論の説明が不十分かのどちらかと思います。

という査読コメントをいただきました。2x3表でp=0.4でしたので各群間での比較は行わずn/aと記載しておいたのですが、このような場合には最初から各群間比較を行うものなのでしょうか?それとも指摘されたとおり間違った分析を行ってしまったのでしょうか。
(なお、有意差が有りそうな場合は、有意水準の調整をした上で個々の群間で総当りのFisher検定を追加で行うことを考えていました)

初歩的な内容の質問で大変恐縮でありますが、何卒よろしくお願いいたします。

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:青木繁伸 2016/06/05(Sun) 07:59 No. 22019
> 単変量解析はフィッシャーの正確確率検定を行ったと記載がありますが、通常この検定は2×2表に適応されると思われます。

不勉強な査読者もいるものだと,驚きますね。どこの学会ですか。

> 年齢層のORがnot availableであることや、p valueが一つしか無い点が、なぜか解りませんでした。

「2x3表でp=0.4でしたので各群間での比較は行わずn/aと記載しておいた」ということをはっきり書いてやらないと,この馬鹿査読者には通じないのでしょう。

しかたないので,よく説明してやってください。

> 有意差が有りそうな場合は、有意水準の調整をした上で個々の群間で総当りのFisher検定を追加で行う

この場合には,検定の多重性に注意してくださいね。

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:下山 2016/06/05(Sun) 08:17 No. 22020
青木先生

ご回答くださりありがとうございます。分析に自信を失いかけて凹んでおりましたが、おかげさまで本文への追加記載とコメントへの返信を自信を持ってできます。

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:鈴木康弘 2016/06/06(Mon) 06:53 No. 22021
 年齢層で2X3分割表を作ったということは、分割表に順序があるんですよね?
 そういう場合にFisherの検定をしていいのかどうか?

 一人一人の年齢データが入手できるなら、t検定をしておしまい。(オッズ比がほしければロジスティック回帰。)
 それが入手できず、分割表データにせざるをえないんだったら、カテゴリーごとに1点2点3点と点数をつけて、マン・ホイットニーをするのがオーソドックスな方法では?

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:下山 2016/06/08(Wed) 13:06 No. 22022
鈴木先生

ご指摘ありがとうございます。確かに年齢は順序尺度という扱いでした。年齢についての生データは生後数日、数か月、○歳と単位がバラバラなのでカテゴライズするのが妥当と考えております。
もし年齢ではなくて主要疾患(呼吸器、循環器、消化器)というような順序がないようなものならFisherの検定でも構わない、という理解でよろしいでしょうか。

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:鈴木康弘 2016/06/09(Thu) 10:21 No. 22023
その通りと思います。

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:scdent 2016/06/10(Fri) 09:43 No. 22024
年齢層の2x3分割表についてですが、年齢層が順序尺度だからといって安易にマン・ホイットニーにするというのは問題があります。

中心位置の相違を検出しようとするならそれでよいのですが、分布の違いを見る場合はχ2(あるいはFisher)検定の方が適していたりします。

例をあげますと、薬剤群では低年齢層と高年齢層で副作用発現率が高く、プラセボ群では中年齢層で副作用発現率が高い場合、マン・ホイットニーではその違いが検出できないのです。

中年齢層では活性が高いので副作用発現率が高いのですが、薬剤の影響で低年齢や高年齢層に副作用が多発するという可能性は十分あります。

安全性を扱うものとして、こういう特性を見逃すことがないようにしなくてはなりません。

Re: 2X3分割表のFisher検定の結果の記述について
  投稿者:鈴木康弘 2016/06/10(Fri) 13:52 No. 22025
なるほど。


R 文字列の生成
  投稿者:明石 2016/05/27(Fri) 22:10 No. 22015
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

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

------
添付の画像をご覧ください。

お示しをしました例題は、説明のためにサイズを小さくしてあります。

実際のデータは大規模であり、0要素が多くあるスパースなデータです。

入力データ("in.csv")を読み込み、以下のルールで、
テキストファイル("out.txt")を生成出力するRプログラムを作成しています。

【ルール】

・表のセルの数値(非負の整数)を読み取り、
 数値が1以上の場合には、
 セルの列名("aaa", "bbb", "ccc", "ddd", "eee")を、
 セルの数値の回数分を繰り返した文字列に置き換える。

・表のセルの数値(非負の整数)を読み取り、
 0の場合には、空文字に置き換える

・文字列のセパレータは、半角スペースとする。

以下のRプログラムを作成しましたが、2つの問題が残ります。
私の力では、解決できません。

良いお手本をお示しいただけると、大変に助かります。
どうぞ、よろしくお願いいたします。

【問題点】

・実際に利用するデータは大規模であることから、
 多重ループの処理を回避して、性能向上を図りたい。

・実際に利用するデータは0要素が多いことから、
 生成出力したテキストファイルには、区切り記号の半角スペースが多く混入します。

 この冗長な半角スペースを、gsub()関数で置換して削除しようと試みましたが、
 私の知識では、できませんでした。

以下、プログラムです。

dat <- read.table("in.csv",header=T, sep=",")

nr <- nrow(dat)
nc <- ncol(dat)

# セルの数字を文字列に置き換える

for(i in 1:nr) {
for(j in 1:nc) {

r <- dat[i,j]

if (r >= 1) {
v <- rep(colnames(dat)[j], r)
v <- paste(v , collapse=" ")
dat[i,j] <- v
}

if (r == 0) {
dat[i,j] <- ""
}
}
}

gsub(" ", "", dat)

write.table(dat, "out.txt", col.names = F, row.names = F, quote = F, sep = " " )


Re: R 文字列の生成
  投稿者:青木繁伸 2016/05/27(Fri) 23:12 No. 22016
テストデータを「生成する」プログラムも付けてくれるといいのだけど。
# テストデータ生成
dat = matrix(0, 300, 52)
colnames(dat) = c(LETTERS, letters)
set.seed(123)
idx = cbind(sample(300, 300*26, replace=TRUE),
sample(52, 300*26, replace=TRUE))
dat[idx] = sample(6, 300*26, replace=TRUE)
dat = data.frame(dat)
# 処理プログラム
names = colnames(dat)
fn = "out.txt"
if (file.exists(fn)) invisible(file.remove(fn))
invisible(apply(dat, 1, function(x) cat(rep(names, x), "\n", file=fn, append=TRUE)))
20倍くらいは速い。ただし,行末に1つ空白が付く。これを回避するには若干の追加処理が必要。

【御礼】 Re: R 文字列の生成
  投稿者:明石 2016/05/28(Sat) 06:52 No. 22017
青木先生、
お忙しいところを失礼いたします、明石と申します。

ご丁寧なご教示をいただき、誠にありがとうございました。
また、助けていただきました。

投稿に際して、テストデータの件、不手際がございましたこと、お詫びいたします。

今後、投稿させていただく際には、お手を煩わせることが少なくなるように注意をいたします。

ご教示くださいました技は、あまりにも研ぎ澄まされていて、にわかには理解できませんので、
この土日、楽しみながら、じっくりと勉強させていただきます。

また、先生がお書きになられたように、確かに、体感的には数十倍の速度があります。

大変に良いお手本を示してくださり、心から御礼を申し上げます。
ありがとうございました。


母平均の区間推定における標準正規分布表の利用について
  投稿者:hase 2016/05/17(Tue) 12:57 No. 22010
青木先生・皆様

統計学を勉強し直しているhaseと申します.

母平均の区間推定における標準正規分布表の利用についてご教示下さい.
なお,手元の教科書付属の正規分布表にはPr(0 ≦ z)の値が小数第4位まで載っております.

99%信頼区間(=有意水準0.01)のとき,0.99 = 0.4950 * 2ということで,表中の0.4950に対応するzを探したところ,0.4950丁度はなく,0.4949(z=2.57)と0.4951(z=2.58)の2つが最も近い値として見つかりました.

Rでqnorm(0.005, lower.tail=FALSE)とすると,2.575829となることから,より値の近い0.4951(z=2.58)を使うべきでしょうか?

それとも,
0.4951(z=2.58)では区間推定した母平均が過大な推定値となることから,0.4949(z=2.57)を使うべきでしょうか?

この例に限らず,信頼区間から出てくる(0.4950に相当する)値に最も近い表中の値を用いるべきなのか,それとも,その値以下の最も大きい値を用いるべきなのか,統計学としてはどちらの考え方が正しいのでしょうか?

※以下のURLでの議論のように,精度の良い正規分布表を使う,あるいはRなどコンピュータを使って計算することが出来ない場合での疑問です.
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc033/06569.html

以上,よろしくお願いいたします.

Re: 母平均の区間推定における標準正規分布表の利用について
  投稿者:太郎 2016/05/17(Tue) 16:41 No. 22011
「正規分布表」、「補間」でググるとこの答えに近いものがでてきました。

Re: 母平均の区間推定における標準正規分布表の利用について
  投稿者:hase 2016/05/17(Tue) 18:49 No. 22012
太郎さん,ありがとうございます.

滋賀大学の熊澤先生の資料で線形補間法についてはよく分かりました(昔習ったような気もします).
http://www.biwako.shiga-u.ac.jp/sensei/kumazawa/toukeigaku/06/0623k.pdf
線形補間法の結果は0.2575となるので,qnorm(0.005, lower.tail=FALSE)から求めた2.575829と比較しても,過大ではありません.
ただ,線形補間には過大推計の可能性があるので,少し怖い気もします.

一方,
http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1233951557
のscairflaimさんの回答にある,
> 表にはP(X<x)≦0.1になる値がないときは、0.1より小さい方の(以下なので、大きくなってはいけない)0.1に近い値を選んでもほぼ問題ないです。なぜなら0.1以下なら条件を満たすし、0.1との差はすごく小さいので、確率に出てくる影響は小さいのでその差は無視できるからです。
> そもそも表を使いなさいと言われたら、そこまでの精度は求められていません。
という考え方なら,少なくとも過大推計とはなりません.

どちらの(あるいは別の)考え方が適切なのでしょうか?
引き続きご教示いただければ幸いです.

Re: 母平均の区間推定における標準正規分布表の利用について
  投稿者:hase 2016/05/18(Wed) 10:05 No. 22013
> どちらの(あるいは別の)考え方が適切なのでしょうか?
この記述が分かりづらいので補足いたします.

1.線形補間法が一般に使われることは分かりましたが,過大推計の恐れがあるように思います.これは気にしないで良いものでしょうか?

2.scairflaimさんの回答の考え方は,適切なものでしょうか?

3.もし,1.と2.どちらも選択できる場合は,どちらを取るのが適切なのでしょうか?

4.あるいは,他により適切な方法があるのでしょうか?

Re: 母平均の区間推定における標準正規分布表の利用について
  投稿者:青木繁伸 2016/05/19(Thu) 11:00 No. 22014
実際問題として,母分散が既知と(仮定)した時点で推計値の精度は落ちているわけで,Z=2.57, 2.58 のいずれを取るとしても,標準誤差(σ/√n)の 0.39 % ほどしか推計値の違いはないわけで,どっちをとっても過大推計・過少推計の問題は生じないでしょう。
試験やレポートの場合であっても,どちらに基づいて答えを書いても○になると思いますが。


ばらつき補正係数
  投稿者:すーさん 2016/05/14(Sat) 18:16 No. 22003
青木先生

重要顧客の担当者から下記質問をメールで受けて困っていたところ、このサイトにたどり着きました。どうぞよろしくお願いします。

少ないデータ数で求めたばらつきは、全体の集団(母集団:n=∞)のσより小さくなるため、データ数を考慮し、ばらつきを補正する必要があることは理解できるのですが、下記の補正係数(1.45、1.17)はどのようにして求めたものか? 質問を受けています。

正規分布の場合、
  n=10で算出したσの1.45倍が母集団のばらつき、
   n=20で算出したσの1.17倍が母集団のばらつき

 n=30、100・・・などの場合、補正係数はどうやって求めるのか?
何卒、ご教示をよろしくお願いいたします。

Re: ばらつき補正係数
  投稿者:青木繁伸 2016/05/14(Sat) 21:39 No. 22004
2010 年頃の質疑応答の記録があります。当方の掲示板の参照もありますが,結論としては「よくわからない」ということになっていますかね。

http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1436950684

Re: ばらつき補正係数
  投稿者:すーさん 2016/05/15(Sun) 11:08 No. 22007
青木先生

早速のご回答を頂き大変ありがとうございます。

記載頂いたサイトは私も2日前にチェックさせて頂いたのですが良く判りませんでした。
補足:先方からは添付の資料を頂いておりますが、この図の出処が不明です。

 カイ2乗分布から求めた90%信頼区間の上限値?、不偏分散と標本分散から求めたシグマの比?、補正係数として1/c4で計算しても一致しません。
(エクセルで乱数を発生させシミュレーションもして見たのですが、標本標準偏差であっても不偏推定量である(=stdev関数)を用いることで、n=100以上であれば、おおよそ母標準偏差に近い値となっている)

1)結論はやはり「よくわりません」という事になりますでしょうか?
2)サンプル数は100以上であれば補正の必要なしと考えて良いでしょうか?

ご面倒をお掛けしますが宜しくお願いいたします。

Re: ばらつき補正係数
  投稿者:青木繁伸 2016/05/15(Sun) 11:14 No. 22008
「出典は某工業製品メーカー作成の統計教育テキスト」ということなのですが,そちらへ問いあわせる...にしても,「某工業製品メーカー」では無理ですね?

Re: ばらつき補正係数
  投稿者:すーさん 2016/05/15(Sun) 19:01 No. 22009
これ以上は調査不可能といった所ですね。
このような質問にお付き合い頂き大変有難うございました。
また何か出てきましたら質問させて頂くかも知れませんが宜しくお願いします。


R str()関数の結果の取り出し、整形出力
  投稿者:明石 2016/05/13(Fri) 21:51 No. 22000
青木先生 様;

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

データフレームの構造は、str()関数で取得することができます。

以下に簡単な例でお示しをしますように、
str()関数の結果から、列の名前と型を、2列のデータフレームとして出力したいと思います。

str()関数の結果の取り出し方法が分からないため、作成することができません。

何卒、ご教示をよろしくお願いいたします。
失礼いたします。

> str(dat)

'data.frame': 1000 obs. of 9 variables:
$ sid : chr "X100111" "X100115" "X100123" "X100143" ...
$ name : chr "安田博美" "木村克己" "桜井真二" "入江優秀" ...
$ gender : Factor w/ 2 levels "女","男": 1 2 2 2 1 2 1 1 1 2 ...
$ BirthDay : Date, format: "1999-12-09" "1981-07-21" ...
$ BloodType: Factor w/ 4 levels "A","AB","B","O": 3 1 3 4 4 4 2 1 4 3 ...
$ class : Factor w/ 3 levels "1","2","3": 3 2 2 3 2 1 3 1 1 1 ...
$ score1 : num 3.3 4.1 3 1.9 3.6 2.6 4.5 4 2 3.2 ...
$ score2 : num 3.3 4.1 2.5 3 4.4 2.8 3.9 3.8 3.1 4.8 ...
$ score3 : num 2.5 3.8 3.9 3.3 3 3.1 4.4 4.1 2.8 4.8 ...

⇒ 以下のような、データフレームを作成したいと思います。

名前,型
sid,chr
name,chr
gender,Factor
BirthDay,Date
BloodType,Factor
class,Factor
score1,num
score2,num
score3,num

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

Re: R str()関数の結果の取り出し、整形出力
  投稿者:青木繁伸 2016/05/14(Sat) 05:12 No. 22001
a = data.frame(names(iris), sapply(iris, class))
colnames(a) = c("名前", "型")
print(a, row.names=FALSE)
などはどうでしょうか

Re: R str()関数の結果の取り出し、整形出力
  投稿者:明石 2016/05/14(Sat) 05:40 No. 22002
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
ご教示をいただき、誠にありがとうございました。

いつもながら、研ぎ澄まされたエッセンス!

勉強になりました。

sapply についても、勉強させていただきます。

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


単変量と多変量で結果の違う説明変数の解釈
  投稿者:りんご 2016/05/04(Wed) 19:57 No. 21993
青木先生
お世話になります。
どうぞよろしくお願い申し上げます。

縦断研究でリスク因子を求めたいと思い、ロジスティック回帰分析を行いました。
単変量のロジスティック回帰分析では有意にならない変数が、多変量で行うと有意になりました。
このような説明変数はどのように解釈すればよいのでしょうか。
何か他の因子との組み合わせでたまたま重要に見える見かけ上の有意な変数で、意味のないものと考えればよいでしょうか。

Re: 単変量と多変量で結果の違う説明変数の解釈
  投稿者:青木繁伸 2016/05/05(Thu) 06:26 No. 21994
この掲示板でも何回も取り上げましたが,単変量のロジスティック回帰分析で有意な変数だけ集めても,多変量ロジスティック回帰分析にはなりません。

回帰分析に限らず,データ解析は分析に用いた変数セットにおける分析結果になります。単回帰分析で有意であっても,他に重要な変数があった場合,その重要な変数と一緒に多変量回帰分析すれば,単回帰で有意であった変数が,取るに足らない変数として捨てられてしまうのは当然でしょう。重要な変数が欠落した多変量解析の結果は,重大な欠陥があるでしょう(単回帰分析はまさにそのような危険性のある分析法。そんな分析結果を重要視する意味がわからない)。
逆の例もあります。例えば,2群で平均値の同じ説明変数だけを使ってを使った判別関数で,は,そんな変数が無意味であるのは自明でしょう。しかし,その平均値が同じになる説明変数を使った方が判別率が高くなる例があるのです。説明図

> 単変量のロジスティック回帰分析では有意にならない変数が、多変量で行うと有意になりました。

まず,単変量の結果を述べて,次に,多変量の結果を述べるというのが医学論文のお手本としてよく取られる筋書きですが,統計学的に問題があるのは明らかなので,駆逐されるべきと思っております。

Re: 単変量と多変量で結果の違う説明変数の解釈
  投稿者:りんご 2016/05/05(Thu) 11:35 No. 21997
青木先生

ご回答いただき、誠にありがとうございます。
単変量の結果を述べて、次に多変量の結果を述べるということがあまり意味がないということについてとても勉強になりました。

今回の質問について、併せて教えていただきたいのですが、
多変量である変数が有意になった場合、単変量で有意でなかったという結果を無視して、多変量の変数セットの中ではその変数が重要な意味を持つとシンプルに解釈してよいでしょうか。

このような単変量:有意でない⇔多変量:有意 という乖離が起きるのは、多変量の変数セットの中に重要な変数が欠落していることを示唆していると考えるべきでしょうか。

すみませんが、ご教示いただけますと幸いです。
どうぞよろしくお願い申し上げます。

Re: 単変量と多変量で結果の違う説明変数の解釈
  投稿者:青木繁伸 2016/05/05(Thu) 17:19 No. 21998
> 多変量である変数が有意になった場合、単変量で有意でなかったという結果を無視して、多変量の変数セットの中ではその変数が重要な意味を持つとシンプルに解釈してよいでしょうか。

> このような単変量:有意でない⇔多変量:有意 という乖離が起きるのは、多変量の変数セットの中に重要な変数が欠落していることを示唆していると考えるべきでしょうか。

繰り返しますが,「データ解析は分析に用いた変数セットにおける分析結果になります」ということに尽きます。

分析に用いた変数セットが妥当かどうかは,あなた以外の第 3 者にはわからない。

Re: 単変量と多変量で結果の違う説明変数の解釈
  投稿者:りんご 2016/05/05(Thu) 20:54 No. 21999
青木先生

ご教示いただき、どうもありがとうございました。
多変量解析の限界について、とても勉強になりました。
ご回答に深謝致します。


Rのソート(複数キー)
  投稿者:奈々瀬 2016/05/04(Wed) 10:26 No. 21992
青木先生、
お世話になります、奈々瀬と申します。
どうぞよろしくお願いいたします。

−−−
複数キーを指定する、Rのソートについてご教示をいただけたら助かります。

ネットで調べましたら、sort.list関数を利用する例が幾つかありました。
私が所望する記載がありませんでしたので、投稿させていただきます。

添付ファイルに、簡単なサンプルデータをお示しをします。

このデータで、例えば、以下の3つの項目
「経済力」「性格」「容姿」の優先順に降順ソートした場合の
Rプログラムを例示していただけると助かります。

つまり、
まずは「経済力」で降順ソートし、
「経済力」の値が同じ場合には、「性格」で降順ソートし、
さらに、「性格」の値が同じ場合には、「容姿」で降順ソートします。

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


Re: Rのソート(複数キー)
  投稿者:青木繁伸 2016/05/05(Thu) 06:57 No. 21995
order 関数を使います(詳しい説明は ? order で)

order の引数で,優先順に指定します。降順でソートする変数は - を付けてやればよいですね。
d = data.frame(
経済力=c(2,1,2,3,2,1,1,4,4,1,3,2,2,1,2),
性格 =c(5,4,5,5,5,5,4,5,5,3,5,5,4,4,5),
容姿 =c(5,2,4,4,3,5,3,2,4,4,4,4,3,4,5))
d[order(-d$経済力, -d$性格, -d$容姿),]
結果
経済力 性格 容姿
9 4 5 4
8 4 5 2
4 3 5 4
11 3 5 4
1 2 5 5
15 2 5 5
3 2 5 4
12 2 5 4
5 2 5 3
13 2 4 3
6 1 5 5
14 1 4 4
7 1 4 3
2 1 4 2
10 1 3 4
plyr パッケージの arrange 関数を使うというのもありますが,わざわざそんなもの使わなくても。

御礼: Re: Rのソート(複数キー)
  投稿者:奈々瀬 2016/05/05(Thu) 08:15 No. 21996
青木先生、
お世話になります、奈々瀬と申します。

お休みのところ、ご教示をいただき、大変に助かりました!

先生がお書きになられたように、plyrパッケージ、dplyrパッケージを
使う記事を検索で見つけましたが、
Rだから、きっとorder関数を利用するだけで対応できるはずだ、
という強い期待感がありましたので、投稿させていただきました。

やはり、Rです。

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

失礼いたします。


御教示下さい。
  投稿者:初心者 2016/05/02(Mon) 10:39 No. 21989
青木先生、よろしくお願い致します。
私は認知症の発症のリスクとして、1000人以上を対象とした疫学調査の結果から、高血圧や糖尿病、うつなどの有無が発症のリスクとなるかどうか、統計学的手法をを用いて結果を出すことを考えています。その場合の統計学的手法ですが、ロジスティック回帰分析を使うのが妥当か、それともハザード比を使えばいいのかで悩んでいます。初歩的な御質問で大変恐縮ですが、御教示の程よろしくお願い致します。

Re: 御教示下さい。
  投稿者:青木繁伸 2016/05/03(Tue) 10:04 No. 21990
「ハザード比」というのは,比例ハザードモデルを使うということでしょうね。

比例ハザードモデルは,イベント発生までの時間が観測されていなければ使えません。
ロジスティック回帰分析は,時間変数は使いません。

Re: 御教示下さい。
  投稿者:初心者 2016/05/03(Tue) 11:40 No. 21991
青木先生、ありがとうございました。
ロジスティックス回帰分析で統計をかけてみます。
今後ともよろしくお願い致します。


R ネットワークグラフのエッジリストの集計
  投稿者:明石 2016/04/29(Fri) 11:46 No. 21986
青木先生 様、
お世話になります、明石と申します。

先日は、Rの、表〜リストの相互変換についてご教示をいただき、
誠にありがとうございました。

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

−−−
今回は、集計について、ご教示をいただければ助かります。
どうぞ、よろしくお願いいたします。

添付の画像ファイルをご覧ください。

右側に、ネットワークグラフがあります。(大規模グラフ)

このネットワークのエッジリストを、左側の3つ組で表しています。
 lhsは、エッジの始点ノード番号
 rhsは、エッジの終点ノード番号
 weightは、エッジの重み

この3つ組(ファイル名は net.csv)を読み込み、
各ノードごとに、接続している全てのエッジの重みの和を計算し、
以下のように出力したいと考えています。


ノード番号 エッジ重み和
1 17.613908
2 26.414958
3 14.171645
4 21.339434
5 22.665279
6 22.433722
7 14.537605
8 10.936302
9 25.582985
10 25.582985
・・・・・・

このRプログラムを作成しましたが、
ご専門家から見れば、Rの利点を活かして、まだまだ改良の余地があると思われます。

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

dat <- read.table("net.csv", header=TRUE, sep=",")

v1 <- unique(dat$lhs)
v2 <- unique(dat$rhs)
v <- union(v1,v2)

n <- length(v)

w <- numeric(n)

for (i in 1:n) {
s1 <- is.element(dat$lhs, v[i])
s2 <- is.element(dat$rhs, v[i])
s3 <- s1 | s2
d <- dat[s3,]
w[i] <- sum(d$weight)
}

df <- data.frame(ノード番号=v, エッジ重み和=w)


Re: R ネットワークグラフのエッジリストの集計
  投稿者:青木繁伸 2016/04/29(Fri) 17:22 No. 21987
lhs と rhs を添え字として,エッジの重みを要素とする行列 mat を作り,行和ベクトルとと列和ベクトルの和が「重み和」になるので,以下のようなプログラムができる
dat = data.frame(lhs=c(3,3,2,5,6), rhs=c(2,5,5,4,2), weight=c(0.1, 0.2, 0.2, 0.3, 0.5))

v = sort(unique(c(dat[, 1], dat[, 2])))
m = max(v)
mat = matrix(0, m, m)
subscript = as.matrix(dat[, 1:2])
mat[subscript] = dat[, 3]
w = colSums(mat) + rowSums(mat)
data.frame(ノード = v, 重み和 = w[v])

mat
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0.0 0 0.0 0.0 0
[2,] 0 0.0 0 0.0 0.2 0
[3,] 0 0.1 0 0.0 0.2 0
[4,] 0 0.0 0 0.0 0.0 0
[5,] 0 0.0 0 0.3 0.0 0
[6,] 0 0.5 0 0.0 0.0 0
結果の出力順は,あなたのプログラムの v <- union(v1,v2) を, v <- sort(union(v1,v2)) とすれば同じになる。

mat 行列の疎の具合によるが,場合によってはかなりの処理速度が得られるはず。

【御礼】 Re: R ネットワークグラフのエッジリストの集計
  投稿者:明石 2016/04/29(Fri) 18:12 No. 21988
青木先生 様、
お世話になります、明石と申します。

お休みのところ、ご丁寧なご教示を頂戴し、誠にありがとうございます。
追試し、調査して、使いこなしたいと思います。

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


共分散分析の「等分散」の意味
  投稿者:Kei 2016/04/20(Wed) 13:07 No. 21976
基本的な質問で恐れ入りますが、教えて下さい。
共分散分析の前提条件である「等分散性」とは、F検定などと同じで各群の目的変数の分散が等しいという意味なのでしょうか?よろしくお願いします。

Re: 共分散分析の「等分散」の意味
  投稿者:鈴木康弘 2016/04/23(Sat) 07:21 No. 21977
そのとおりだと思います。

Re: 共分散分析の「等分散」の意味
  投稿者:Kei 2016/04/26(Tue) 14:57 No. 21984
お返事ありがとうございます。これまであまり統計に触れてこなかったもので、調べていくうちに混乱してしまいました。
もう一つ質問をさせて下さい。
共分散分析を進めて全群を合わせた共通の回帰式の傾きを計算した後に、群ごとの回帰式の傾きが共通の回帰式の傾きに合うように各郡内でデータを調整する。次に郡内で調整したデータで群ごとに回帰分析を行う。ここで、その残差の散らばりが均一(等分散性)である必要があると思ったのですが、この理解は間違っているのでしょうか。よろしくお願いします。

Re: 共分散分析の「等分散」の意味
  投稿者:鈴木康弘 2016/04/27(Wed) 08:14 No. 21985
過去の青木先生のご発言
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc040/03957.html

によれば、SPSSでは回帰直線からの「誤差の分散の等質性検定」をするようですね。

 でも普通に共分散分析の「等分散性」というと、単なる目的変数の分散じゃないのかな..?(だんだん自信がなくなってきた。)


【R】 表 〜 リスト の相互変換
  投稿者:明石 2016/04/24(Sun) 10:50 No. 21978
青木先生、
お忙しいところを失礼いたします、明石と申します。

過日は、Rについて、丁寧にご教示をいただきまして、誠にありがとうございました。

添付の画像ファイルについて、ご教示をいただきたく、よろしくお願いいたします。

左側は、2部グラフ(人〜本)の、表データです。

人(行)が本(列)を5段階で評価した値(セル)であり、大規模データです。

この表で、評価値(セル)が0でない、人(行)と本(列)の組み合わせを取り出して、
右側の、2部リスト(人〜本)を作成する処理を(1)とします。

逆に、
右側の2部リスト(人〜本)から、左側の2部グラフ(人〜本)を作成する処理を(2)とします。

(1)、(2)の処理プログラムを作成しましたが、
ご専門家から見れば、まだまだ改良の余地が多いにあると思われますので、
改良の観点、プログラム例をお示しいただけると、勉強になります。

お手数をお掛けいたします。
どうぞ、よろしくお願いいたします。

#------------------------------------------------
# (1) 2部グラフ → 2項リスト
#------------------------------------------------

dat <- read.table("hito_book_table.csv", header=TRUE, sep="," , row=1)

dat <- as.matrix(dat)

m <- nrow(dat)
n <- ncol(dat)

hito <- numeric(m*n)
book <- numeric(m*n)
score <- numeric(m*n)

cnt <- 0

for (i in 1:m) {
for (j in 1:n) {
if (dat[i,j] != 0) {
cnt <- cnt + 1
hito[cnt] <- rownames(dat)[i]
book[cnt] <- colnames(dat)[j]
score[cnt] <- dat[i,j]
}
}
}

length(hito) <- cnt
length(book) <- cnt
length(score) <- cnt

df <- data.frame(hito, book, score)

write.csv(df, "hito_book_list.csv")

#------------------------------------------------
# (2) 2項リスト → 2部グラフ
#------------------------------------------------

dat <- read.table("hito_book_list.csv", header=TRUE, sep=",", stringsAsFactors=F)

r <- nrow(dat)

hito <- dat$hito
s1 <- unique(hito)
m <- length(s1)

book <- dat$book
s2 <- unique(book)
n <- length(s2)

mx <- matrix(0, nrow=m, ncol=n)

rownames(mx) <- s1
colnames(mx) <- s2

for (i in 1:r) {
x <-which(is.element(s1, dat$hito[i]))
y <-which(is.element(s2, dat$book[i]))
mx[x, y] <- dat$score[i]
}

write.csv(mx, "hito_book_table.csv")


Re: 【R】 表 〜 リスト の相互変換
  投稿者: 2016/04/24(Sun) 19:58 No. 21979
#------------------------------------------------
# (1) 2部グラフ → 2項リスト
#------------------------------------------------

data.frame(
hito = rep(rownames(dat), each = ncol(dat)),
book = rep(colnames(dat), nrow(dat)),
hyouka = as.vector(dat)
)

#------------------------------------------------
# (2) 2項リスト → 2部グラフ
#------------------------------------------------

xtabs(hyouka ~ hito + book, data = dat)

でいかがでしょうか?

Re: 【R】 表 〜 リスト の相互変換
  投稿者:明石 2016/04/24(Sun) 20:28 No. 21980
荒 様;

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

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

コードは、あまりにも研ぎ澄まされているため、にわかには理解できないでおります。

追試しながら、調べながら、理解させていただきます。

まずは、御礼まで。
ご親身にありがとうございました。

Re: 【R】 表 〜 リスト の相互変換
  投稿者: 2016/04/24(Sun) 20:48 No. 21981
明石さん

1番目のコードですが、行と列を入れ替えないと駄目でした。

data.frame(
hito = rep(rownames(dat), each = ncol(dat)),
book = rep(colnames(dat), nrow(dat)),
hyouka = as.vector(t(dat)) # <= t(dat) に修正
)

評価が0の行が不必要な場合にはこの後にsubset関数で除去して下さい。

rep関数の使用方法さえ分かれば理解できると思います。
また、datの行と列の両方ともラベルがついていないと失敗します。
おそらく他にもっと簡単なやり方があるかもしれません。

Re: 【R】 表 〜 リスト の相互変換
  投稿者:青木繁伸 2016/04/25(Mon) 19:02 No. 21982
(1) 2部グラフ → 2項リスト は,以下のようにも書けますね。
速いかどうかは見ていません。

dat2 = t(dat)
index = which(dat2 > 0, arr.ind=TRUE)
data.frame(人=colnames(dat2)[index[,2]], 本=rownames(index), 評価=dat2[index])

【御礼】 Re: 【R】 表 〜 リスト の相互変換
  投稿者:明石 2016/04/25(Mon) 21:33 No. 21983
青木先生 様、荒 様;

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

今回も、親身に、ご相談にのってくださり、誠にありがとうございました。

荒様からご教示いただきました内容は、追試し、調査して、理解できました。
今後、使いこなしができます。
ありがとうございました。

青木先生からご教示いただきました内容は、にわかには理解できませんので、
じっくりと勉強させていただき、技をマスターしたいと思います。
ありがとうございました。


R-TipsのHPについて
  投稿者:さくらとお城 2016/04/18(Mon) 12:59 No. 21975
こことは直接関係のない内容ですが,もしご存じの方がいれば,ぜひ教えてください。

R-Tipsの下記URLにアクセスできなくなったようです。
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html

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


Rのメモリ空間のサイズを明示的に設定するコマンド
  投稿者:奈々瀬 2016/04/16(Sat) 10:47 No. 21972
青木先生、
奈々瀬と申します、
Rに関する投稿が多く見られますので、ご教示を頂戴できれば助かります。

Rの64bit版(Windows)を利用しています。
メモリは32GBを積んでいますので、実メモリは十分にあると思いますが、
メモリ不足のエラーが出ることがあります。

Rでも、Javaのように、メモリサイズを明示的に指定できればよいなと
思って調べていますが、見つかりません。

知人に聞いても、誰も知りません。

Rのメモリ空間のサイズを明示的に設定するコマンドをご存知でしたら、
ご教示を頂戴できれば助かります。

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

Re: Rのメモリ空間のサイズを明示的に設定するコマンド
  投稿者:青木繁伸 2016/04/17(Sun) 09:12 No. 21973
Windows は使用できないので,よくわかりませんが,「R windows メモリ」で検索すると色々出てくるようですね。
http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1474631834
などはいかが?

Re: Rのメモリ空間のサイズを明示的に設定するコマンド
  投稿者:奈々瀬 2016/04/17(Sun) 18:06 No. 21974
青木先生、
奈々瀬と申します、ご教示を頂戴し、誠にありがとうございます。

コマンドを探すというだけでなく、
問題を広げることで色々な観点があることも分かりました。

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

私が、メモリエラーに遭遇するのは、
NMF(非負行列分解)という手法を適用する場合です。

大規模な行列を利用することが根本原因ですが、
手法の性格上、行列を分割して計算し、その計算結果を統合するという
ことはできませんので、苦慮しております。

ご教示をありがとうござました。


Spearman順位相関係数
  投稿者:Kazuki 2016/04/03(Sun) 12:26 No. 21963
いつも丁寧にご回答いただきありがとうございます。
非常に助かっております。

Spearmanの順位相関係数に関してなのですが、どこを探しても順位相関係数に関する検定力分析の方法が見当たりません。
そもそも存在しないものなのか、同じ相関分析と仮定してPearsonの積率相関係数と同様に算出すべきなのかご教授いただければ幸いです。
効果量は「ロー」であろうとは思うのですが...

お忙しいところ恐れ入りますがよろしくご教授下さい。

Re: Spearman順位相関係数
  投稿者:青木繁伸 2016/04/03(Sun) 19:27 No. 21965
spearman correlation coefficient power で検索すると色々出てきます

最初に出ているのは
https://www.statisticssolutions.com/spearman-correlation-2-tailed/

Re: Spearman順位相関係数
  投稿者:Kazuki 2016/04/04(Mon) 13:14 No. 21966
青木先生、お忙しい中、お返事ありがとうございます。
青木先生にご指定いただきましたサイトの他にpearsonとspearmanでは係数に9%ほどの差があるという趣旨の記述があるサイト(下記2つ目)も検索されたのですが、種々意見があるということでしょうか?

https://www.statisticssolutions.com/spearman-correlation-2-tailed/

The power analysis was conducted in G-POWER using an alpha of 0.05, a power of 0.80, and a large effect size (ρ = 0.5) for a two-tailed test. Because Spearman’s rank correlation coefficient is computationally identical to Pearson product-moment coefficient, power analysis was conducted using software for estimating power of a Pearson’s correlation.

https://www.statstodo.com/CorReg_Exp.php

StatToDo to date cannot find a calculation for sample size or power for Spearman's correlation. However, according to Siegal's book of non parametric statistics, Spearman's is about 91% as efficient as Pearson's in terms of power. In other words, for the same power of 0.8 in a Pearson's correlation, the sample size required for Spearman must use a power of 0.8 / 0.91 = 0.88 to calculate sample size. Similarly, after analysis, the power, as calculated for Pearson, must be multiplied by 0.91 to obtain the power for Spearman.

Re: Spearman順位相関係数
  投稿者:青木繁伸 2016/04/04(Mon) 18:42 No. 21967
所詮,近似・見積なので,どこまで正確に細かくやるかはあまり重要ではなくなってしまうということでしょう。たとえば,サンプルサイズが100か1000かは意味があるけど,100か110かの差は意味がないとかいうことでしょう。

Re: Spearman順位相関係数
  投稿者:Kazuki 2016/04/05(Tue) 19:51 No. 21971
青木先生

お返事ありがとうございます。
また、わかりやすい例えで説明していただき感謝致します。
精進して勉強させていただきます。
今後もまたどうぞよろしくお願い致します。


対数正規分布の期待値
  投稿者:Saito 2016/04/05(Tue) 13:11 No. 21970
お世話になっております。
対数正規分布の期待値が理論解とあわないのですが、どこが間違っているかご指摘頂けないでしょうか。

> ###データの生成###
> set.seed(3)
> a <- 0.2
> b <- -0.2
> id <- 100
>
> x <- seq(1, 10, length=id)
> y <- exp(a + b*x + rnorm(id, 0, sd=4))
>
> ###対数正規分布を仮定して回帰してみる###
> func <- function(par, x, y) {
+ a <- par[1]
+ b <- par[2]
+ sig <- par[3]
+ return(
+ -1*sum(
+ log(
+ (1/(sqrt(2*pi*sig^2)*(y)))*exp(-((a+b*x)-log(y))^2/(2*sig^2))
+ )
+ )
+ )
+
+ }
> res_l <- optim(func, par=c(-0.2, -0.2, 4), x=x, y=y)
>
> ###対数正規分布なのでexp(u+s^2/2)で期待値となるはずだがならない?###
> mean(exp(res_l$par[1] + res_l$par[2]*x + res_l$par[3]^2/2))
[1] 130.4608
> mean(y)
[1] 15.99745

###普通のGLMでやっても合わない###
> res <- glm(log(y) ~ x)
> mean(exp(predict(res)+summary(res)$dispersion/2))
[1] 146.5174

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


Rパッケージ
  投稿者:赤岳 2016/04/01(Fri) 14:54 No. 21958
いつもお世話になっています。
このサイトでRのバグを聞くのはお門違いとは思いますが、ご容赦ください。
MCMCpackをインストールしたいのですが、「graphが利用できない」とか、「MCMCregressを見つけられない」といったメッセージが出て、うまく動きません。
どなたか解決方法を教えていただけませんでしょうか。
使用のR Versionは、R version 3.2.4 Revisedです。
OSは、Windows7です。

> install.packages("MCMCpack")
警告: 依存対象 (dependency) ‘graph’, ‘Rgraphviz’ が利用できません
URL 'https://cran.ism.ac.jp/bin/windows/contrib/3.2/MCMCpack_1.3-5.zip' を試しています
Content type 'application/zip' length 3014557 bytes (2.9 MB)
downloaded 2.9 MB
パッケージ ‘MCMCpack’ は無事に展開され、MD5 サムもチェックされました
ダウンロードされたパッケージは、以下にあります
C:\Users\・(省略)・・\downloaded_packages
> library(MCMCpack)
Error in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]) :
‘graph’ という名前のパッケージはありません
エラー: ‘MCMCpack’ に対するパッケージもしくは名前空間のロードが失敗しました

> y <- c(1,3,3,3,5)
> x1 <- c(1,2,3,5,5)
> x2 <- c(2,3,1,5,2)
> x <- matrix(c(x1,x2),ncol=2,nrow=5)
>
> result3.sim <- MCMCregress(y~x, data = parent.frame(), burnin = 1000, mcmc = 10000)
エラー: 関数 "MCMCregress" を見つけることができません

よろしくお願いします。

Re: Rパッケージ
  投稿者:鈴木康弘 2016/04/03(Sun) 08:37 No. 21962
 MCMCpackなるものは知らないのですが、エラーメッセージからすると、Rgraphvizというパッケージも必要みたいですね。

http://d.hatena.ne.jp/Tanakky/20070727

などを見ると、そのパッケージインストールも一筋縄ではいかないみたい?

Re: Rパッケージ
  投稿者:赤岳 2016/04/03(Sun) 12:34 No. 21964
鈴木先生
レスいただきありがとうございます。
CRANを探してみましたが、見当たらず、下記の文章がありました。
Package ‘Rgraphviz’ was removed from the CRAN repository.
まだ、解決していませんが、取り急ぎお礼申し上げます。

Re: Rパッケージ
  投稿者:青木繁伸 2016/04/04(Mon) 18:49 No. 21968
解決するかどうかはわからないし,古い記事もあるのだけれど,
「Rgraphviz インストール」で検索すると,たくさん出てきますけど
(ちなみに,Mac だと,何の問題もなく動くようだが??)

Re: Rパッケージ
  投稿者:赤岳 2016/04/05(Tue) 12:54 No. 21969
青木先生、鈴木先生
できました!
ありがとうございます。
「Rgraphvizインストール」で検索し、
source("http://bioconductor.org/biocLite.R");biocLite("Rgraphviz")
でインストールできました。
また、MCMCpackも下記のようにバグることなく動きました。
ほんとうに感謝します。ありがとうございました。

> library(MCMCpack)
##
## Markov Chain Monte Carlo Package (MCMCpack)
## Copyright (C) 2003-2016 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
##
## Support provided by the U.S. National Science Foundation
## (Grants SES-0350646 and SES-0350613)


Rの論理ベクトル
  投稿者:明石 2016/03/31(Thu) 21:33 No. 21955
青木先生 様;

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

過日は、Rプログラミングについて丁寧なご教示をいただき、誠にありがとうございました。
改めて御礼を申し上げます。

今回は、Rの論理ベクトルについて、ご教示をいただきたいと思います。
どうぞよろしくお願いいたします。

−−−

以下の2つのファイルを用いた、条件抽出の問題です。

"dat.csv"  ;データファイル

 m 行 n 列のデータフレーム
 行; 事例
 列; 項目(カテゴリカル変数)
    項目名; item_1, item_2, item_3, , , , , item_n
    項目値; "C_1", "C_2", "C_3" などの文字列

"rule.csv"  ;条件抽出ルールを記述したファイル

 r 行 1 列のデータフレーム

 ルールの記述例;
 item_1="C_1" & item_2="C_4"
  item_3="C_2" & item_5="C_4" & item_8="C_8"
 ...............

処理手順は、以下の通りです。

(1) データファイル "dat.csv" の読み込み

dat <- read.table("dat.csv", header=TRUE, sep="," , row=1)

(2) ルールファイル "rule.csv" の読み込み、条件抽出

ループを回して、1つづつルールを取り出し、条件抽出式に変換します。

最終的に、データファイルの各事例に、 r 個のルールを適用した結果の論理ベクトルを得ます。
(「星取表」と呼びます)

(3) 結果ファイル

上記の星取表を出力します。

 m 行 r 列のデータフレーム
 行; 事例
 列; ルール名

先生にお聞きしたいことは、(2) のプログラムです。

(Q1)

以下にお示ししますように、ルールから条件抽出式を組み立てる際に、
どのような文字列操作を行えばよいのか。

ルールに含まれている "" を、どのように扱えばよいのか、苦慮しています。

【ルール】
item_1="C_1" & tem_2="C_4"

【条件抽出式】
dat$item_1=="C_1" && item_2=="C_4"

(Q2)

文字列操作により、条件抽出式を得たとします。

このままでは文字列ですので、論理ベクトルを得ることはできません。

"dat$item_1=="C_1" && item_2=="C_4""

文字列 "dat$item_1=="C_1" && item_2=="C_4"" を、
条件抽出式とした扱うためには、どのようにすればよいのか、苦慮しています。

以上2点、ご教示をいただければ助かります。

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

Re: Rの論理ベクトル
  投稿者:青木繁伸 2016/04/01(Fri) 00:59 No. 21956
dat.csv のテストデータ
item_1,item_2,item_3,item_4,item_5
C_2,C_5,C_5,C_5,C_1
C_4,C_1,C_2,C_5,C_4
C_3,C_4,C_4,C_4,C_3
C_1,C_3,C_5,C_4,C_2
C_5,C_1,C_4,C_1,C_1
C_1,C_4,C_4,C_3,C_1
C_3,C_2,C_3,C_4,C_2
C_5,C_1,C_3,C_2,C_3
C_3,C_2,C_2,C_2,C_2
C_3,C_5,C_1,C_2,C_5
rule.csv のテストデータ
item_1="C_1" & item_2="C_4"
item_3="C_2" & item_5="C_4" & item_2="C_1"
のように
& は && ではなく 1 文字
= も == ではなく 1 文字
で入力されているとする

プログラム
dat <- read.table("dat.csv", header=TRUE, sep=",") # row 引数とは?
m <- nrow(dat)
rule <- readLines("rule.csv") # ベクトルとして読み込む
r <- length(rule)
result <- matrix(logical(m*r), m)
for (i in seq_len(r)) {
r <- gsub("=", "==", rule[i])
result[,i] <- with(dat, eval(parse(text=r)))
}
result
出力
> result
[,1] [,2]
[1,] FALSE FALSE
[2,] FALSE TRUE
[3,] FALSE FALSE
[4,] FALSE FALSE
[5,] FALSE FALSE
[6,] TRUE FALSE
[7,] FALSE FALSE
[8,] FALSE FALSE
[9,] FALSE FALSE
[10,] FALSE FALSE
with(dat, eval(parse(text=r))) が,肝!

【御礼】 Re: Rの論理ベクトル
  投稿者:明石 2016/04/01(Fri) 06:06 No. 21957
青木先生 様;

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

今回も丁寧なご教示をいただき、誠にありがとうございました。

大変に勉強になりました。
技を勉強させていただきます。

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

Re: Rの論理ベクトル
  投稿者:明石 2016/04/01(Fri) 21:28 No. 21959
青木先生 様;

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

ご丁寧にご教示いただき、誠にありがとうございました。
追試を行い、理解できました。

"rule.csv" のテストデータは、以下のようでした。
item_1="C_1" & item_2="C_4"
item_3="C_2" & item_5="C_4" & item_2="C_1"

別の形式として、
Rの連関ルール(aRuleパッケージ)の出力
lhs(条件部) → rhs(結論部)
のlhs(条件部)の場合を読み込んだ場合の検討をしています。

item_1=C_1, item_2=C_4
item_3=C_2, item_5=C_4, item_2=C_1

これを、"rule2.csv" とします。

"rule2.csv" を読み込んで、青木先生のプログラムが利用できるように、
前処理のプログラムを作成しました。

動作は確認できましたが、文字列操作の改善の余地があると思われます。

シンプルに書ける方法をご教示いただければ、勉強になります。
どうぞ、よろしくお願いいたします。

# データファイルの読み込み

dat <- read.table("dat.csv", header=TRUE, sep=",")
m <- nrow(dat)

# ルールファイルの読み込み(連関ルールの lhs部)

v <- readLines("rule2.csv") # ベクトルとして読み込む

r <- length(v)

# 青木先生のプログラムが利用できるように、前処理の追加

qr <- NULL

for(i in 1:r) {
  b <- strsplit(v[i], ", ")
  b <- unlist(b)
  n <- length(b)
  for(k in 1:n) {
    a <- b[k]
    d <- strsplit(a, "=")
   d <- unlist(d)
   d1 <- d[1]
   d2 <- d[2]
   d <- paste(d1, "=", "\"", d2, "\"", sep="")
    b[k] <- d
  }
b <- paste(b, sep="", collapse=" & ")
qr <- rbind(qr, b)
}

rule <- as.vector(qr)

# 青木先生〜

r <- length(rule)

result <- matrix(logical(m*r), m)

for (i in seq_len(r)) {
r <- gsub("=", "==", rule[i])
result[,i] <- with(dat, eval(parse(text=r)))
}

result

# 〜青木先生

result <- result + 0

write.csv(result, "星取表2.csv")

Re: Rの論理ベクトル
  投稿者:青木繁伸 2016/04/02(Sat) 09:11 No. 21960
結果として,動けばどのようにプログラムしてもよいのですが,調べることも含めてプログラムを書く時間が短く,結果としてのプログラムが短ければうれしいでしょうね。

ということで,for の中の r を定義する部分を書き換えます(gsub を使えばよいことが示唆されていましたね)。

r <- gsub("=", "==", gsub('(C_[0-9])', '"\\1"', gsub(",", " &", rule[i])))

, を &
C_* を "C_*" にしてから(* は 0~9)
= を == にします。
gsub ごとに 3 行に分けて書くとわかりやすいかも。

【御礼】 Re: Rの論理ベクトル
  投稿者:明石 2016/04/02(Sat) 12:25 No. 21961
青木先生 様;

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

今回も丁寧なご教示をいただき、誠にありがとうございました。

大変に勉強になりました。
正規表現を使う技を勉強させていただきます。

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


タイトル行の名前の変更
  投稿者:コロン 2016/03/27(Sun) 10:17 No. 21951
いつもお世話になっております。相変わらず基本的な質問で恐縮しております。

タイトル行に

a, b, c, d

という文字列が入ったデータ(data)があるとします。

これを

あ,い,う,え

に変更するのは,colnames(data) <- c("あ", "い", "う", "え")

ですが,a だけを「あ」に変えるやりかたはあるのでしょうか?

いろいろとネットで検索してみましたが,ヒットしませんでしたので,こちらに書かせていただきました。

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

Re: タイトル行の名前の変更
  投稿者:青木繁伸 2016/03/27(Sun) 17:08 No. 21952
個別代入ができますので,
colnames(data)[1] <- "あ"
など
colnames(data) <- sub("a", "あ", colnames(data))
でも目的は達せられるが,やり過ぎ
そのほかいろいろな解があるだろうから,探索してみるのも勉強になるかな?

Re: タイトル行の名前の変更
  投稿者:コロン 2016/03/28(Mon) 15:19 No. 21953
青木先生

お返事ありがとうございます。あるんですね! 
[ ] どんなときに使えるのか,もう少し勉強が必要です。
ありがとうございました。

Re: タイトル行の名前の変更
  投稿者:青木繁伸 2016/03/28(Mon) 17:30 No. 21954
> [ ] どんなときに使えるのか

ベクトルや行列(データフレーム)や配列の要素を指定するときに使います。
右辺にあれば引用,左辺にあれば定義。
colnames(data) は文字ベクトルですので,colnames(data)[1] は,ベクトルの最初の要素を表します。


分散分析のpower analysisについて
  投稿者:Kei 2016/03/23(Wed) 20:41 No. 21941
分散分析で有意差が出なかった研究結果について、post hocでパワーアナリシスをするように求められており質問させていただきます。

標準パッケージに含まれているpower.anova.test(群間分散、群内分散を使用)と
pwrパッケージに含まれているpwr.anova.test(効果量を使用)の結果が大きく異なり
使い方が間違えているのではないかと思い質問させていただきます。

例えば
group <- c("a","a","a","a","a","b","b","b","b","b","c","c","c","c","c")
data <- rnorm(15)
testdata <- data.frame(group,data)
anova(aov(data~group,data=testdata))
として結果で表示される

Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.0252 1.5126 1.1924 0.337
Residuals 12 15.2231 1.2686

にあるgroupのMean Sqを群間分散、ResidualsのMean Sqを群内分散として
用いてよいのでしょうか。この結果から

a <- anova(aov(data~group,data=testdata))
power.anova.test(group = 3, n = 5, between.var = a[1,3], within.var = a[2,3], sig.level = 0.05, power = NULL)

として計算されるPowerと、以下の方法の方法で計算される結果が大きく異なってしまいます

# 効果量es = sqrt(群間の平方和/群内の平方和)

es <- sqrt(a[1,2]/a[2,2])
pwr.anova.test(k = 3, n = 5, f = es, sig.level = 0.05, power = NULL)

ご多忙のところ恐縮ですがご教授いただけると助かります。
よろしくお願いいたします。

Re: 分散分析のpower analysisについて
  投稿者:青木繁伸 2016/03/24(Thu) 13:47 No. 21942
power.anova.test の中で,F 分布の非心度 lambda の計算をするところがあります

lambda <- (k - 1) * n * (between.var/within.var)

pwr.anova.test では,

f <- sqrt(between.var / within.var)
lambda <- k * n * f^2

となっております

(k-1) を掛けるか k を掛けるかの違いになっているわけです

どちらが正しいのか,参考書を全て処分したので,私にはわかりません

Re: 分散分析のpower analysisについて
  投稿者:Kei 2016/03/24(Thu) 18:27 No. 21943
青木先生

早々のお返事ありがとうございます。先生のご返信を見て気づいたのですが、

within.var(群内分散)と
Analysis of Variance Tableにある
residualsのMean Sq(群内平均平方)、Sum Sq(群内平方和)
を私は混同してしまっているかもしれません。

power.anova.testを使うとき私は
between.varに群間のMean Sqを、within.varに群内のMean Sqを
代入しました。

また、効果量を求めるときに
青木先生はf <- sqrt(between.var / within.var)とされていますが
私はf <- sqrt(群間のSum Sq/群内のSum Sq)としました。

このどちらかが間違っていると思うのですが、ご教授いただけませんでしょうか。

Re: 分散分析のpower analysisについて
  投稿者:青木繁伸 2016/03/24(Thu) 19:15 No. 21944
> また、効果量を求めるときに
> 青木先生はf <- sqrt(between.var / within.var)とされていますが
> 私はf <- sqrt(群間のSum Sq/群内のSum Sq)としました。

f <- sqrt(between.var / within.var) は,プログラムにそう書かれていて,pwr.anova.test で実際に行われている演算だからそう書いたまで。

ところで,ここでいうところの between.var,within.var は,データから計算された分散分析表での MSb, MSw とは本来別のものであるということがそもそもの勘違い。

k 群あって,各群のサンプルサイズが同一の n であるとき,各群の平均値を μi,全体の平均値を μ としたとき,between.var は Σ(μi-μ)^2 / k。すなわち,
f = sqrt{ (Σ(μi-μ)^2 / k) / within.var } 。
within.var は見積もられるものですから,いろいろな値をとりうるけど,MSw ではない。

example(power.anova.test) の最初の例で,
power.anova.test(groups=3, n=4, var(c(59, 66, 42)), 12^2, sig.level=0.05)
が挙げられている。12^2 が within.var。これは天下り的に来ているわけではないが MSw とは無関係。
between.var は var(c(59, 66, 42)) とされているが,上の f の定義でいうと var(c(59, 66, 42))*(k-1)/k
lambda = k * n * var(c(59, 66, 42))*(k-1)/k / 12^2
= (k-1) * n * var(c(59, 66, 42)) / 12^2

となる。
一方,pwr.anova.test は,引数に f を与える。
繰り返すが,f = sqrt(var(c(59, 66, 42))*(k-1)/k / 12^2) で,内部では lambda の計算に f^2 が使われる。
lambda = k * n * f^2
= k * n * var(c(59, 66, 42))*(k-1)/k / 12^2
= (k-1)*n*var(c(59, 66, 42)) / 12^2


結局の所,以上を理解した上で使えば,power.anova.test も pwr.anova.test も同じ結果になる。
> power.anova.test(groups=3, n=4, var(c(59, 66, 42)), 12^2, sig.level=0.05)

Balanced one-way analysis of variance power calculation

groups = 3
n = 4
between.var = 152.3333
within.var = 144
sig.level = 0.05
power = 0.5846256

NOTE: n is number in each group

>
> f <- sqrt(var(c(59, 66, 42))*(k-1)/k / 12^2)
> f
[1] 0.8397898
> pwr.anova.test(k=3, n=4, f=f, sig.level=0.05)

Balanced one-way analysis of variance power calculation

k = 3
n = 4
f = 0.8397898
sig.level = 0.05
power = 0.5846256

NOTE: n is number in each group

Re: 分散分析のpower analysisについて
  投稿者:青木繁伸 2016/03/24(Thu) 19:37 No. 21945
以下のように,MSw = n * between.var なので,within.var = MSe / n と指定してやると,データに基づく power の予測になるでしょう
この場合は比例するわけだから,MSb,MSw をそのまま使っても結果は同じになる
> k <- groups <- 3
> n <- 10
> x <- c(42, 44, 58, 48, 61, 62, 47, 52, 46, 50,
+ 43, 40, 56, 49, 54, 39, 42, 47, 42, 54,
+ 37, 62, 59, 37, 58, 27, 56, 50, 26, 50)
> g <- rep(factor(1:k), each=n)
>
> summary(aov(x ~ g))
Df Sum Sq Mean Sq F value Pr(>F)
g 2 141.9 70.93 0.786 0.466
Residuals 27 2436.0 90.22
>
> #######
>
> ( between.var <- var(by(x, g, mean)) )
[1] 7.093333
> ( within.var <- 16)
[1] 16
> power.anova.test(groups=groups, n=n, between.var=between.var, within.var=within.var)

Balanced one-way analysis of variance power calculation

groups = 3
n = 10
between.var = 7.093333  ★★ MSb の 1/n
within.var = 16
sig.level = 0.05
power = 0.713057

NOTE: n is number in each group

>
> ( f <- sqrt(var(by(x, g, mean))*(k-1)/k / within.var) )
[1] 0.5436502
> pwr.anova.test(k=k, n=n, f=f)

Balanced one-way analysis of variance power calculation

k = 3
n = 10
f = 0.5436502
sig.level = 0.05
power = 0.713057

NOTE: n is number in each group

>
> #######
>
> ( within.var <- 9.022)
[1] 9.022
> power.anova.test(groups=groups, n=n, between.var=between.var, within.var=within.var)

Balanced one-way analysis of variance power calculation

groups = 3
n = 10
between.var = 7.093333   ★★ MSb の 1/n
within.var = 9.022 ★★ MSw の 1/n を指定
sig.level = 0.05
power = 0.9285187

NOTE: n is number in each group

>
> ( f <- sqrt(var(by(x, g, mean))*(k-1)/k / within.var) )
[1] 0.7239826
> pwr.anova.test(k=k, n=n, f=f)

Balanced one-way analysis of variance power calculation

k = 3
n = 10
f = 0.7239826
sig.level = 0.05
power = 0.9285187

NOTE: n is number in each group

Re: 分散分析のpower analysisについて
  投稿者:Kei 2016/03/25(Fri) 21:57 No. 21947
ありがとうございます。

f = sqrt{ (Σ(μi-μ)^2 / k) / within.var } と定義したとき、
power.anova.testとpwr.anova.testが求めるlambdaが同一であり、結果も同じになることは理解できました。

また、
データ個数をn、グループ数をkとすると
群間平方和が
Σ((群内平均 - 全体平均)^2 * n)
それを群間の自由度 k – 1で割ったものが群間平均平方なので
群間平均平方 = Σ( (群内平均 - 全体平均)^2 * n) / (k - 1 )

ここで
between.var = var (群内平均) = Σ( (群内平均 - 全体平均)^2) )/ (k - 1 )
なので
群間平均平方 = n * between.var つまり between.var = 群間平均平方 / n
となるところまではなんとか理解できました。

ここで質問なのですが、
within.var に群内平均平方 / n を指定するのはなぜなのですか。
within.var は見積もられるものなので、仮に群内平均平方 / nを指定するということでしょうか。

あと、もう1つ質問なのですが、
いくつかのサイトで異なるeffect sizeの定義を見つけました。
http://www.unt.edu/rss/class/Jon/ISSS_SC/Module009/isss_m91_onewayanova/node8.html
https://en.wikipedia.org/wiki/Effect_size

Eta^2 = SSb /SSt
f^2 = Eta^2 / (1 - Eta^2) もしくは R^2 / (1- R^2) など

例えば以下の式で
k <- groups <- 3
n <- 10
x <- c(42, 44, 58, 48, 61, 62, 47, 52, 46, 50,43, 40, 56, 49, 54, 39, 42, 47, 42, 54,37, 62, 59, 37, 58, 27, 56, 50, 26, 50)
g <- rep(factor(1:k), each=n)
summary(aov(x ~ g))
の結果から

SSb <- 141.9
SSt <- sum((x - mean(x))^2)
Eta.squared <- SSb / SSt
f <- sqrt(Eta.squared / (1 - Eta.squared))

ここから求まったfは
within.var <- 9.022 
f <- sqrt(var(by(x, g, mean))*(k-1)/k / within.var)
で求めたfとかなり異なります。
f <- sqrt(Eta.squared / (1 - Eta.squared))
で求めたfはpwr.anova.testに使用するべきではないということでしょうか。

お忙しいところ恐縮ですがまたご指導下さいますようお願い申し上げます。

Re: 分散分析のpower analysisについて
  投稿者:青木繁伸 2016/03/26(Sat) 05:04 No. 21949
> within.var に群内平均平方 / n を指定するのはなぜなのですか。
> within.var は見積もられるものなので、仮に群内平均平方 / nを指定するということでしょうか。

「within.var = MSe / n と指定してやると,データに基づく power の予測になるでしょう」と書いておきましたが。

> f <- sqrt(Eta.squared / (1 - Eta.squared))
> で求めたfはpwr.anova.testに使用するべきではないということでしょうか。

要するに,oneway ANOVA の power を計算するには,非心 F 分布が必要で,その非心度(power.anova.test, pwr.anova.test では lambda)をどのように設定するかということです。

それを計算するための統計量は between.var, within.var を使う(power.anova.test),および,それを使って計算される f を使う(pwr.anova.test)以外にも,sqrt(Eta.squared / (1 - Eta.squared)) でも,F 値でも何でも計算の元にしてよいのだけれど,非心度として計算される lambda は同じでなければならない。逆に言えば,同じ lambda が得られように計算するなら,どのようなものを与えてもよいということ。pwr.anova.test はそのような要望を満たすために作られた。そして,power.anova.test と同じ答えが出る。

sqrt(Eta.squared / (1 - Eta.squared)) を引数に与えて power.anova.test と同じ答えが出るように lambda を計算して答えを出すプログラム(lambda を計算する部分だけが異なるプログラム)を書いて,それを pwr2.anova.test なりの名前で公表すればよいだけの話(自分で使うだけでも良いのだけど)。

Re: 分散分析のpower analysisについて
  投稿者:Kei 2016/03/27(Sun) 08:48 No. 21950
青木先生

ご指導ありがとうございます。
データに基づく power の予測にはwithin.var = MSe / n と指定するとよい点は理解しました。
今回の研究に関してはこれでパワーアナリシスを行うようにします。

また、
f <- sqrt(Eta.squared / (1 - Eta.squared))
をそのままpwr.anova.testに代入してしまうとlambdaの値が正しく求められないので適切でないということは理解できました。
勉強のためにf <- sqrt(Eta.squared / (1 - Eta.squared))から正しいlambdaを求めてパワーアナリシスを行うようなプログラムも書いてみるようにします。


比率の差の検定統計量
  投稿者:kamo 2016/03/14(Mon) 19:20 No. 21937
抽象的な質問で恐れ入ります。

Z(p.x p.y )=(p.x-p.y)/√((1/n.x +1/n.y )(((x+y)(n-(x+y)))/n^2 ) )と紹介されていましたが、

比率の差の信頼区間の導出式が
P(p.x-p.y-Z_α/2 √((p.x (1-p.x ))/n.x +(p.y (1-p.y ))/n.y )
&#8804;p.x-p.y
&#8804;p.x-p.y-Z_α/2 √((p.x (1-p.x ))/n.x +(p.y (1-p.y ))/n.y ))
=1-α
であることから
検定統計量は、
((p.x )-(p.y ))/√((p.x (1-p.x ))/n.x +(p.y (1-p.y ))/n.y )
であると思うのですが、この理解は正しいでしょうか?
また、一般的に
P{検定統計量の下側α/2点 )≦θ≦検定統計量の上側α/2点)}≧1−α
が成立すると理解しているのですが、こちらの理解も正しいでしょうか?

Re: 比率の差の検定統計量
  投稿者:青木繁伸 2016/03/14(Mon) 19:41 No. 21938
どのページの話ですか?私のページについてのクレームなのか,他のページに書いてあったことについての確認なのか?
「抽象的な質問」ではないにもかかわらず,数式の記述がとてもわかりにくく,理解不能(理解しようとすると頭が痛くなりますが)変数名が何を表すかとかもちゃんと書くべきでしょうし。
数式を示すならば,LaTeX や Word で数式を組んで(LaTeX ならタイプセットできるように。もっともそれができる人ならそうするんだろうけど),そのスクリーンショットを取ってページに張り込んでください。

Re: 比率の差の検定統計量
  投稿者:青木繁伸 2016/03/15(Tue) 11:21 No. 21939
以下のように...
クリックで拡大


Re: 比率の差の検定統計量
  投稿者:kamo 2016/03/15(Tue) 23:18 No. 21940
青木様

言葉足らずの質問から意図をくみ取って頂き、また、丁寧にご回答頂きまして誠にありがとうございます。こちらの掲示板で画像を貼り付けられることに気付かず、元々数式エディタで書いた数式を2次形式に落とし込んで記述させて頂きました。青木様及び閲覧者の方々に見づらい記事となってしまい、申し訳ございませんでした。

また、上で質問させて頂きました検定統計量の式は、複数のテキスト及びWebページで紹介されていたものでございます。この点も記述不足でございました。

質問の内容は、青木様が作成された画像データの通りで、回答にも納得しました。
自分の理解が誤っていないことが分かり、すっきり致しました。
ありがとうございました。


AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:明石 2016/03/06(Sun) 20:27 No. 21922
青木先生 様、
お忙しいところを失礼いたします、明石と申します。

「AIC による,ヒストグラム(度数分布表)の最適階級分割の探索」
 http://aoki2.si.gunma-u.ac.jp/R/AIC-Histogram.html
について、ご教示をいただきたいと思います。
よろしくお願いいたします。

−−−
例示では、ヒストグラムのビンが19本あります。

(1)
最初の2本⇒併合
次の5本 ⇒併合
次の5本 ⇒併合
次の7本 ⇒併合

結果として、4階級になるという理解でよろしいでしょうか?

(2)
ans$df[1:10,]の読み方が分かりません。

AICの値の小さい順に上位10個を表示していると理解していますが、
c1 r c2 の読み方が、分かりません。

出力結果の見方を、教えていただけないでしょうか?
    c1 r c2      AIC
78 2 5 7 487.5261
60 2 5 2 488.7143
105 3 5 1 490.2904
72 2 6 5 490.6599
130 3 4 8 491.1111
59 2 3 2 491.6356
123 3 5 6 492.0030
116 3 4 4 492.4315
55 2 4 1 492.9979
129 3 2 8 493.2800

(3)
ヒストグラムと ans$df[1:10,] を一緒に見ています。

ヒストグラム(19本のビンを、4階級に併合)の区間を知りたい場合には、
ansのどの要素を参照すればよろしいのでしょうか?

(4)
数量データの離散化については、
値等分割、レコード数等分割(四分位など)、平均と分散
などの方法があります。

本方式「AIC による,ヒストグラム(度数分布表)の最適階級分割の探索」
は、数量データの離散化の一つの手法であるという理解でよろしいでしょうか?

ご教示を、よろしくお願いいたします。

Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:青木繁伸 2016/03/07(Mon) 09:31 No. 21923
(1) 結果として、4階級になるという理解

その通り

(2) c1 r c2 の読み方が、分かりません。

プログラム中の注釈を参照
for (c1 in 1:(c-2)) { # 左端で併合する階級数を探索
for (c2 in 1:(c-2)) { # 右端で併合する階級数を探索
:
for (r in 1:(c-c1-c2)) { # 中央付近で併合する階級数を探索

(3) ヒストグラム(19本のビンを、4階級に併合)の区間を知りたい

df 要素ではなく,breaks 要素を見ます。ans$breaks

(4) 数量データの離散化の一つの手法

その通り

【御礼】 Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:明石 2016/03/07(Mon) 09:37 No. 21924
青木先生 様、
お忙しいところを失礼いたします、明石と申します。

ご教示いただき、大変に助かりました。

今まで、EMアルゴリズム(混合分布)を用いて、離散化していましたが、
うまくいかないことが多く、苦慮していました。

心から、心より、御礼を申し上げます。

Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:明石 2016/03/07(Mon) 21:50 No. 21928
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
お世話になります。

追加の質問がございます。
ご教示を何卒どうぞ、よろしくお願いいたします。

(4)
ans$n

> [1] 16 22 11 6 7 9 11 4 3 2 3 2 1 1 0 0 0 0 2

ビン(19本)の高さ、と理解しました。

(5)

ans$breaks

> [1] 0.015000 5.369737 10.724474 16.079211 21.433947 26.788684 32.143421 37.498158 42.852895
[10] 48.207632 53.562368 58.917105 64.271842 69.626579 74.981316 80.336053 85.690789 91.045526
[19] 96.400263 101.755000

ビン(19本)の区間の左端、右端の値、と理解しました。

(6)
.frame': 416 obs.

⇒ 収束計算の繰り返し回数と理解しました。

(7)
併合した結果の4階級の「4」という値を取得するには、どうすればよいのでしょうか?

(8)
4階級を、左から、1,2,3,4 と番号付けしたといたします。
ベクトルxの要素(要素数は100)が、どの階級(1〜4)に属するのかを出力する
Rプログラムを考えていますが、作成できません。
誠に申し訳ございませんが、Rコードを例示していただけると助かります。

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

Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:青木繁伸 2016/03/07(Mon) 22:52 No. 21929
併合後のヒストグラム(の長方形)を描くのは,以下の部分で,その長方形の X 座標は,lo+cmg[i]*delta,lo+cmg[i+1]*delta
なので,
	sapply(1:(m+2), function(i) rect(lo+cmg[i]*delta, 0,            # 併合後のヒストグラム
lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
border="red", col="red", density=15))
それを,cat でコンソールに書き出すとすれば<
sapply(1:(m+2), function(i) { rect(lo+cmg[i]*delta, 0, # 併合後のヒストグラム
lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
border="red", col="red", density=15)
cat(lo+cmg[i]*delta, lo+cmg[i+1]*delta, "\n") # これを追加しますね。カテゴリー数は m+2 かな?
})

> ans <- AIC.Histogram(x, 0.01)
0.015 10.72447
10.72447 37.49816
37.49816 64.27184
64.27184 101.755


【Special Thanks !】 Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の...
  投稿者:明石 2016/03/08(Tue) 05:11 No. 21930
青木先生 様、

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

ご教示いただき、大変に助かりました。

ご教示いただきました内容を基にして、自分の頭でトレースをいたします。

心から、心より、御礼を申し上げます。

Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:明石 2016/03/11(Fri) 16:29 No. 21934
青木先生 様、
お忙しいところを失礼いたします、明石と申します。

お世話になります。
ご教示いただきました内容を追試して、理解をしました。

最適階級の区間(左端の値、右端の値)を、コンソールに出力できるようになりましたが、
関数の戻り値に、加えたいと思い、
先生からご教示いただきましたコードに下記のように加筆しましたが、
うまくいきません。

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

k <- m+2 # 最適階級数

mx <- matrix(0, nrow=k, ncol=2)

sapply(1:(m+2), function(i) { rect(lo+cmg[i]*delta, 0, # 併合後のヒストグラム
lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
border="red", col="red", density=15)

r1 <- lo+cmg[i]*delta
r2 <- lo+cmg[i+1]*delta
cat(r1, r2, "\n")

mx[i,1] <- r1
mx[i,2] <- r2

})

invisible(list(df=df, n=n, breaks=brks, k=k, mx=mx)) # 結果を返す

以上、よろしくお願いいたします。

Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:青木繁伸 2016/03/11(Fri) 16:58 No. 21935
一番少ない修正量で動かすには
mx[i,1] <<- r1
mx[i,2] <<- r2
しかし,<<- はなるべく使わないほうがよい

正統的には
#       k = m+2                          # 使わない     
# mx <- matrix(0, nrow=k, ncol=2) # 使わない
cmg <- cumsum(c(0, c1, rep(r, m), c2))
mx <- sapply(1:(m+2), function(i) { rect(lo+cmg[i]*delta, 0, # 結果を mx に代入
lo+cmg[i+1]*delta, mean(p[(cmg[i]+1):(cmg[i+1])]),
border="red", col="red", density=15)
r1 <- lo+cmg[i]*delta
r2 <- lo+cmg[i+1]*delta
cat(i, r1, r2, "\n")
return(c(r1, r2)) # sapply から返される値
})
mtext(paste("AIC =", round(AIC, 2)), side=3, line=-1.2,
at=lo+0.8*c*delta)
invisible(list(df=df, n=n, breaks=brks, k=ncol(mx), mx=t(mx)))


【御礼】 Re: AIC による,ヒストグラム(度数分布表)の最適階級分割の探索
  投稿者:明石 2016/03/11(Fri) 20:40 No. 21936
青木先生 様、

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

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

ご教示いただきました内容を基にして、自分の頭でトレースをいたします。

心から、心より、御礼を申し上げます。


期待されるカプランマイヤー曲線?
  投稿者:Gtaka 2016/03/10(Thu) 02:10 No. 21931
いつも参考にさせていただいております。
初めて質問致します。

ある前癌病変の経過観察中の発癌を見るための後ろ向きコホート研究を考えています。コントロールとなる自前のデータがないため、地域がん登録の罹患データを用いたいと考えました。
SIRを計算するといったことはできるようになったのですが、例えば以下の論文にあるように、一般人口の罹患率から期待されるカプランマイヤー曲線を描いて自データの曲線と重ねて表示したいと思いました。
自データからのnumber at riskと罹患率を掛けあわせてイベント数を1年毎に出して、累積発癌率を計算していくのかなと何となく想像できるのですが、具体的な方法が分からず困っています。Rなどを使ってできるものなのでしょうか、または参考となる文献などありましたら御教示頂けませんでしょうか?
分かりにくい質問でしたら申し訳ありません。よろしくお願い致します。


"Time-specific expected events, again based on rates from the Rochester general population, were calculated for each 1-yr follow-up interval, accounting for age and calendar period
as described above. A modified Kaplan-Meier approach was then used to cumulate expected CRC incidence over these 1-yr periods and the resulting expected event curve was smoothed using linear interpolation."

Paul J. Limburg, et al. Clinically Confirmed Type 2 Diabetes Mellitus and Colorectal Cancer Risk: A Population-Based, Retrospective Cohort Study, Am J Gastroenterol 2006;101:1872-1879

Re: 期待されるカプランマイヤー曲線?
  投稿者:青木繁伸 2016/03/10(Thu) 22:09 No. 21932
その論文は未だ見ることは出来ないのですが,

> 自データからのnumber at riskと罹患率を掛けあわせてイベント数を1年毎に出して、累積発癌率を計算

というのは,単なる指数曲線ということでいいのですか?

Re: 期待されるカプランマイヤー曲線?
  投稿者:Gtaka 2016/03/11(Fri) 03:19 No. 21933
青木先生

すみません、変なことを書きました。

自データ:経過観察を始めた暦年、打ち切り・イベントが発生した暦年、生年月日
がん登録データ:年齢階級別、毎年の罹患率

例えば、自データに2000年から経過観察を開始した50歳の人が10人いたとします。
この10人のコホートについて、がん登録データから1年毎の期待罹患数を計算します。
2000年50歳の期待罹患数A= 10*(2000年50歳の罹患率)
2001年51歳の期待罹患数B= (10-A)*(2001年51歳の罹患率)
2002年52歳の期待罹患数C= (10-A-B)*(2002年52歳の罹患率)
・・・
同様にして、2001年から開始した60歳の人が10人
2001年60歳の期待罹患数D= 10*(2001年60歳の罹患率)
2002年61歳の期待罹患数E= (10-D)*(2002年61歳の罹患率)

2002年で終了、これで全員分だとすると
1年目の期待罹患数=A+D
2年目の期待罹患数=B+E
3年目の期待罹患数=C
最初のnumber at riskを10+10人として、累積罹患率を求める。
という具合に計算できるのかな、と考えた訳です。

この考え方が正しいのか、そもそもこれが一般的な方法なのか、もしそうであればこれを計算するプログラムがあるのか(自分はプログラミングができません)。
拙い質問で恐縮ですがお知恵を貸していただけないでしょうか。
よろしくお願い申し上げます。


相関係数について
  投稿者:平林 2016/03/05(Sat) 07:12 No. 21919
自動的に頭部MRI画像の病変体積を求めるツールと手でトレースして求めた病変体積の相関係数をもとめています。病変体積は、1〜100まで幅がありますが、平均は10程です。例えば体積値が小さいところでは、マニュアルトレースが1で、自動が0.5と50%の正答率ということはよくあります。逆に100に近いところでは、100と90といったように大きくはずれません(絶対値としては大きくずれているといえるかもしれませんが)。相関係数というのは、この場合体積が大きいところと小さいところのどちらにより左右されるのでしょうか?大きいところで大きくずれていなければそれにひっぱられて相関係数もよくなるのか、小さいところでずれているので、相関係数も悪く出るのか?ご教授よろしくお願いいたします。

Re: 相関係数について
  投稿者:青木繁伸 2016/03/05(Sat) 09:31 No. 21920
r = Σ(Xi-mean.X)(Yi-mean.Y) / sqrt(Σ(Xi-mean.X)2 Σ(Yi-mean.Y)2) ですから,0 を基準にして大きい小さいというのではなく,平均を基準にして大きいか小さいかを考えるのです。

あなたのデータは,1〜100の範囲といいながら,平均値は10程度ということで,正規分布から大きく外れている。
マニュアルと自動の測定値が (1, 0.5) なら平均値からの差は (-9, -9.5) であるのに対し,(100, 90) なら平均値の差は (90, 80) ということなので,後者の方が相関係数には大きい影響を与えることになるでしょう。

Re: 相関係数について
  投稿者:かい 2016/03/07(Mon) 12:42 No. 21925
上記の質問に相乗りで質問させていただきます.

もしも,分布の形状が対数正規分布とした場合,両方の変数を対数変換して相関係数を計算したほうが妥当性があるということになるのでしょうか?

Re: 相関係数について
  投稿者:平林 2016/03/07(Mon) 13:26 No. 21926
青木先生ありがとうございます。正規分布していないので、スピアマンの相関係数でやれば、あるデータにひっぱられずに相関係数をもとめられるということになりますかね?

Re: 相関係数について
  投稿者:青木繁伸 2016/03/07(Mon) 18:37 No. 21927
> 分布の形状が対数正規分布とした場合,両方の変数を対数変換して相関係数を計算したほうが妥当性がある

多分そういうことになるでしょう


比率の検定の問題
  投稿者:高木洋介 2016/02/17(Wed) 11:01 No. 21910
比率の検定における質問です。

ある同一LOT製品を別のラインで作ったときに

Aライン 25/7580 不良発生
Bライン 0/1500  不良発生)(不良0)

の結果。ラインの人間は「1%の有意水準でラインに差があるとはいえない、5%では有意であるが」と言い張ります。

確かに検定をすると、χ2でもフィシャーの正確検定でも α=0.01〜0.05の間になり、言っていることは合っているのですが、

αリスクを低くしても、βリスクが高いのではと考えます。

実際は物理的に工程差があり不良が多い理由があるのですが、この検定の妥当性はないことが計算で示せるでしょうか?

Re: 比率の検定の問題
  投稿者:鈴木康弘 2016/02/19(Fri) 07:10 No. 21911
 ここの

http://aoki2.si.gunma-u.ac.jp/R/power_prop_test2.html

 で実際に検出力を計算してみればいかがかと。

Re: 比率の検定の問題
  投稿者:高木洋介 2016/02/19(Fri) 13:47 No. 21912
ありがとうございます。

実行して見ましたが

Nc 7605
Nt 1500
Pc 0.003298153
Pt 0



sig.level 0.01・・・power=0.216214
sig.level 0.05・・・power=0.726417
sig.level 0.1 ・・・power=0.9051945

有意水準1%の結果で有意差あるとはいえないが、事後の検定ではpowerが小さく(βリスク 80%)有意差あるものを見逃すリスクが大きかった。有意水準をもう少し高くして検定をすべきであったという解釈でしょうかね。

ちなみに片側検定もできるのでしょうか?

Re: 比率の検定の問題
  投稿者:鈴木康弘 2016/02/20(Sat) 07:10 No. 21913
 検出力が0.7以上の結論を言いたいなら、有意水準が0.05あたり(かもっと甘く)でないといけない、とは言えるでしょうね。

 > ちなみに片側検定もできるのでしょうか?

 そういう難しい話は、作者の青木先生に丸投げ。

Re: 比率の検定の問題
  投稿者:青木繁伸 2016/02/20(Sat) 18:48 No. 21914
原理的に考えましょう。たとえばね,
> power.prop.test(p1=0.4, p2=0.3, n=100)

Two-sample comparison of proportions power calculation

n = 100
p1 = 0.4
p2 = 0.3
sig.level = 0.05
power = 0.3155744
alternative = two.sided

NOTE: n is number in *each* group

> power.prop.test(p1=0.4, p2=0.3, n=100, sig.level=0.025, alternative="one.sided")

Two-sample comparison of proportions power calculation

n = 100
p1 = 0.4
p2 = 0.3
sig.level = 0.025
power = 0.3155744
alternative = one.sided

NOTE: n is number in *each* group

ということで,両側検定での sig.level=0.05 は片側検定で sig.level=0.025 を指定するのと同じこと。(図に描けば明らか)
ということで,
両側検定
> power.prop.test2(Nc=7580, Nt=1500, Pc=25/7580, Pt=0, sig.level=0.01)
[1] 0.216613
に対して,片側検定
> power.prop.test2(Nc=7580, Nt=1500, Pc=25/7580, Pt=0, sig.level=0.01/2)
[1] 0.09619114
でよいのかな?(プログラムが正しく作られていればということだけど)
いや,間違えた
> power.prop.test2(Nc=7580, Nt=1500, Pc=25/7580, Pt=0, sig.level=0.01*2)
[1] 0.4118825
だな?

Re: 比率の検定の問題
  投稿者:高木洋介 2016/02/24(Wed) 11:14 No. 21918
御丁寧なご説明ありがとうございました。実際問題での考え方に活用させていただきます。


検出力分析について
  投稿者:Kazuki 2016/02/22(Mon) 13:16 No. 21915
私の知識不足かもしれませんが質問させていただければと思います。
統計に関しては扱うのみで数学的な内容に関しては詳しくありません。

事後の検出力分析に関してなのですが、様々な著書を読むと事後の検出力分析でもαエラーの値は仮説上の有意水準(α<0.05)で入力するように記載されています。
1つの統計学検討(1つの研究)に関する事後の検出力を分析するのであればαエラーは実際に出力された値を入力することで、その検討に関する実際の検出力を算出できるものと思うのですが...それは間違った考え方でしょうか。

完結に述べますと
 検定方法:Pearson correlation coefficient(有意水準p<0.05)
 結果:n=20, r=0.63, p<0.005
に対して事後の検出力分析を行う際に、αエラーの値は0.05(仮説)とするべきか0.005(実際の値)とするべきか、どちらが正しいのかお聞きしたいです。

お忙しいところお手数ですが、ご教授いただければ幸いです。

Re: 検出力分析について
  投稿者:青木繁伸 2016/02/22(Mon) 14:32 No. 21916
「0.005(実際の値)」は有意水準ではなく有意確率です。
また,「そのような条件で検定を繰り返したときに...」ということからいっても,有意水準を書くべきでしょう。

Re: 検出力分析について
  投稿者:Kazuki 2016/02/22(Mon) 19:30 No. 21917
青木先生

お忙しいところどうもありがとうございます。
とても簡潔で理解しやすいご説明をいただきありがとうございました。
用語の学習からもう一度見直してみようと思います。

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