青木先生　　皆様

お世話になります。いろいろ調べたのですが、自信がないので、質問させていただきたいと存じます。

添付のサイトのP182 表88のような検定を行いたいのですが、検定方法がわかりません。

質問1：　直接法とは何でしょうか？

表88には、注釈に、注１）年齢階級（20−29歳，30−39歳，40−49歳，50−59歳，60−69歳，70歳以上の６区分）と世帯員数（１人，２人，３人以上世帯の３区分）での調整値。割合に関する項目は直接法，平均 値に関する項目は共分散分析を用いて算出。 む

と書いてあります。

また、群間差の検定について、注２）世帯の所得について，多変量解析（世帯の所得額を当該世帯員に当てはめて，割合に関する項目はロジスティック回帰分析，平均値に関する項目は共分散分析）を用いて600万 円以上を基準とした他の２群との群間比較を実施。 とあり、

質問2: 群間差の検定は、200万円未満VS600万円以上、200-600未満VS600万円以上の2回、群間差の検定を行えばよろしいでしょうか？T検定ではご法度なので、躊躇しております。　　

また、報告書の冒頭には、
（１）所得と生活習慣・食生活に関する分析 　年齢（20−29歳，30−39歳，40−49歳，50−59歳，60−69歳，70歳以上 の６区分）と世帯人数（１人，２人，３人以上世帯の３区分）で調整を行った後， 割合に関する項目は多変量ロジスティック回帰分析，平均値に関する項目は共分 散分析を行った

群間比較について，年齢調整を行い分析を行ったものについては，割合に関す る項目はCochran-Mantel-Haenszel検定，平均値に関する項目は，共分散分析を 行った。

とあります。

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

サイト情報がありません（家マークでのリンクではなく，本文中に記載してください。）

というわけで，詳細は分かりませんが，直接法については以下のごとく。

直接法とは年齢調整死亡率の計算に用いられる方法。ほかに，間接法もある。
死亡率に限らず，率の調整に使われることもある。
直接法についての説明は，
http://www.takenet.or.jp/~hayakawa/u-tanqa04.htm
http://www.kenko.pref.yamaguchi.lg.jp/img/kenko21/download/h26map/h26kaisetsu-new.pdf
など，
標準化死亡比の検定は，
https://www.niph.go.jp/soshiki/07shougai/datakatsuyou/data/h22smr-kinki.pdf
などを参照のこと。

ここで質問しても詳細な（適切な）回答が得られるとは限りません。
著者に直接お問い合わせになるのが筋でしょう。

青木先生

お返事をどうもありがとうございました。

頂いたリンク、勉強してみます。

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

青木先生　様；

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

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

---------------------------
受験者の生年月日と試験日から、試験日での年齢を計算したいと思います。

受験者データはデータフレームdatです。
生年月日の変数名はbirthです。日付型です。

試験日の変数名はtestです。日付型です。
test <- as.Date("2018-9-13")

以下の繰り返し処理で、受験者の試験日での年齢yearを計算できることを確認しました。
for (i in 1:nrow(dat)) {
dat$year[i] <- length(seq(dat$birth[i], test, by = "year")) - 1
}

受験者データが大規模ですので、繰り返し処理を高速化したいと思います。

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

50000 人分のテストデータを作ってやってみました

test <- as.Date("2018-9-13")
dat <- data.frame(birth=as.Date(Sys.Date()-1:50000)) # 50000 人分のテストデータ

system.time({
for (i in 1:nrow(dat)) {
	dat$year[i] <- length(seq(dat$birth[i], test, by = "year")) - 1
}
})
   ユーザ   システム       経過  
    11.711      6.689     18.254 

年齢の計算は，簡単なのです。試験年から生まれ年を引くだけ。
ただし，「生まれ月が試験月より後（大きい）なら」，
または，「生まれ月が試験月と同じで，生まれ日が試験日より後（大きい）なら」
1を引く。

わざわざ as.Date となっているけど，そこから年月日（スカラー y, m, d）を取り出す。
試験年月日は yyyy, mm, dd の整数で与える（as.Date から取りだしてもよいけど）

age2 <- function(x, yyyy = 2018, mm = 09, dd = 13) {
	y <- as.integer(substr(x, 1, 4))
	m <- as.integer(substr(x, 6, 7))
	d <- as.integer(substr(x, 9, 10))
	val <- yyyy - y
	if (m > mm || (m == mm && d > dd)) {
		val <- val-1
	}
	val
}
system.time({
for (i in 1:nrow(dat)) {
	dat$year2[i] <- age2(dat$birth[i])
}
})

でも，これでは for を含むし，年月日の取り出しコストが掛かるのでかえって遅い。

   ユーザ   システム       経過  
    13.662      6.667     20.312 

一応答えは合っているようだ。

> all(dat$year == dat$year2)
[1] TRUE

for を避けるためにベクトル演算を行うようにする。
関数には dat$birth ベクトルを渡す。
年月日の取り出し，それに基づく年齢（val）の計算も全部がベクトル演算。

age3 <- function(birth, yyyy = 2018, mm = 09, dd = 13) {
	y <- as.integer(substr(birth, 1, 4))
	m <- as.integer(substr(birth, 6, 7))
	d <- as.integer(substr(birth, 9, 10))
	val <- yyyy - y
	ifelse(m > mm | (m == mm & d > dd), val-1, val)
}
system.time(dat$year3 <- age3(dat$birth))

27 倍ほど速くなった。

   ユーザ   システム       経過  
     0.424      0.023      0.448 

答えも合っているようだ。

> all(dat$year == dat$year3)
[1] TRUE

なお，val を計算している 2 行を（ifelse を使うのが気持ち悪いからといって）

	yyyy - y - (m > mm | (m == mm & d > dd))

にしても，ほとんど速度に違いはない。

年月日が3つのベクトルで与えられれば as.Date から取り出す必要がないのでもっと速い。
birth を "2018-08-15" とわざわざ入力しているのか，2018, 8, 15 という3つの整数を入力しているのかによる。
後者なら，その3つの整数から文字列を作って as.Date していることになる。
3つの整数値を入力してそれをそのまま使うなら，人間もコンピュータも無駄なことをしなくてもよいということになる。

青木先生　様；

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

多面的にご教示いただき、
特に、for を避けるためにベクトル演算を行うようにするところは、
大変に勉強になります。

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

素朴な疑問ですが、ここのところ青木先生が代入演算子を

<-

ではなく、

=

と書かれるようになったのには、何か意味があるのでしょうか？

※書名は忘れましたが = ではなく <- の方がいいと書いてあった記憶があります。

最近になってソースコード表現の定石も変わってきているのでしょうか、、、

ソースコード表現の定石は変わっていないと思います
「GoogleのRスタイルガイドは、割り当てのために「=」を禁止する」ということもあり，
? "<-" で表示される R のマニュアルでも，The operator <- can be used anywhere, whereas the operator = is only allowed at the top level (e.g., in the complete expression typed at the command prompt) or as one of the subexpressions in a braced list of expressions. とありますので，<- を使うのが妥当であることに変わりありません。

= を使うようになったのは些細な理由からですが，= を使うようになっての感想は，何の不便もなく，タイプ数も（わずか1文字の違いですが）少なく，他の言語を使ってきた経験からも = の方が違和感がないなあというだけです。

公の場では <- を使うようにしたいと思います。

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

私はRが経験年数的には一番長いので　<-　に慣れてしまっていますが、多言語習得ユーザーとしては（明示的な問題が生じない限り）　=　でも不都合はなさそうと思います。

※代入演算子とは別に、最近のR使いの方はパイプを多用した記述をされるので、コードを読むのに一苦労することが多いです。

いつも勉強させていただいています。
Rの基本操作だと思いますが、わかりません。
今、次のようなデータが1000個あるとします。
x<- rnorm(1000, 15.5, 3.3)

このデータの95％区間にあるデータを取り出したいのですが、どのように取り出せばよいのでしょうか。
たとえば、
> y<- quantile(x, c(0.025, 0.975))
とすると、yのlength=2となり、950個のデータが取り出せません。
どなたか、yが950個のベクトルになるようなコマンドを教えてもらますでしょうか。

x <- rnorm(1000, 15.5, 3.3) のベクトル x から，10 以上， 15 以下のデータを取り出すときはどのようにしますか？

x[x >= 10 & x <= 15] ですね

では y <- quantile(x, c(0.025, 0.975)) から，下限と上限を取り出すのはどのようにしますか？

y[1], y[2] ですね

それをあわせれば

x[x >= y[1] & x <= y[2]] です

青木先生、
ありがとうございました。
ここ2日くらいそのような関数というかコマンドばかり探していました。
今後ともよろしくお願いします。

いつもお世話になっております。

基本的なことで恐縮なのですが、ご教示いただけますか。

手元に、テスト得点データがあります。実験群、統制群の2グループで、サンプルサイズはそれぞれ、30、25です。平均点の差の検定を行いたい場合は、普通はt検定を行うのではと考えます。

これを、サンプルサイズが小さいという理由で、ノンパラメトリック検定で行うことは大丈夫でしょうか。なぜこのような質問をしたかと言いますと、ある論文の査読をしているのですが、サンプルサイズが小さいとありますが、GPowerで必要な被験者数を求めると上のサイズは小さいとは言えず、テスト得点ですので、母分布は正規分布をなしていると考えます。また、本論文では、正規分布をなしていない根拠として、手元のデータに対して正規分布をなしているかどうかの検定を行い、正規分布をなしていないとも言っています。

そもそもとしまして、ノンパラにするかどうかは手元のデータが正規分布をなすかどうかで決めるものなのか、以前、この掲示板のどこかで指摘ぐあったような記憶があります。

どういう視点で、ノンパラか否かを決めるべきか、また、決定のプロセスの中でやってはいけないこと、ぜひご教示ください。

> サンプルサイズが小さいという理由で、ノンパラメトリック検定で行う

あまり説得力はないでしょう
しかし，

> GPowerで必要な被験者数を求めると上のサイズは小さいとは言えず

とはいっても，マン・ホイットニーの U 検定の検定力は，条件を満たしているときの平均値の差の検定（ｔ 検定）の 95% ほどもあるので，劣っている検定手法を採用したと批難するのはちょっと厳しすぎるかも。

> 手元のデータに対して正規分布をなしているかどうかの検定を行い、正規分布をなしていない

検定の多重性の問題が出てくるでしょう
しかし，理論的には正規分布に従うと考えられるものの，標本データでそれを確認できないとすれば，正規性の検定を行っても批難はされないでしょう

ノンパラメトリック検定は種々の仮定を置かないので，ノンパラメトリック検定で有意ならばややこしいことを避けることができるという長所があります。

> 決定のプロセスの中でやってはいけないこと

「データを見てから検定法を選択すること」はやってはいけないことと言えるでしょう。

青木先生

とても詳細なご説明をありがとうございます。勉強になりました。

今後ともよろしくお願いいたします。

失礼いたします。

青木先生　様；

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

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

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

半角、全角、大文字、小文字が混在する文字列を、全角大文字に統一したいのです。

例示ですが、

moji <- "ボーイング787 ボーイング７８７ boeing787 Boeing787 ｂｏｅｉｎｇ７８７ ＢＯＥＩＮＧ７８７"
⇒
"ボーイング７８７　ボーイング７８７　ＢＯＥＩＮＧ７８７　ＢＯＥＩＮＧ７８７　ＢＯＥＩＮＧ７８７　ＢＯＥＩＮＧ７８７"

大文字に置換は、toupper(moji)を使うえばよいことは分かりましたが、
誠に情けないのですが、半角→全角の関数を見つけることができません。

ご教示いただけましたら大変に助かります。
お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

文字の変換には基本的には chartr 関数を使えば良いです
例題の場合には，大文字変換に toupper を使えば，少し簡単に書けます

> moji = "ボーイング787 ボーイング７８７ boeing787 Boeing787 ｂｏｅｉｎｇ７８７ ＢＯＥＩＮＧ７８７"
> moji = toupper(moji)
> chartr("A-Z0-9", "Ａ-Ｚ０-９", moji)
[1] "ボーイング７８７ ボーイング７８７ ＢＯＥＩＮＧ７８７ ＢＯＥＩＮＧ７８７ ＢＯＥＩＮＧ７８７ ＢＯＥＩＮＧ７８７"

> moji = "ボーイング787 ボーイング７８７ boeing787 Boeing787 ｂｏｅｉｎｇ７８７ ＢＯＥＩＮＧ７８７"
> chartr("a-zA-Zａ-ｚ0-9", "Ａ-ＺＡ-ＺＡ-Ｚ０-９", moji)
[1] "ボーイング７８７ ボーイング７８７ ＢＯＥＩＮＧ７８７ ＢＯＥＩＮＧ７８７ ＢＯＥＩＮＧ７８７ ＢＯＥＩＮＧ７８７"

青木先生　様；

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

テキスト分析で、語のゆらぎを吸収する際に必要な前処理であり、大変に助かりました。

つまらない質問で、ご高名な青木先生に失礼にならないか、
小心者ですので、躊躇していましたが、
勇気をふりしぼってご相談して良かったと思います。

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

対応のある2標本の母分布に違いがあるかを検証したいのですが、
どのような分析手法を採用すべきでしょうか。

2標本コルモゴロフ・スミルノフ検定は独立2群の比較なので
この場合は使えないと考えております。

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

「母分布に違い」というのは，位置の母数の違いなのか，散布度の違いなのか？それとも，両方含めての違いなのか？

統計手法の選択ガイド
http://aoki2.si.gunma-u.ac.jp/FlowChart/Tutorial.html
をやってみるとどうなりました？

青木先生
コメントありがとうございます。

「母分布に違い」とは散布度（分布の形）の違いを意図しておりました。説明不足であり、
申し訳ありません。なお、データは離散型の比率尺度データです。

統計手法の選択ガイド
http://aoki2.si.gunma-u.ac.jp/FlowChart/Tutorial.html
では、「検定→その他」で適合度の検定のページ
http://aoki2.si.gunma-u.ac.jp/lecture/tekigoudo.html
にたどり着きますが、対応のある2群の散布度（分布の形）に違いがあるかを検証する分析が
見つけられず、質問させていただきました。

また、
http://aoki2.si.gunma-u.ac.jp/FlowChart/kentei_heikinti_taiouari.2.html
のページで紹介されている検定は平均値（代表値）の比較が目的であり、
散布度（分布の形）の違いを検証するものではないと思っており、判断に迷っております。

青木先生　様；

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

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

添付ファイルでお示しをします。　ご覧ください。

表1から表2を作成したいと思います。

表1のデータフレームは、以下のようになっています。
・socreは、すべてNA
・score1とscore2の組み合わせは、以下の３ケースのいずれか
　　・score1がNA,score2は数値
　　・score1が数値,score2はNA
　　・score1,score2両方ともにNA

score1、score2で、数値がセットされている側の値をscoreにセットする。

データフレームを dat として、Rコードを書きました。

jun <- which(! is.na(dat$score1))
dat$score[jun] <- dat$score1[jun]

jun <- which(! is.na(dat$score2))
dat$score[jun] <- dat$score2[jun]

score1、score2を、各々、NAの判定をして、scoreに代入しています。

これらを１つにまとめる書き方があれば、ご教示をお願いいたします。

お手数をおかけいたします。

いろいろやり方はあるでしょうが，もっとも真っ正直に書くならば

> dat = data.frame(
+   score=NA,
+   score1=c(NA, NA, 2.174, 6.042, NA, 2.656, NA, NA, NA, NA, 3.153, 4.344),
+   score2=c(1.481, 4.113, NA, NA, 2.674, NA, 2.500, 2.276, 4.162, NA, NA, NA))

> dat$score = ifelse(is.na(dat$score1) & !is.na(dat$score2), dat$score2, ifelse(!is.na(dat$score1) & is.na(dat$score2), dat$score1, NA))
> dat
   score score1 score2
1  1.481     NA  1.481
2  4.113     NA  4.113
3  2.174  2.174     NA
4  6.042  6.042     NA
5  2.674     NA  2.674
6  2.656  2.656     NA
7  2.500     NA  2.500
8  2.276     NA  2.276
9  4.162     NA  4.162
10    NA     NA     NA
11 3.153  3.153     NA
12 4.344  4.344     NA

トリッキーな書き方ならば，出現する可能性のある数値よりも小さい値を dummy にわりあてて，以下のようにすることも可能

> dummy = -99999 # 出現する可能性のある数値が正ならばこのような値を仮定する
> dat[is.na(dat)] = dummy
> dat$score = ifelse(dat$score1 > dat$score2, dat$score1, dat$score2)
> dat[dat == dummy] = NA
> dat
   score score1 score2
1  1.481     NA  1.481
2  4.113     NA  4.113
3  2.174  2.174     NA
4  6.042  6.042     NA
5  2.674     NA  2.674
6  2.656  2.656     NA
7  2.500     NA  2.500
8  2.276     NA  2.276
9  4.162     NA  4.162
10    NA     NA     NA
11 3.153  3.153     NA
12 4.344  4.344     NA

どちらが効率がよいか？

青木先生　様；

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

今回も、ご丁寧なご教示をいただき、大変に良い勉強をさせていただきました。

２つの方法とも理解できましたが、ifelse()関数を使わせていただきます。

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

サンプリング誤差について質問です。

サンプリング誤差について、視聴率での以下のサイトでの考え方や計算方法は分かったのですが、回答選択肢が3以上ある場合、例えば、非常に好き、やや好き、どちらでもない、やや嫌い、非常に嫌い、などでも、誤差について同様に考えて良いのでしょうか？

https://toukeigaku-jouhou.info/2015/08/25/audience-rate/


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

　サンプリング誤差って何かな？と思ったら、母集団の比率の95%信頼区間の求め方のことですか。
　選択肢がたくさんあっても同じだと、僕は思います．．？

鈴木康弘 先生

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

視聴率や知名度、購入経験のように、yes no いずれかから一つ選ぶ場合は二項分布であって、それが正規分布に近似することから、95％の信頼区間推定が出来るとテキストにあったのですが、

非常に好き、やや好き、普通、やや嫌い、非常に嫌いとかの、
多項分布の場合、誤差をどう考えたらよいかご教示下さい。

例えば、サンプル数600、非常に好き10％、やや好き15％、ふつう30％
やや嫌い20％、非常に嫌い25％の場合、
95％信頼区間で

非常に好き、非常に好き＋やや好きは

どうなるのでしょうか？

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

各カテゴリーは排他なので，
「非常に好き」と「それ以外」で，sqrt(0.1*0.9/600)
「非常に好き＋やや好き」と「それ以外」で，sqrt(0.25*0.75/600)
というように考えるということです。

青木教授、アドバイス、大変、ありがとうございました。

よろしくお願いします。

RでKuder RichardsonのKR20を計算するためのスクリプトは公開されているのでしょうか。教えていただけますでしょうか。

検索しましたか？

https://rdrr.io/github/cddesja/validateR/man/kr20.html
などもあります。

ここに書かれていることが難しければ，ともかく以下のようにすればよいでしょう。

devtools がまだインストールされていないなら，以下の 1 行を一度だけ

install.packages("devtools")

インストルされていた（インストールした）なら，以下の 2 行を一度だけ
library(devtools)
install_github("cddesja/validateR")
そして，kr20 を使うときには R を立ち上げるたびに，以下の 1 行を 一度だけ
library(validateR)
さて，準備が整ったら，データを以下のように準備する
例では 7 人の被検者に 6 項目の変数を 0/1 で準備
x = matrix(c(
  1,1,0,1,1,0,
  1,0,0,1,1,1,
  0,1,1,0,1,1,
  0,0,0,1,1,0,
  0,1,0,1,0,1,
  0,0,0,1,1,0,
  0,0,1,0,0,0), byrow=TRUE, ncol=6)
その後，
kr20(x)
とすれば，結果は
>  kr20(x)
[1] 0.1382488
となる。

＝＝＝＝

なお，やっていることは簡単なので，自分で以下のような関数を書けば，面倒なことは一切なしで結果が得られる。
計算方法につては Wikipedia にさえ書かれているが，例えば，
https://www.radford.edu/.../KR21%20&%20KR20%20examples.xls
には，Excel での計算例が書かれているので，計算方法を理解すれば以下のような関数を書くのは容易である（R がそこそこできればではあるが）
KR20 = function(x) {
  cm = colMeans(x)
  k = ncol(x)
  k/(k-1)*(1-sum(cm*(1-cm))/var(rowSums(x)))
}

> KR20(x)
[1] 0.1382488

青木先生

大変詳細な説明をいただきありがとうございます。

無事解決できました。

いつもお世話になっております。

2クラス（実験群と統制群）で、授業を行う前と行った後の平均値の差を見たいのですが、分散分析（被験者間×被験者内）で分析をしたところ、交互作用は非有意、主効果がそれぞれの要因で有意、となりました。

解釈に苦しむところがあり、t検定でpreの平均を調べたところ、差がない、という結果となりました（等分散性の検定もOK）。

この結果を利用して、事前の能力が同じだと仮定して、postのみをｔ検定で処理してもよろしいでしょうか？実際に行ったところ、有意でした。

アドバイスをお願いいたします。

　単に授業を行った後の2クラスの比較をしたいならｔ検定でいいんでしょうが、
2要因分散分析ではいけないんですか？

鈴木康弘先生

お返事、ありがとうございます。簡単に説明させていただきますと以下のようになります。

・（要因A）被験者間要因（実験群と統制群）×（要因B）被験者内要因（preとpost）の2要因分散分析（ともに2水準）
・preにはなく、postにおける平均値に差があることが実験仮説
・分散分析の結果、交互作用がなく、要因AおよびBそれぞれに主効果があった（グラフに書きますと、ともに右肩上がりの平行線）

これを、実験仮説に即してうまく解釈できず...つまり、要因Aの主効果に差があると言うことは、「実験群の方が統制群よりも平均が高い」ということですよね？これだと「postで差が出ている」と解釈できないのではないかと考えた次第です。なので、ｔ検定で行ったところ、うまい具合に出てくれたので、上のような質問をさせていただいた次第です。

もしよろしければ再度コメントをいただけますと幸いに存じます。

　授業前に有意差はありませんでした、授業後には有意差はありました、と言うことはできるでしょう。

　でも、分散分析での交互作用はどうなのか？と問われたら、ありませんでした、と白状するしかないのでは。

青木先生
お世話になっております。

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

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

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

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

フィッシャーの正確検定を使うなら，ボンフェローニ法を適用するかな？

カイ二乗検定を使うならば
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

青木先生

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

取り急ぎお礼まで

予測区間の式について教えて頂けましたら幸いです。

検索したところ、下記サイトがあり、式がすべて違い、値も違います。
どれを使えばよいのかわかりません。
コメント頂けましたら幸いです。
https://qiita.com/ksksk/items/75ba95337ccdb32e7cb1
https://ameblo.jp/cabay/entry-12184838325.html
http://heartland.geocities.jp/ecodata222/ed/edj1-1-4.html
https://www.monodukuri.com/gihou/article/946

それはそれとして、今回やりたいことは、
重回帰分析をした後に予測区間をプロットすることです。
お勧めのやり方をご教示頂けましたら助かります。
Ｒは少し使えます。

お手数おかけいたしますがよろしくお願いいたします。

2 番目のページと 3 番目のページは，標本平均値の信頼区間と予測区間の話です
1 番目のページと 4 番目のページは，単回帰分析における信頼区間と予測区間の話です

1 番目と 4 番目は，基本的に同じですが，

1 番目のページの式で t_n-p-1((1_0.95)/2) の部分は，(1_0.95)/2 は，下に「自由度 n-p-1 の t 分布の (1+0.95)/2(1+0.95)/2 パーセンタイル」と書かれているので誤解を生じるのでしょう。1 番目のページの結果は正しいです。

4 番目のページの式は，t(n-2, α) を「α は有意水準で 5%，(1-α) が信頼水準で，ここへ 95% が入る。」と紛らわしいことを書いているが，信頼水準を 95% にするには t(n-2, α/2) でなければならない。つまり，上側確率が 0.025 である t 値を求める必要がある。αを使ってしまうと間違った答えが出る。

> 重回帰分析をした後に予測区間をプロットする

予測区間の求め方については，
http://sayalabo.com/labo4.html
に説明があります。

ただ，プロットは難しいでしょうね。

早速のご回答ありがとうございます。
自分の勉強不足を痛感いたしました。
昨日一日勉強いたしました。
まだ理解できていませんが、取り急ぎお礼とご報告をさせて頂きます。
理解が進んだ後、また質問させて頂くこともあると思います。
今後ともよろしくお願いいたします。

ロジスティック回帰分析を行う際、説明変数として用いたい変数が、３種類以上の値を含む名義尺度データであった場合に、それを複数のダミー変数に置き換える方法が、いろいろあることについて悩んでいます。

インターネットで調べられた限りですが、おもに次の(1)〜(3)の方法があるようでした。
これらはどのように使い分ければよいでしょうか。

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

説明変数Ａに、"x","y","z" の３種類の名義が含まれている場合、

(1)ダミー変数 { d1, d2, d3 } を作り、変数Ａが、

"x" なら { 1, 0, 0 }
"y" なら { 0, 1, 0 }
"z" なら { 0, 0, 1 }

　に置き換える方法。

(2)ダミー変数 { d1, d2 } を作り、変数Ａが、

"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { 0, 0 }

　に置き換える方法。

(3)ダミー変数 { d1, d2 } を作り、変数Ａが、

"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { -1, -1 }

に置き換える方法。

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

上に挙げたうち、どれか１つの方法を紹介しているサイトはいくつか見つかるのですが、それぞれの方法の比較や、どう使い分ければよいか、という点についての情報がなかなか見つからず、質問させていただきました。

また、この問題が体系的に解説されている書籍などで、もしおすすめのものがありましたら、ご紹介いただければ幸いです（それぞれの方法の結果の解釈、オッズ比への逆変換の方法、ステップワイズ法で変数を増減する場合にダミー変数はどう扱うのか、などについてもよく知りたいと思っています）。

少なくとも（１）は多重共線性でダメなのでは

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

多重共線性は、互いに相関が強い複数の変数が説明変数の中に含まれていると、起こる問題ですよね。
（１）の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…？


（追記）
すみません、次のように修正しました
修正後「 d1 が 1 なら d2 + d3 は必ず0に決定されてしまうので、」
修正前「 d1 が 1 なら d2 とd3 は必ず0に決定されてしまうので、」

> １）の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…？

そういうこと

R では，変数Ａが名義尺度であるということを
FA = factor(A)
としておけば，
glm(result ~ FA)
でちゃんと処理してくれます。
変数 A は 2 個のダミー変数で表現できます。
3 個以上の変数を使っても，R では 3 個目からは分析に使われません。他の統計ソフトではエラーメッセージを吐いて死んでしまうこともあるでしょう。
独立な2変数変数であれば，どんな値であっても(0 と 1 でなくても -6.23123 と 578.21234 でもなんでも），その変数に対する係数は見かけは違いますが，その係数を使って計算される予測値は全く同じになります（そうでないと困ります。よい結果が出るような変数値のセットがある！なんてことになると，客観性がなくなります）。
ではなぜ，（3）のような変数を作らない！！かというと，出てくる値が3種類になるからです。0と1だと計算も簡単だしわかりやすい。
それにしても 0,1,-1 を与えるというのは，どんなメリットがあるのかなあ？

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

(3)の意義については、この方法が載っていたpdf資料の内容を後でアップいたします。
（私のパソコンのDLフォルダの整理が悪く、見つかり次第）

通常は (2) の置き換えを行い、(1)はダメ、ということですね。

さっそく、R で、名義尺度の変数を factor型 にし、glm メソッド にかけてみました。
（スクリプトは、別に投稿します）

このスクリプトは次のようなロジスティック回帰分析を行うつもりで書きました。

・データセット: birthwt　低体重出生とそのリスク因子のデータ
・目的変数: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
・説明変数: 母親の年齢、母親の体重、母親の人種、...etc

母親の人種 race {1:白人 2:黒人 3:その他 } が名義尺度なので、スクリプトの中であらかじめ次のように変換しています。

# race を factor型 に変換する
DF1$race <- factor(DF1$race)

glmメソッド の summary をみると、以下に引用したように race2, race3 という２つの変数が生成されているのを確認しましたが、これらは race { 1:白人 2:黒人 3:その他 } をどのように変換した変数なのかを確認する方法がわかりません。

---------------------------------------------------------------
# 結果表示
summary(out1.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8946  -0.8212  -0.5316   0.9818   2.2125  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.48062    1.19689   0.402  0.68801   
age         -0.02955    0.03703  -0.798  0.42489   
lwt         -0.03397    0.01524  -2.229  0.02580 * 
race2        1.27226    0.52736   2.413  0.01584 * 
race3        0.88050    0.44078   1.998  0.04576 * 
smoke        0.93885    0.40215   2.335  0.01957 * 
ptl          0.54334    0.34540   1.573  0.11571   
ht           1.86330    0.69753   2.671  0.00756 **
ui           0.76765    0.45932   1.671  0.09467 . 
ftv          0.06530    0.17239   0.379  0.70484   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 201.28  on 179  degrees of freedom
AIC: 221.28

Number of Fisher Scoring iterations: 4
---------------------------------------------------------------

ここでの {race2, race3} は、たとえばですが、内部的に

race=1 なら { 1, 0 }
race=2 なら { 0, 1 }
race=3 なら { 0, 0 }

のように変換されているのだと思いますが、それをどのように確認できるのでしょうか？　

str(out1.glm)

としてみたりもしましたが…

以下、スクリプトです。

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

###################################################
#
# ロジスティック回帰分析例
#
# https://www.gixo.jp/blog/3200/
# より一部改変
#

library(MASS)

str(birthwt)

# *** 目的変数 ***
# low:新生児の体重が2.5kg未満か否か(2.5kg未満が1)
#
# *** 説明変数 ***
# age:母親の年齢
# lwt:母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke:母親の喫煙有無(喫煙有りが1)
# ptl:母親の早産経験の有無(経験有りが1)
# ht:母親の高血圧の有無(有りが1)
# ui:母親の子宮神経過敏の有(有りが1)
# ftv:妊娠第１期における医師の訪問回数
# 
# *** 使わない ***
# bwt:新生児の体重(グラム単位)

# ■ 前処理 

# bwtは使わないので除く
DF1 <- birthwt[,1:9] 

# lwt をポンド単位からグラム単位に変換する
DF1$lwt <- DF1$lwt * 0.454 

# race を factor型 に変換する
DF1$race <- factor(DF1$race) 

str(DF1)

# ■ ロジスティック回帰

out1.glm <- glm(
  low ~ . ,
  family = binomial(link = "logit"),
  data = DF1)

# 注: モデル式 low ~ . は、目的変数に low、説明変数に low 以外の全ての変数を使うという意味 

# 結果表示
summary(out1.glm)

# 各説明変数の係数を指数変換し、オッズ比を表示
exp(out1.glm$coefficients)

例としてあげられたデータ birthwt$race は整数型で 1, 2, 3 をとる変数でしたが，factor(birthwt$race) で factor型の変数にすると，1, 2, 3 というレベルを持つ factor 変数になります。

> DF1$race
  [1] 2 3 1 1 1 3 1 3 1 1 3 3 3 3 1 1 2 1 3 1 3 1 1 3 3 1 1 1 2 2 2 1
  途中省略
[161] 1 3 3 3 2 1 3 3 1 1 2 2 2 3 3 1 1 1 1 2 3 3 1 3 1 3 3 2 1
Levels: 1 2 3　← ここ

このような factor 変数を多変量解析などに使用すると，一番最初の Level が基準値 0 になります。つまり，必要なダミー変数において，取る値の全てが0
一般的には，変数を factor 変数に変換すると，"辞書順の最初の値が基準値 0 になります。
詳しくは factor 関数のオンラインヘルプを参照してもらいたいですが，ちょっと見以下のような結果になりますが，訳ワカメに陥らないようにしっかり理解してください。

> sex = c("male", "female", "U.K.", "male", "male", "female")
> factor(sex, labels=c(1,2,0))
[1] 2 1 0 2 2 1
Levels: 1 2 0
> factor(sex, levels=c("male", "female", "U.K."))
[1] male   female U.K.   male   male   female
Levels: male female U.K.
> factor(sex, levels=c("male", "female", "U.K."), labels=c(1,2,0))
[1] 1 2 0 1 1 2
Levels: 1 2 0

さて，このようなことから，R が作っているダミー変数を自分で作るとすると以下のようになります

> DF2 = DF1
> DF2$race = NULL # factor である race を削除
> Race = birthwt$race # 元のデータは整数型 1, 2, 3
> DF2$Race2 = (Race == 2)+0 # Race が 2 のときのみ整数値 1，それ以外なら 0
> DF2$Race3 = (Race == 3)+0 # Race が 3 のときのみ整数値 1，それ以外なら 0
> head(DF2) # Race = 1 のときは DF2$Race2 = DF2$Race3 = 0
   low age    lwt smoke ptl ht ui ftv Race2 Race3
85   0  19 82.628     0   0  0  1   0     1     0
86   0  33 70.370     0   0  0  0   3     0     1
87   0  20 47.670     1   0  0  0   1     0     0
88   0  21 49.032     1   0  0  1   2     0     0
89   0  18 48.578     1   0  0  1   0     0     0
91   0  21 56.296     0   0  0  0   0     0     1
> out2.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF2)
> summary(out2.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8946  -0.8212  -0.5316   0.9818   2.2125  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  0.48062    1.19689   0.402  0.68801
age         -0.02955    0.03703  -0.798  0.42489
lwt         -0.03397    0.01524  -2.229  0.02580
smoke        0.93885    0.40215   2.335  0.01957
ptl          0.54334    0.34540   1.573  0.11571
ht           1.86330    0.69753   2.671  0.00756
ui           0.76765    0.45932   1.671  0.09467
ftv          0.06530    0.17239   0.379  0.70484
Race2        1.27226    0.52736   2.413  0.01584
Race3        0.88050    0.44078   1.998  0.04576

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 234.67  on 188  degrees of freedom
Residual deviance: 201.28  on 179  degrees of freedom
AIC: 221.28

Number of Fisher Scoring iterations: 4

ということで，DF1 を使ったときと同じ結果になりました。

なお，reorder という関数がありますが，この存在を知ると，きっと悪夢を見ます（^_^;）

複雑怪奇な変換規則を理解するよりは，自分で「適切な」ダミー変数を作成する方が幸せかも知れません。

factor関数をテストしてみると、自動的に定まる Levels の順序が独特ですね…。

> factor(c(100,1,3,1,-5))
[1] 100 1 3 1 -5
Levels: -5 1 3 100 → 数値は昇順なのはまあ分かりますが、

> factor(c("b","a","A","c","A","b","C"))
[1] b a A c A b C
Levels: a A b c C → 文字列は辞書順だが小文字が先…？

最初のLevelがどの値になるか、いまいち確信が持てないので、ご説明にあるように、名義尺度をfactor型に変換する際は、 Levels パラメータ と Labels パラメータをちゃんと与えて、最初のLevelになる値を明示すればいい、と考え、次のように進めました。

> DF3 = DF1
> DF3$race = NULL
> DF3$race = factor(birthwt$race,levels=c(1,2,3),labels=c("白人","黒人","その他"))
> out3.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF3)

> summary(out3.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF3)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4806 1.1969 0.40 0.6880
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race黒人 1.2723 0.5274 2.41 0.0158 *
raceその他 0.8805 0.4408 2.00 0.0458 *

この結果が、お示しいただいた out.glm2 の結果と一致するのを見て、やっと気付いたのですが、out1.glm に含まれるダミー変数 race2 , race3 の "2"や"3"は、2:黒人、3:その他 のラベルなんですね。
（No.22543を投稿した時点では、通し番号のようなものかな、と誤解していました）

ダミー変数 { race黒人, raceその他 } は、おそらく

race=1 について { 0, 0 }
race=2 について { 1, 0 }
race=3 について { 0, 1 }

という対応になっていて、

> exp(out3.glm$coefficients)
(Intercept) age lwt smoke ptl
1.62 0.97 0.97 2.56 1.72
ht ui ftv race黒人 raceその他
6.44 2.15 1.07 3.57 2.41

から、

・人種が「黒人」であることは、白人である場合と較べ、低体重児の出生リスクが 3.57倍になる
・人種が「その他」　　〃　　2.41倍になる

という解釈に行き着くのか、と。

投稿No. 22537に挙げた（３）の方法
（ダミー変数を -1, 0, 1 などとおく方法）

についてですが、掲載されているのは、ある大学病院で開催された「ＥＢＭ・統計セミナー」という研修会のPDF資料でした（現在はDLできなくなっていました）。
この資料には、次の２つの変換方法が紹介されており、

ア）基準とするカテゴリ（reference group）を仮定する場合
イ）カテゴリ全体の平均をreferenceとする場合

-------------------
ア）基準とするカテゴリ（reference group）を仮定する場合

例　名義尺度データが｛ "薬剤なし", "薬剤Ｐ", "薬剤Ｑ" } のような場合

ダミー変数 D1, D2 を

　"薬剤なし"　=　{0,0}
　"薬剤Ｐ"　　=　{1,0}
　"薬剤Ｑ"　　=　{0,1}

とおく。ロジスティック回帰の結果、D1, D2 について、

　b1 : D1 に対応する回帰係数
　b2 : D2 に対応する回帰係数

が推定され、

　exp(b1) は"薬剤なし"に対する"薬剤Ｐ"のオッズ比、
　exp(b2) は"薬剤なし"に対する"薬剤Ｑ"のオッズ比、

という解釈となる。

--------------------
イ）カテゴリ全体の平均をreferenceとする場合

例　名義尺度データが｛ "地域Ａ", "地域Ｂ", "地域Ｃ" ｝のような場合

ダミー変数 D1, D2 を

　"地域Ａ"　=　{1,0}
　"地域Ｂ"　=　{0,1}
　"地域Ｃ"　=　{-1,-1}

とおく。
　
　b1: "地域Ａ"の回帰係数
　b2: "地域Ｂ"の回帰係数
　b3: "地域Ｃ"の回帰係数

とし、モデル上で、

　b1 + b2 + b3 = 0

という制約を置いたうえで推定する。
計算上は、b1, b2 しか推定されないが、b3 = - b1 - b2 として計算可能である。
exp(b1)をとればオッズ比が得られる。

--------------------
引用ここまでです。

最後の一文にある、「exp(b1)をとれば」は、「exp(-b1-b2)をとれば」ではないのかな…？ という疑問がありますが、原文ママです。

なお、この資料の末尾の参考文献には次が挙げられていました。

青木繁伸（2009）「Ｒによる統計解析」
大橋靖雄（2011）「医師のための臨床統計学基礎編」
丹後俊郎（1996）「路地すてぃく回帰分析SASを利用した統計解析の実際」
豊田秀樹（1992）「原因をさぐる統計学」
服部環・海保博之（1996）「Ｑ＆Ａ心理データ解析」
森實敏夫（2004）「入門　医療統計学」

No.22545でお示しいただいた、自分でダミー変数を作る方法を真似して、（３）の方法でダミー変数を作成する処理を追加してみました。

###################################################
#
# ロジスティック回帰分析例 2
# https://www.gixo.jp/blog/3200/
# より一部改変

library(MASS)
str(birthwt)
#
# low: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
# age: 母親の年齢
# lwt: 母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke: 母親の喫煙有無(喫煙有りが1)
# ptl: 母親の早産経験の有無(経験有りが1)
# ht: 母親の高血圧の有無(有りが1)
# ui: 母親の子宮神経過敏の有(有りが1)
# ftv: 妊娠第１期における医師の訪問回数
# bwt: 新生児の体重(グラム単位)

options(digits = 4)

### DF1 ダミー変数を、Rで自動生成する

# bwtは使わないので除く
DF1 <- birthwt[,1:9]

# lwtをポンド単位からg単位に変換する
DF1$lwt <- DF1$lwt * 0.454

# raceをfactor型に変換する
DF1$race <- factor(DF1$race,levels = c(1,2,3),labels = c("白人","黒人","その他"))
str(DF1)

# ロジスティック回帰
out1.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF1)
summary(out1.glm)
exp(out1.glm$coefficients)

#### DF4 ダミー変数を、白人{1,0} 黒人{0,1} その他{-1,-1} とする方法で自作する

DF4 <- DF1

# ダミー変数に変換する
race = birthwt$race
DF4$race白人 = ifelse(race==1, 1, ifelse(race==2, 0, -1))
DF4$race黒人 = ifelse(race==1, 0, ifelse(race==2, 1, -1))
DF4$race = NULL
str(DF4)

# ロジスティック回帰
out4.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF4)
summary(out4.glm)
exp(out4.glm$coefficients)

# "raceその他" の偏回帰係数、オッズ比
b3 = - coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"]
names(b3) = c("raceその他_coef")

b3_exp = exp(- coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"])
names(b3_exp) = c("raceその他_exp(coef)")

-----------------
スクリプトここまでです。
以下、結果を貼ります。
-----------------

> summary(out4.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF4)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1982 1.1509 1.04 0.2978
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race白人 -0.7176 0.2699 -2.66 0.0079 **
race黒人 0.5547 0.3232 1.72 0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.3

Number of Fisher Scoring iterations: 4

> exp(out4.glm$coefficients)
(Intercept) age lwt smoke ptl ht ui ftv race白人
3.3142 0.9709 0.9666 2.5570 1.7217 6.4450 2.1547 1.0675 0.4879
race黒人
1.7414

> b3
raceその他_coef
0.1629
> b3_exp
raceその他_exp(coef)
1.177

------------
結果は以上です。
資料の説明からすると、

「低体重児の出生確率は、全人種の平均出生確率を基準にすると、
　・白人の場合、0.4879倍
　・黒人の場合、1.7414倍
　・その他の人種の場合、1.177倍」

というような解釈になるのかな、と思いましたが…。
これ以上の詳細は、資料を作った方に直接質問するのが筋かと思いますが、具体的にやってみたことの報告のつもりで投稿しました。

すいません、教えて下さい。確率変数X,Yがあった時に、条件付き確率　P(X|Y)を考えます。P(X=2|Y=1)とすれば、何らかの確率が定数として求まると思うのですが、P(X=2|Y)というものを考えると、Yは確率変数なので、このP(X=2|Y)も確率変数になるのでしょうか？

ベイズの定理に関わる式 P(X|Y)は、Yがある値に確定した後での確率P(X)ですから、先にXが確定するなら、YはPに関与しないと考えるべきでしょう。従って、P(X=2)と等価の確定値でしょう。

どうもありがとうございました。
P(X=2|Y=1)もP(X=2|Y=2)も値としては同じになるのですかね。というか、そもそも、この「Y=」の部分は、関与しないわけですから、記法としての意味がないということなのでしょうか。

あくまで、Yが事前に確定してこその、P(X|Y)です。おっしゃる通り、Xが先に確定するなら、記法として成り立ちません。

よろしくお願いします。

指導法の効果を見るために，実験群と統制群のクラスで，それぞれテストを３回行いました。

初回のテストは実験前に，２回目は授業終了直後，３回目は遅延テストとして行いました。

データ分析の段階で，初回のテストの等分散性が確保できませんでした。この場合，２要因分散分析ではなく，共分散分析をしなければならないところまではなんとなく理解しています。

【質問１】分散分析の場合は，すべてのデータを一括で処理できるのですが，今回のデータは，

（１）初回と２回目の共分散分析
（２）初回と３回目の共分散分析

と分析しなければいけないのでしょうか？

【質問２】２回目と３回目の差についても興味があるのですが，この場合はどういう方法で見ることが可能なのでしょうか？

よろしくお願いします。

　あまり自信はないのですが．．

　等分散性が確保出来ない場合に共分散分析をしなければならないって、本当でしょうか？
　そういう場合にはいい方法はないので、2要因分散分析でごまかすしかないような．．．

お世話になります。

ログランク検定のサンプルサイズを計算する際、以下の理論式（プログラム）が広く使われていると学習いたしました。理解を深めるために、この計算結果で得られる理論値（サンプルサイズ）、と検出力の関係をシミュレーションで再現したいと考えております。
例えば、以下の例ですと、受け入れR年、追跡期間T年とした場合、各群57例という結果が得られます。ここで、群ごとに57例分の生存時間を乱数で発生させて、1000回か10000回程度、検定を回したときに、検出力（有意になった回数）が実際に8割を超えることを確認したいと考えております。

生存時間は、ハザードをパラメータとした指数分布に従うことを仮定していることは理解しております。しかし、打ち切り症例が発生する確率と、打ち切り例の場合の打ち切りまでの時間に、どのような分布とパラメータを仮定すべきなのか分からず、煮詰まっております。この点について、ご教授頂けないでしょうか？

お手数をおかけして恐縮でございますが、ご返信頂けますと幸いです。また、質問に不備がございましたら、ご指摘頂けますと幸いです。

#  R      : 受け入れ期間
#  Tr     : 追跡期間
#  TIME   : TIME年生存率
#  A      : 100A(%)生存率(<1)
#  B      : 100B(%)生存率(<1)
#  ALPHA  : 有意水準(<1)
#  BETA   : 検出力（1-BETA）
#  S      : 例数比(Na/Nb)

  R <- 3
  Tr <- 5
  TIME <- 5
  A <- 0.2
  B <- 0.05
  ALPHA <- 0.05
  BETA <- 0.2
  S <- 1

	LA <- -log(A)/TIME
	LB <- -log(B)/TIME
	LBER <- (LA + LB)/2
	QA <- 1/(1 + S)
	QB <- S/(1 + S)
	ZA <- qnorm(1 - ALPHA/2)
	ZB <- qnorm(1 - BETA)
	P00 <- LBER * LBER
	P01 <- LA * LA
	P02 <- LB * LB
	P10 <- LBER^2/(1 -exp(-LBER*Tr))
	P11 <- LA^2/(1-exp(-LA*Tr))
	P12 <- LB^2/(1-exp(-LB*Tr))
	P20 <- LBER^3*Tr/(LBER*Tr-1+exp(-LBER*Tr))
	P21 <- LA^3*Tr/(LA*Tr-1+exp(-LA*Tr))
	P22 <- LB^3*Tr/(LB*Tr-1+exp(-LB*Tr))
	P30 <- LBER^2*(1-(exp(-LBER*(Tr-R))-exp(-LBER*Tr))/(LBER*R))^(-1)
	P31 <- LA^2*(1-(exp(-LA*(Tr-R))-exp(-LA*Tr))/(LA*R))^-1
	P32 <- LB^2*(1-(exp(-LB*(Tr-R))-exp(-LB*Tr))/(LB*R))^-1
	N <- ((ZA*sqrt(P00*(1/QB+1/QA)) + ZB*sqrt(P01/QA+P02/QB))/abs(LA-LB))^2/(1+S)
	M <- ((ZA*sqrt(P10*(1/QB+1/QA)) + ZB*sqrt(P11/QA+P12/QB))/abs(LA-LB))^2/(1+S)
	O <- ((ZA*sqrt(P30*(1/QB+1/QA)) + ZB*sqrt(P31/QA+P32/QB))/abs(LA-LB))^2/(1+S)
	L <- ((ZA*sqrt(P20*(1/QB+1/QA)) + ZB*sqrt(P21/QA+P22/QB))/abs(LA-LB))^2/(1+S)

	c(S=S, 受け入れ期間=R, 生存率=TIME, ハザードA=LA, ハザードB=LB, 
		ALPHA=ALPHA, ZA=ZA, BETA=BETA, ZB=ZB, 
		"R=0,T=Inf"=N, "R=0, T>0"=M, "0<R<T"=O, "O<R=T"=L)

参考：[新版]実用SAS生物統計ハンドブック[SAS®9.4/R3.2.0対応]（臨床評価研究会）

> 打ち切り症例が発生する確率と、打ち切り例の場合の打ち切りまでの時間に、どのような分布とパラメータを仮定すべきなのか

シミュレーションデータの各例は，以下の順で生成されるのでは？
(1) 受入期間内に一様分布で研究に参入する
(2) エンドポイントまでの期間は指数分布で決まる
(3) エンドポイント発生が追跡期間後になるなら「打切り例」，前なら「故障例」
(4) 「打切り例」なら，研究対象であったのは研究参入から追跡期間終了までの期間
(5) 「故障例」なら，研究対象であったのは研究参入からエンドポイントまでの期間
これで，「打ち切り症例が発生する確率と、打ち切り例の場合の打ち切りまでの時間に、どのような分布とパラメータを仮定すべきなのか」というのは無関係，意味がないことになる。

なお，あなたがいっている「打ち切り症例」が「脱落例（追跡不能例）」を含んでいるなら，そもそもそういうものは考える必要がない。例えば，2群の平均値の差のサンプルサイズの決定などでも，測定不能，測定不全，拒否などによる例数減は考えないでしょう？

受入開始　　　　　　　　　受入終了　　　　　　　　追跡終了
   +..........................+.......................+
          o--------------x     受入終了前に死亡
        o-------------------------------x     追跡終了前に死亡
                o-------------------------------------@=======x     打切り例: 追跡終了後に死亡
- と = の和が指数分布で決定されるエンドポイントまでの期間
- が研究対象であった期間

青木 繁伸 先生

ご返信頂きましてありがとうございます。

打ち切り例は、脱落例を含むという前提で質問申し上げておりました。
二群の平均値の検定の場合は、（補完などの特別な処理をしなければ）脱落例は解析対象に含められませんが、生存時間解析においては、そういった対象も情報として利用するようですので、
サンプルサイズの計算においても実験途中の脱落を考慮していると理解しておりました。

ご教授頂いた内容を踏まえ、Rコードを書いてみました。
しかし、各群57例の場合の検出力は70~75%程度にとどまり、上の計算結果と整合しませんでした。
どこか誤りがありますでしょうか。
恐れ入りますが、改めてご指導賜りたく、何卒よろしくお願い申し上げます。

#初期化
sign <- NULL
pval <- NULL

#ハザードを計算
B <- 0.05
A <- 0.2
LB <- -log(B)/(5*365.25)
LA <- -log(A)/(5*365.25)

n <- 57*2

accrualtime <- 3*365.25#登録期間
followuptime <- 5*365.25#観察期間※登録期間を含む

for(N in 1:1000){
  #初期化
  x1 <- NULL#群
  time <- NULL#イベントまでの時間
  timec <- NULL#打ち切りまでの時間
  randomized.at <- NULL
  dt <- NULL
  status <- NULL
  i <- 0
  for(cal in 1:accrualtime){#登録期間の初日を1として、
    entry <- rbinom(size=1,n=1,prob=n/accrualtime)#当該時点に登録があるか（登録期間を通じて一様に集積されることを仮定する
    if(entry==1){                                 #当該時点に登録ありの場合
      randomized.at <- c(randomized.at,cal)       #登録日
      i <- i+1                                    #例数カウンター
      x1[i] <- rbinom(size=1,n=1,prob=0.5)        #ランダム化 
      
      if(x1[i]==0){time[i] <- rexp(n=1,rate=LB)}#対照群の場合の生存時間 
      else        {time[i] <- rexp(n=1,rate=LA)}#治療群の場合の生存時間
      
      timec <- followuptime-randomized.at[i]#試験終了による打ち切り
      
      if(time[i] < timec){        #打ち切り前にイベント発生
        status[i] <- 1            #イベント
      }                          
      
      if(timec < time[i]){        #イベント発生前に打ち切り
        status[i] <- 0            #打ち切り
        time[i] <- timec          #打ち切り症例には打ち切りまでの時間を代入置換
      }
      
    } 
    
    else {}#登録なしの場合はレコード発生させない
    
    if (i==n) break#予定症例数に達したら終了
    
  }

  #ログランク検定
  result <- survdiff(Surv(time)~x1)
  #p値を直接取り出せないので、検定統計量を格納し、そこからp値を導出
  pval[N] <- pchisq(result$chisq, length(result$n)-1, lower.tail = FALSE)
  if (pval[N] < 0.05){sign[N] <- 1} else{sign[N] <- 0}
}

sum(sign)/N
hist(pval)
> sum(sign)/N
[1] 0.724

一番の間違いは，
survdiff(Surv(time) ~ x1)
としたところ。これは，
survdiff(Surv(time, status) ~ x1)
でしょう。
もう一つは，
x1[i] <- rbinom(size=1,n=1,prob=0.5) #ランダム化
ですが，これでは群ごとのサンプルサイズが同じになりません。
その他も，無駄を省き，分かりやすく書いて，ベクトル演算による高速化を図ったので，参考にしてみてください。

f = function(n.each, B, A, accrualtime, followuptime, REPETITION = 1000) {
  library(survival)
  LB = -log(B) / followuptime
  LA = -log(A) / followuptime
  n = n.each * 2
  x1 = rep(0:1, each = n.each)
  pval = numeric(REPETITION)
  for (N in 1:REPETITION) {
    # 研究に登録される時刻（研究開始時刻を 0 とする）
    entry = runif(n, min = 0, max = accrualtime) # 一様分布に従う
    # 追跡可能如何に関わらず，登録時刻からエンドポイントまでの時間
    # 前半に B 群，後半に A 群をまとめて生成しても一般性を損なわない
    surv.time = c(rexp(n.each, rate = LB), rexp(n.each, rate = LA))
    # 研究開始時刻（0）からエンドポイントまでの時間は entry + surv.time
    # その時間が観察期間より短い場合は故障例（status = 1)
    status = (entry + surv.time) < followuptime
    # 故障例ならエンドポイントまでの時間は surv.time そのもの
    # 打切り例 status = 0 なら観察期間終了までの時間
    time = ifelse(status, surv.time, followuptime - entry)
    # ログランク検定
    # result = survdiff(Surv(time) ~ x1) # ここが間違い
    result = survdiff(Surv(time, status) ~ x1) # こちらが正しい
    # p値を直接取り出せないので、検定統計量を格納し、そこからp値を導出
    pval[N] = pchisq(result$chisq, 1, lower.tail = FALSE)
  }
  mean(pval < 0.05)
}

n.each = 57
B = 0.05
A = 0.2
# 以下の2項目の単位は何でもよい
accrualtime = 3 * 365.25 # 登録期間
followuptime = 5 * 365.25 # 観察期間 ※登録期間を含む
REPETITION = 1000 # シミュレーション試行回数

> f(n.each, B, A, accrualtime, followuptime, REPETITION) # 実行結果は毎回異なる
[1] 0.823
> f(n.each, B, A, 3, 5, 10000)
[1] 0.8253

青木 繁伸 先生

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

>一番の間違いは，
>survdiff(Surv(time) ~ x1)
>としたところ。これは，
>survdiff(Surv(time, status) ~ x1)
>でしょう。
基本的なところを間違っておりました。大変失礼いたしました。
この点を修正したところ、理論式と概ね整合のとれる結果が得られました。

>もう一つは，
>x1[i] <- rbinom(size=1,n=1,prob=0.5) #ランダム化
>ですが，これでは群ごとのサンプルサイズが同じになりません。
乱数なので同じサンプルサイズになるとは限らない、といういことですね。ご指摘ありがとうございます。

シミュレーション回数10000回で、頂いたプログラムを何度か回してみましたが、検出力は全て82%と少し、という感じでございました。
80~80.5%程度の範囲に収まるのかなと勝手に想像しておりましたが、これはやはり、シミュレーションと理論式の違い（誤差?）と考えてよいのでしょうか？

他のサンプルサイズ決定については，
http://www.st.nanzan-u.ac.jp/info/ma-thesis/2008/TANAKA/m07mm028.pdf
が，参考になるのではないかと思います。
私はとても追い切れませんが，必要とするあなたらならそれぞれを試してみるとよいのではないでしょうか？
その他にも，いくつかあるかも知れませんね（追求する情熱が私にはない）
何かわかったら，ここで教えてください。

青木先生　様；

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

青木先生からご教示いただいたノウハウを、読み返して整理しています。

[22299] Ｒ　ftable()関数 　　　　　投稿者：明石 投稿日：2017/03/11(Sat) 19:03

タイタニックのデータを題材にして、多重クロス表をftable()関数で整形出力する方法について、ご教示を頂戴しました。
大変に重宝しています。
ありがとうございました。

ftable(dat, row.vars=c("Sex", "Age", "Class"), col.vars="Alive")

                   Alive 死亡 生還
Sex  Age    Class                 
女性 子ども 1等             0    1
            2等             0   13
            3等            17   14
            乗務員          0    0
     大人   1等             4  140
            2等            13   80
            3等            89   76
            乗務員          3   20
男性 子ども 1等             0    5
            2等             0   11
            3等            35   13
            乗務員          0    0
     大人   1等           118   57
            2等           154   14
            3等           387   75
            乗務員        670  192

上記コンソール出力と同じようなレイアウトで、csv出力したいと思って調べましたが、分かりませんでした。

出力したcsvファイルをEXCELで開くとコンソールへ出力と同じようなレイアウトになることを期待してます。

ご教示を頂戴できれば、大変に助かります。
毎々、お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

私は，Excel がないので，LibreOffice で確認したのですが，

a = ftable(dat, row.vars=c("Sex", "Age", "Class"), col.vars="Alive")
write.table(format(a, quote=FALSE), "temp.csv", sep=",", row.names=FALSE, col.names=FALSE)

で良いと思います。

青木先生　様；

お忙しいところを失礼いたします、明石と申します。
今回も助けていただきました。
ありがとうございました。

出力したcsvファイルをEXCELで開いて、所望したフォーマットであることを確認しました。

本体の機能にちゃんと備わっているのですね。
これを機に、関数の引数（ふだんはデフォルト使用で意識しません）をヘルプで調べるようにします。

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

青木先生　様；

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

約2年前、青木先生から以下の有難いご教示を頂戴し、大変に重宝しています。

Re: Ｒ　str()関数の結果の取り出し、整形出力
　　投稿者：青木繁伸 2016/05/14(Sat) 05:12 No. 22001

a = data.frame(names(dat), sapply(dat, class))
colnames(a) = c("名前", "型")
print(a, row.names=FALSE)

などはどうでしょうか

ーーーーーーーーーーーーーーーーーー
名前、型の２列に、欠損値NAの個数 も追加した３列の出力 を検討しています。

まず、欠損値NAの個数 KAZU を計算します。
kazu <- numeric(ncol(dat))
for (i in 1:ncol(dat)) {
kazu[i] <-sum(is.na(dat[,i]))
}

次に、３列からなるデータフレームを作成します。
a = data.frame(names(dat), sapply(dat, class), kazu)
colnames(a) = c("名前", "型", "NA数")
print(a, row.names=FALSE)

data.frame()の中に、
欠損値NAの個数を算出する式を直接に書きたいと思っていますが、力不足で書くことができません。

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

どうしても data.frame 関数中に書きたいならば，

a = data.frame(names(dat), sapply(dat, class), colSums(is.na(dat))
colnames(a) = c("名前", "型", "NA数")
print(a, row.names=FALSE)

青木先生　様；

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

お示しをしてくださいましたＲコードは、よく理解できました。

どうしても data.frame 関数中に書きたいと思いますので、
勉強になりました。
ありがとうございました。

御世話になります。

3人が行った視覚的判定のバラツキについて、統計学的な検定方法を御教示いただけましたらと思い、投稿いたしました。

サンプルを30回自由に粉砕した後、その粉砕片を肉眼的に確認して評価しました(スコア：0-9）。3人がスコアを付けし、1000パターンのサンプルに対して行いました。なお粉砕片はその後溶媒に溶かし、物質の溶出濃度を測るため、実際の粉砕程度を後で確認しています。
スコアを付ける際は、スコア0-9の粉砕片の見本写真があり、その見本写真と実際の粉砕片を見比べてスコアを付けるため、バラツキは小さいということを証明したいのですが、それを証明する方法が分かりません。

級内相関とCV値を使用して、バラツキは小さかったと結論づけしましたが、スコアの正確性と再現性を証明するために十分なのか自信がありません。

恐れ入りますが御教示いただけましたら幸いです。よろしくお願い申し上げます．

クラスカル-ウォリス検定を考えております。しかし、1群はサンプル数が非常に小さい(20くらい)ことが予想され、もう1群も60前後と予想しています。このような場合、3群目は、200人など多いほうが良いのでしょうか。それとも、60前後など近い人数が良いのでしょうか。お教えいただければありがたいです。よろしくお願いいたします。

実験デザインの段階では，各群のサンプルサイズは成るべく同じにするようにする
一般の統計調査の場合には，層別の母集団サイズが異なる場合などは特に，サンプルサイズを揃えることが難しい或いは意味がないこともある
標本が得られたあとで群を分けて違いを検討する場合には，サンプルサイズは揃えようがない

結論として，群ごとのサンプルサイズが違っても，そのまま検定する
間違っても取ったデータを捨てて，サンプルサイズを揃えたりしないように

青木先生、
ご回答をありがとうございました。
確かに母集団のサイズが異なることを考慮しなければなりませんね。その点は、失念いたしておりました。
今回は第１群の母集団が小さく、サンプル数を増やすのが難しいと思われますので、第２群と第３群は揃えるように設計して見ます。
ご教示どうもありがとうございました。

お世話になっております。

２つのグループの、ある会議への参加回数の差に有意な差があるかどうかを確認する方法はt検定でよろしいでしょうか？参加回数がパラメトリックかノンパラメトリックなのかで悩んでおりまして。

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

参加回数は0を含むし，負の値を取らないので，比尺度でないことは確実です
参加回数4回の人と参加回数2の人と，参加回数10の人と参加回数8の人は，それぞれ2回の差があるということです。したがって，このデータは間隔尺度であるということがわかりますね。
パラメトリックかノンパラメトリックかといわれれば，パラメトリックでも問題ないでしょう。参加するかしないかというパラメータがあるなら，二項分布に従うのかも知れないし，もうすこし確率が低ければポアソン分布に従っているのかも知れない。そのいずれでもないかも知れないけど。
分布型，母数（パラメータ）が特定できない場合には，ノンパラメトリックと仮定して進めればよいのではないでしょうか？

青木先生

お忙しい中、詳細なご説明、ありがとうございました。今回は参加する（しない）のパラメータがありませんので、間隔尺度として処理したいと思います。

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

お世話になります。相変わらずつまらない質問で申し訳ございませんが、

m1 <- mean(d[class==1, 2])

の値をグラフで表示すると小数第10位くらいまで表示され格好悪いです。例えば、3桁にする場合の方法についてなのですが、私は

m1 <- mean(d[class==1, 2])
m1 <- round(m1, digits=3)

と書き、下のm1が使われるようにしましたが、meanの中で引数として入れることは不可能なのでしょうか？もしくは、もっとスマートな書き方はないのでしょうか？

平均のオブジェクトがm1からm7まであり、コピペして数字の修正をすれば問題ないのですが、洗練された書き方がないのかなと悩んでおります。

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

mean2 = function(x, digits=3) round(mean(x), digits)
として，mean2(d[class==1, 2]) などとするだけでしょう。

> set.seed(123)
> mean2 = function(x, digits=3) round(mean(x), digits)
> mean2(rnorm(100), digits=4)
[1] 0.0904
> mean2(rnorm(100), digits=6)
[1] -0.107547
> mean2(rnorm(100), digits=8)
[1] 0.1204651

あるいは必要に応じて options(digits=xxx) を使う

> options(digits=3)
> pi
[1] 3.14
> mean(rnorm(100))
[1] 0.106
> options(digits=5)
> pi
[1] 3.1416
 mean(rnorm(100))
[1] -0.0423

青木先生

お礼が遅くなり申し訳ございません。関数を自分で作成するんですね！

無事、きれいなコードができました。ありがとうございました。

はじめて投稿させて頂きます。
宜しくお願い致します。

現在、A, B, Cの3群の測定値について平均値を比較したいと考えております。

対立仮説としては
A>B>Cだと言う仮説を考えてます。
それぞれの群の具体的な平均値、SDも過去の文献で調べてあります

この様な状況では、測定値が正規分布に従う場合、
統計解析するには、以下の流れで進めると思います
(1) 分散分析
(2)その後事後検定(Scheffe等)

では、上記の仮設で試験を行う場合、
この試験に必要なサンプルサイズはどのように見積もればよいでしょうか？

(2)の事後検定用の例数設計でしょうか？
(1)の分散分析用の例数設計でしょうか？
それとも(1), (2)の両方の例数設計を行い、両方で十分な検出力を持つ例数にすべきでしょうか？

ご教示いただければ幸いです。

そもそも，例数設計は，結論の妥当性を保証するためのものですから，何を結論とするかがはっきりしていれば自明のことではないですか？
事後検定の結果も含めるなら，事後検定にも十分なサンプルサイズが必要でしょう。
事後検定のサンプルサイズの方が大きくなるかな。
いずれにしろ，どちらも十分なサンプルサイズを取るには，両方の例数設計をして，多い方を採らないと，小さい方を採ってしまうと他方の結論の妥当性に疑義が生じるでしょう？

ご返信ありがとうございます。
「それとも(1), (2)の両方の例数設計を行い、両方で十分な検出力を持つ例数にすべき」ということで承知しました。

少し話はそれますが、もう一点相談させて下さい
分散分析の例数設計の方法があるのは存じており、
寄与率、群数、有意水準、検出力を指定すれば、
例数設計できるのも存じてますが、

仮説検定で、分散分析の例数設計を実施して、
Scheffe等の事後検定の例数設計は実施しないという状況は
現実的にありえるのでしょうか？

G*Powerという検出力計算や例数設計のソフトウェアがありまして、
それを使うと分散分析の例数設計はできるのですが、
Scheffe等の例数設計はできないので、ご相談させて頂きました次第です。
http://www.gpower.hhu.de/

後学のためにご教示いただければ幸いです

例数設計法が開発されていない，あるいは，存在するけど特定のソフトには実装されていないということもあるでしょう。

> A>B>Cだと言う仮説を考えてます

　例数設計とは無関係なコメントになりますが、これは片側検定ですので、分散分析の適用は適切ではないと思います。

　有意水準を調整して、対立仮説がＡ＞Ｃ、Ａ＞Ｂ、Ｂ＞Ｃの片側ｔ検定を実施してはいかがでしょうか？

分散分析というのが広義の名称とすれば，Jonckheere 検定（ヨンキー検定）であれば片側検定も行えますね。ただ，Jonckheere 検定の例数設計理論が存在するかどうかは，よほど難しくなりますかね？

　分散分析もヨンキー検定も行わず、片側ｔ検定だけを行うという選択肢もあると思います。
　 ただ、この場合には　Ａ＞Ｂ＞Ｃの対立仮説の設定が困難になるかもしれません。
　
　ヨンキー検定を行う場合には、それに引き続いては検定手続きが全く異なる片側ｔ検定は使えません。

　その場合にはウィルコクソンの順位和検定（片側検定）で対比較を行うことは可能だと思います

ヨンキー検定は2群の場合にも適用可能と思います（しらんけど）

２群のヨンキー検定と片側ウィルコクソン検定は同一の検定になります。

どちらでも可ということです。

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

とりあえず、分散分析の例数設計を実施して、
Scheffe等の事後検定の例数設計は実施しないという状況は
現実的にはないということでしょうか？

＞ 例数設計は，結論の妥当性を保証するためのものですから，何を結論とするかがはっきりしていれば自明のことではないですか？

以外，なんの議論が必要か？

ご助言ありがとうございます
参考になりました。

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

Ｒについてお聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

-----------

添付画像にお示しをしました。

表１に示す、０／１のデータフレーム（変数名 df）から、
列ベクトル同士の共起度数（１−１対のカウント）を計算して、
表２に示す、行列（変数名 mx）に格納します。

二重ループで計算すると以下のようになります。

計算結果は対称行列ですので、上三角行列のみ計算しました。

nc <- ncol(df)

mx <- matrix(0, nc, nc) ###### 結果を格納する行列の確保、０クリア

for(i in 1:(nc-1)) {

v1 <- df[,i]

for (j in (i+1):nc) {
v2 <- df[,j]
mx[i,j] <- sum(v1*v2)
}

}

二重ループをapply()関数を使って書き直した場合について、
例示いただければ、大変に勉強になります。

計算結果は対称行列ですので、
二重ループの場合には、上三角行列のみ計算しましたが、
apply()関数を使う場合には、これは条件ではありません。

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

プログラムする前に，「何をやっているのだろうか？」と考えるとよいと思います。
今の場合，データは 0 か 1 で，両方とも 1 の対を数えるわけです。

set.seed(123)
df = as.data.frame(matrix(sample(0:1, 500, replace = TRUE), ncol = 5))

nc <- ncol(df)
mx <- matrix(0, nc, nc) ###### 結果を格納する行列の確保、０クリア

for (i in 1:(nc - 1)) {
  v1 <- df[, i]
  for (j in (i + 1):nc) {
    v2 <- df[, j]
    mx[i, j] <- sum(v1 * v2)
  }
}
mx

m = as.matrix(df)
(t(m) %*% m)

結果は，

> mx
     [,1] [,2] [,3] [,4] [,5]
[1,]    0   22   20   20   21
[2,]    0    0   21   24   24
[3,]    0    0    0   24   20
[4,]    0    0    0    0   22
[5,]    0    0    0    0    0

と

> (t(m) %*% m)
   V1 V2 V3 V4 V5
V1 47 22 20 20 21
V2 22 50 21 24 24
V3 20 21 48 24 20
V4 20 24 24 44 22
V5 21 24 20 22 46

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

大変に勉強になりました。

ループ　→　apply()関数　と、安易に、思考停止していたことを恥ずかしく思います。

(t(m) %*% m)
の計算方法は、興味深く、大変に勉強になりました。

先生がお書きになられたこと
> プログラムする前に，
> 「何をやっているのだろうか？」と考えるとよいと思います。

よいクリスマスプレゼントとして、肝に銘じます。

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

お世話になっており誠に有難うございます。

以下についてご教授頂けないでしょうか。
各種資源の利用回数を二群で比較しました。
一回資源を利用すれば１になります。
一群はわずか、４例でございます。
平均値±SDは資源種類によりまちまちで、低いもので0.5±0.6、高いもので107.5±26.7です。
一群わずか４例なので、正規分布になっているかどうかはとても見れない状態です。

そこでお教え頂きたいのですが、
１）この２群をunpaird t検定で検定するのはおかしいでしょうか？
２）やはり、ノンパラメトリクであり、マンホイットニーU検定を行うべきでしょうか？
３）もしマンホイットニーU検定を行うべきでしたらその理由は、単にノンパラメトリックだからでよいのでしょうか？
決定的な理由がありましたらご教授頂ければ幸甚でございます。

基本的なことがわからず恐縮でございますが、何卒宜しくお願い申し上げます。

　理論的に分布が正規分布に従うはずだ、と言えればｔ検定でよい、と過去に青木先生は言われてます。

　正規分布でない場合でも、U検定よりもWelchのｔ検定の方がよい場合もあると、このサイト内のどこかに書かれていたと思います。どこだったっけ。

鈴木様

ご指導有難うございました。
取り急ぎお礼まで

お世話になります。

2つのグループの比較についておうかがいしたく、投稿いたしました。

2つのチームA、Bがあり、それぞれが10個ずつアイデアを出す課題があるとします。
各アイデアを(A1, A2, .., A10)(B1, B2, .., B10)として、そのアイデアの良さについて複数の評定者に1〜10点で点数をつけてもらいます。
A1とB1など、数字には対応はありません。

最終的にどちらのチームのアイデアがより優れていたかを知りたいと考えています。

この場合、どのような分析を行えばよいでしょうか。
A1とB1などに対応がないので、2要因分散分析は違いそうですし、各アイデアの平均値でt検定をすると、それぞれの評定者の情報が落ちるように思います。

お忙しいところ恐れ入りますが、ご教授いただけますと幸いです。

　平均値の有意差を出したいなら、各アイデアの平均値でｔ検定、で良いのではと思いますが。

　評定者によってどう違う、と言いたいなら、2要因分散分析でよいと思います。（対応はなくてかまわない。）

鈴木康弘様

さっそくのご返答をいただきまして、ありがとうございます。
返信が遅くなり、恐れ入ります。

t検定でよいとのこと、なんとなくですが違和感があります・・・
例えば、各アイデアの評定の平均値が、10人の評定者から得たものなのか、100人の評定者から得たものかによって、平均値の差の統計的な意味が異なってくるように思うのですが、いかがでしょうか。
うまく説明できなくて申し訳ありません。

なお、評定者による違いについては、検討しなくてよいと考えておりますので、2要因分散分析ではなさそうです。

何度もお手数をおかけして恐れ入りますが、ご教授いただけますと幸いです。

> 例えば、各アイデアの評定の平均値が、10人の評定者から得たものなのか、100人の評定者から得たものかによって、平均値の差の統計的な意味が異なってくるように思うのですが、いかがでしょうか。

サンプルサイズを考慮して結論づけるのが「検定」です。

> 評定者による違いについては、検討しなくてよいと考えております

評定者が同じということですか？

なぜ，「検討しなくてよいと考える」のか？

適切な回答を得るには，諸条件を明示されるほうがよいでしょう。

青木先生

返信をいただき、ありがとうございます。

> 適切な回答を得るには，諸条件を明示されるほうがよいでしょう。
失礼いたしました。
評定の状況としては、下記のようになります。
評定者S1が、A1〜B10までの20個についてそれぞれ1〜10点で評定する。
評定者S2も同様に、A1〜B10までの20個についてそれぞれ1〜10点で評定する。
以下すべての評定者で同様。
ですので、「アイデア20個×評定者数」個の評定値が得られることになります。

知りたいのは、A群（A1〜A10）とB群（B1〜B10）の平均値に差があるかどうかです（差がないという帰無仮説を棄却したい）。
なお、A1とB1などに対応はありません（単にA群のアイデア10個、B群のアイデア10個、という感じです）。

私がt検定でイメージしていたのは、A1〜B10のそれぞれについて、各アイデアの評定値の平均値をいったん求めて、そこからA1〜A10の10個、B1〜B10の10個で平均値を比較する、というものでした。
この場合ですと、例えばA1の平均値が10人から得られたものか、100人から得られたものかの情報が落ちてしまいますので、違和感があったのです。

拙い説明で何度もお手数をおかけし、まことに恐れ入ります。
ご教示いただけますと幸いです。

お世話になります。

現在、既に調査が終わったデータがあり、論文を書くことになっております。

例えば、約5000人の回答があり、その中で、500人が介護をしていると答えました。
その人達と、介護をしていない人たちの健康状態の違いを解析したいのですが、コントロールとしては、5000人の中から、500人の介護をしていないと答えた人々をコントロール群としたいと考えております。

特に、性別、年齢をマッチさせたいと考えておりますが、クロス集計などで、このような方法で対照群を設定してよろしいでしょうか？

また、対照群は、介護していないと答えた人の中から、それぞれ性・年齢のマッチする人を、SPSSの無作為抽出の機能で選択しようと考えておりますが、何か他に方法がございますでしょうか？

お忙しい中大変申し訳ないのですが、教えていただければ幸いです。お手数をお掛け致しますが、どうぞよろしくお願い致します。

大変いいにくいことですが，あなたの調査方針・分析方針は不適切です。はっきりいって，涙が出ます。

仮説を立てて，その可否を確かめるために行うのが調査（データ収集）です。
調査をしてしまった後で，解析デザインを作り上げるのは本末転倒です。
対象が500人なら，対照も500人にすべきで，何のために無駄に4000人調査したのでしょうか。あなたは，500＋500の調査データだけを利用しようとしているのですよ。
5000人の回答が得られる調査をするマンパワーがあるなら，介護をしている人，していない人それぞれを選んで（特に，なんとしても最初から性・年齢をマッチさせて）2500人ずつの調査をすればどれだけ素晴らしかったか

> 性別、年齢をマッチさせたいと考えておりますが

調査開始時に，調査対象を選ぶときにその条件を発動すべきでしたね。

そのような調査デザインならば，ケース・コントロール研究として分析できたのに（4000人の無駄もなかったのに）

まあ，ああだこうだ，死んだ子の年を数えても仕方ないので，横断研究ということで割り切って，全データを使って分析しましょ（別に，それが悪いことではなく，データ収集そのものが横断研究だったのだから仕方ないでしょう）
それしかいえません。
あえて，マッチングするためにデータを捨てるということだけはおすすめできません
だって，どのデータを捨てるか決められますか？
いい結果が出なかったら，データを選び直す？そんなこと，してはだめですからね。

お返事、どうもありがとうございました。

たら、れば、で、2500名で性・年齢をマッチしたケース・コントロールができるとさえ考えられるお立場にいらっしゃる方には想像できない程のことが、現実には生じることがあり、その対応に苦慮する者の気持ちなどはご理解いただけないのでしょう。教科書どおりの理想的な調査ができる幸せな環境で研究をされているようで、とてもうらやましく思います。

ご指摘頂きましたように泣いてもどうしようもなく、また、同僚のやったこととはいえ5000名の方々に大変申し訳なく思い、こちらでもアドバイス頂ければと思い投稿させていただいた次第です。

運営上の問題で，No. 22477 の発言者名が別人になっておりましたことを，お詫びいたします。発言者は，青木繁伸 でした。申し訳ありません。
誤った発言者名が表示されてしまった方には，深くお詫び申し上げます。。
また，その記事に対して，発言者が違ったことにより本来の意図は違った発言を誘発してしまったかもしれないことについては，まことに申し訳ありません。
発言を消去いただいても差し支えありません。申し訳ありませんが，対処方よろしくお願いします。
不快な思いをさせてしまったたことについてはお詫び申し上げます。

　マッチングというのは原則としてデータ収集の段階でやるものですね。

　データを集めてしまってから性・年齢の交絡因子を調整したい、という話なら、
多変量解析じゃないでしょうか。

青木先生

お返事、どうもありがとうございました。

「いい結果が出なかったら，データを選び直す？そんなこと，してはだめですからね。」という助言であっても、5000名の人々に大変申し訳ないという思い、ネットですから詳細は書けないなどの制約の中で、それでも、何か、アドバイスが頂きたかったのです。

鈴木康弘様

コメント、どうもありがとうございました。
多変量解析は検討しております。詳細は説明できないのですが、年齢の幅が15歳から80歳とひろく、一般的な調査データではない為、また、介護を例として考えますと、介護している　500 (女性・年齢が50〜70歳台に偏る)と　していない　4500(性別問わず15歳から80歳まで)という間抜けなデータになっており、横断研究として、ロジスティックなどで性・年齢をとにかく調整してもいいものやら、(実際には、その通りの横断研究ではありますが、)そうするしかないのかと思い始めておりますが、扱いに苦慮しております。

再度、検定力分析について調べ直しましたがまだ十分に理解できておりません。
またirr パッケージのオンラインヘルプに引用されている文献ですが取り寄せることができておらず、まだ読んでおりません。
このような状況で質問するのは無礼だと承知していますが、ご教授いただけると幸甚です。

irrパッケージ のN2.cohen.kappaを使用したいのですが

オンラインヘルプの
mrg a vector of marginal probabilities given by raters
ということがわかりません。
そのため
Examples
require(lpSolve)
# Testing H0: kappa = 0.4 vs. HA: kappa > 0.4 (=0.6) given that
# Marginal Probabilities by two raters are (0.2, 0.25, 0.55).

0.2と0.25と0.55の数字が何を現しているかわかりません。
勉強不足で理解が及ばないためご教授頂けると大変助かります。
何卒よろしくお願いいたします。

marginal probabilities というのは，「周辺確率」ということですが，統計学を知らないと難しいかも知れませんね。
集計表で marginal sum というのは「周辺和」といいますが，一般用語では「合計欄」とか「行和」，「列和」といったほうが分かりやすいのでしょう。
irr の kappa2 の example で扱われている diagnose データセットで例を挙げると

> data(diagnoses)
> # Unweighted Kappa for categorical data without a logical order
> kappa2(diagnoses[,2:3])
 Cohen's Kappa for 2 Raters (Weights: unweighted)

 Subjects = 30 
   Raters = 2 
    Kappa = 0.631 

        z = 7.56 
  p-value = 4.04e-14 

この分析例では，rater4 と rater5 を対象としています

> head(diagnoses[,c(4, 5)])
            rater4           rater5
1      4. Neurosis      4. Neurosis
2         5. Other         5. Other
3 3. Schizophrenia 3. Schizophrenia
4         5. Other         5. Other
5      4. Neurosis      4. Neurosis
6 3. Schizophrenia 3. Schizophrenia

集計表は table(diagnoses[,c(4, 5)]) で求められますが，簡単に表示すると

> t=table(diagnoses[,c(4, 5)])
> colnames(t) = rownames(t) = 1:5
> t
      rater5
rater4  1  2  3  4  5
     1  1  0  0  0  1
     2  0  1  0  0  0
     3  0  0  6  0  1
     4  0  0  0 12  1
     5  0  0  0  0  7

これに「周辺和」を付加すると，

> addmargins(t)
      rater5
rater4  1  2  3  4  5 Sum
   1    1  0  0  0  1   2
   2    0  1  0  0  0   1
   3    0  0  6  0  1   7
   4    0  0  0 12  1  13
   5    0  0  0  0  7   7
   Sum  1  1  6 12 10  30

となります。Sum と書かれた行と列にあるのが「周辺和」です。
「周辺確率」はそれをサンプルサイズ（今の例では 30）で割ったものです

> rowSums(t)/30 # rater4 に対する周辺確率
         1          2          3          4          5 
0.06666667 0.03333333 0.23333333 0.43333333 0.23333333 
> colSums(t)/30 # rater5 に対する周辺確率
         1          2          3          4          5 
0.03333333 0.03333333 0.20000000 0.40000000 0.33333333 

それぞれの周辺確率は本来は同じ理論的な確率にしたがうものでしょう。N2.cohen.kappa で指定するのは，この「理論的な（周辺）確率」です。本来は得られたデータから決めるものではありませんが，取りあえずは 2 人の評価者の平均を参考にするということもあるでしょうし，このデータセットでは6人の評価者の平均を使うということもあるでしょう。先行研究の結果を使うこともあるでしょう。

青木先生

　ありがとうございます。
　返信が遅くなり申し訳ございません。
　また周辺確率という新しい言葉を聞いて再度調べ直しています。
　
　先のスレッドを見て、私自身もサンプルサイズなどを検討せず研究を進めていたため反省しています。

　研究方法、研究倫理を見直します。
　統計方法以前の問題ですね。

　
　ありがとうございました。
　またご相談させて頂きます。

Cohen's Kappa のsample sizeを調べています。
Rの'irr'パッケージを取り込むことができたのですが、

> # Testing H0: kappa = 0.7 vs. HA: kappa > 0.7 given that
> # kappa = 0.85 and both raters classify 50% of subjects as positive.
> N.cohen.kappa(0.5, 0.5, 0.7, 0.85)

この0.5と0.7と0.85の数字の意味が分からなく、自験例ではどの数字を入れて良いか判断できません。

初歩的な質問で申し訳ないですがご教授頂けると幸甚です。

オンラインヘルプは見ましたか？

Usage

 N.cohen.kappa(rate1, rate2, k1, k0, alpha=0.05, 
               power=0.8, twosided=FALSE)
Arguments

rate1	
the probability that the first rater will record a positive diagnosis

rate2	
the probability that the second rater will record a positive diagnosis

k1	
the true Cohen's Kappa statistic

k0	
the value of kappa under the null hypothesis

青木先生
ありがとうございます。
オンラインヘルプを確認しました。

再度、質問をさせて頂きたいのですが
rate1
the probability that the first rater will record a positive diagnosis
というのは2択であれば50％→0.5という意味でしょうか。
7択であれば14%→0.14となるのでしょうか。

また
# Testing H0: kappa = 0.7 vs. HA: kappa > 0.7 given that
# kappa = 0.85 and both raters classify 50% of subjects as positive

2行目のkappa=0.85は080以上となるように求めたいので0.80としますが、
1行目のTesting H0:kappa =0.7とありますが、この0.7は何を指しているのでしょうか。

内容を理解しておらず、質問も分かり難いかと思いますがご教授頂けると幸いです。

> rate1
> the probability that the first rater will record a positive diagnosis
> というのは2択であれば50％→0.5という意味でしょうか。
> 7択であれば14%→0.14となるのでしょうか。

N.cohen.kappa のオンラインヘルプの1行目に，

sample size estimator for the Cohen's Kappa statistic for a binary outcome
                                                      ~~~~~~~~~~~~~~~~~~~~

とありますから前者ですね。(N2.cohen.kappa なら more than one category）

「ネガティブな反応」が多い状態とはどういうことか，つまりそのような状態は，評価者がやる気がない（反抗的）ということで，信頼できないデータである可能性があるということ（逆にポジティブな反応が多すぎても，いい加減な肯定と言うことで信頼できないであろう）。50-50 というのは，ポジティブ・ネガティブの基準が平均的（中央値！）ということで客観的で妥当な評価が行われていることを意味する。

> # Testing H0: kappa = 0.7 vs. HA: kappa > 0.7 given that
> # kappa = 0.85 and both raters classify 50% of subjects as positive
> 2行目のkappa=0.85は0.80以上となるように求めたいので0.80としますが、

「0.80以上となるように求めたい」のではなく，個々で指定する値は，実際にあなたのデータで得られたκの値

> 1行目のTesting H0:kappa =0.7とありますが、この0.7は何を指しているのでしょうか。
文字通り帰無仮説のκ。これもあなたが決めるもの。
つまり，あなたは「母数がこれ以上（0.7 以上）かどうかを検定したい」というわけでしょ。「いやいや，0.6 以上あれば十分」ということなら 0.6 を指定するわけです。

「κ統計量の検定・検定力分析」ということではなく，検定，検定力分析の「一般理論」をおさらいすべきですね。基礎がわかれば，応用はおのずとわかります。基礎がわかっていないと，応用はできません。
irr パッケージのオンラインヘルプに引用されている
Cantor, A. B. (1996) Sample-size calculation for Cohen's kappa. Psychological Methods, 1, 150-153.
Flack, V.F., Afifi, A.A., Lachenbruch, P.A., & Schouten, H.J.A. (1988). Sample size determinations for the two rater kappa statistic. Psychometrika, 53, 321-325.
を読むことも必要でしょう。

> N.cohen.kappa(0.5, 0.5, 0.7, 0.85)
[1] 96

これが例題

> N.cohen.kappa(0.5, 0.5, 0.6, 0.85)
[1] 38

帰無仮説の kappa を下げる(0.6)といことは，サンプルサイズは少なくてよいということになる

> N.cohen.kappa(0.5, 0.5, 0.8, 0.85)
[1] 753

kappa を上げる(0.8)ということは，サンプルサイズはより多く必要ということになる

> N.cohen.kappa(0.5, 0.5, 0.7, 0.75)
[1] 1142

実際のデータから得られたκが小さいと（帰無仮説のκとの差が小さくなる），サンプルサイズはより多く必要となる

> N.cohen.kappa(0.2, 0.2, 0.7, 0.85)
[1] 153
> N.cohen.kappa(0.8, 0.8, 0.7, 0.85)
[1] 153

ポジティブ（ネガティブ）率が 0.5 から離れるとサンプルサイズは多く必要

オンラインヘルプには思った以上に必要（有用）な情報が書いてあります。隅から隅まで，読むと良いでしょう。
例えば，デフォルトが twosided = FALSE であるのは，うっかり見逃すことも多いのではないかな？

青木先生

　基礎，内容も十分に理解していないまま使用方法をご質問し申し訳ございませんでした。
　先生の丁寧な解説にて使用方法は理解できましたが、やはり内容を十分に理解できていません。先生のご指摘の通り一般理論から勉強します。ご助言ありがとうございました。
　分からない場合はまたご質問させて頂きます。

再度、検定力分析について調べ直しましたがまだ十分に理解できておりません。
またirr パッケージのオンラインヘルプに引用されている文献ですが取り寄せることができておらず、まだ読んでおりません。
このような状況で質問するのは無礼だと承知していますが、ご教授いただけると幸甚です。

irrパッケージ のN2.cohen.kappaを使用したいのですが

オンラインヘルプの
mrg a vector of marginal probabilities given by raters
ということがわかりません。
そのため
Examples
require(lpSolve)
# Testing H0: kappa = 0.4 vs. HA: kappa > 0.4 (=0.6) given that
# Marginal Probabilities by two raters are (0.2, 0.25, 0.55).

0.2と0.25と0.55の数字が何を現しているかわかりません。
勉強不足で理解が及ばないためご教授頂けると大変助かります。
何卒よろしくお願いいたします。