No.17251　Rを使ったTable集計　　【赤岳】　2012/08/01(Wed) 14:44

本サイトには，大変お世話になっています。
Rの初歩的な使い方について質問させてください。
下表のような表があります。
これを
TREATMENT別に，Day1〜5までの1と0の頻度，1の割合，その割合の95%CIを求めたいのですが，どうしてもうまくできません。すなわち，
CAFFのDay1ですと，1の割合は，76.5%（13/17例），95%CI？？をDay1〜5それぞれについて求めたいと思っています。
できればTableの形できたらよいのですが。
ご教示いただけませんでしょうか。

TREATMENT	PT	Day1	Day2	Day3	Day4	Day5
CAFF	513	1	1	1	1	1
CAFF	515	1	1	1	1	1
CAFF	516	1	0	0	0	0
CAFF	520	1	1	1	1	1
CAFF	523	1	1	1	1	1
CAFF	524	0	1	1	1	1
CAFF	601	1	1	1	1	1
CAFF	603	1	1	0	0	0
CAFF	605	1	1	1	1	1
CAFF	607	1	1	1	1	1
CAFF	608	1	1	1	1	1
CAFF	609	1	1	1	1	1
CAFF	614	0	0	0	0	0
CAFF	616	1	1	1	1	1
CAFF	617	0	0	0	0	0
CAFF	703	1	1	0	0	0
CAFF	706	0	1	1	1	1
PLA	101	1	1	1	1	1
PLA	114	1	1	1	1	1
PLA	115	0	0	0	0	0
PLA	116	1	1	1	0	0
PLA	202	0	0	0	0	0
PLA	203	0	1	1	1	1
PLA	206	1	1	1	1	1
PLA	209	1	1	1	1	1
PLA	210	0	0	0	0	0
PLA	212	1	1	0	0	1
PLA	214	0	0	0	0	0
PLA	301	0	0	0	0	0
PLA	304	1	1	1	1	1
PLA	308	0	1	0	0	0
PLA	309	0	1	1	1	0
PLA	311	0	0	0	0	0
PLA	316	1	1	1	1	1

No.17253　Re: Rを使ったTable集計　　【青木繁伸】　2012/08/01(Wed) 21:27

自分でちょっとしたプログラムを書かねばなりません。
とはいっても，既存の関数をつなぎ合わせるだけですが（既存の関数を使わないで書くこともできますが，面倒です）。
より汎用性を持ったプログラムにすることもできますが，汎用プログラムに値しないでしょうからこのまま。
データフレーム名も d に固定して書きました。
func <- function(y) {
	func2 <- function(x, n) {
		ans <- prop.test(x, n)
		return(c(p=unname(ans$estimate), lcl=ans$conf.int[1], ucl=ans$conf.int[2]))
	}
	a <- addmargins(table(d$TREATMENT, y))
	return(c(func2(a[1,2], a[1,3]), func2(a[2,2], a[2,3])))
}
round(sapply(d[,3:7], func), 3)
実行結果は，
> round(sapply(d[,3:7], func), 3)
     Day1  Day2  Day3  Day4  Day5
p   0.765 0.824 0.706 0.706 0.706
lcl 0.498 0.558 0.440 0.440 0.440
ucl 0.922 0.953 0.886 0.886 0.886
p   0.471 0.647 0.529 0.471 0.471
lcl 0.239 0.386 0.285 0.239 0.239
ucl 0.715 0.847 0.761 0.715 0.715
上3行が CAFF, 下3行が PLA。
あとは，必要ならば細かな調整をすればよいでしょう（プログラムのどこで何をしているかを把握できれば）。
No.17255　Re: Rを使ったTable集計　　【赤岳】　2012/08/01(Wed) 23:36

青木先生，いつもご指導ありがとうございます。
すごい複雑なのに驚いています。
xtabsなどの既成の関数で，もっと簡単に書けるものと思っていました。
先生の式を元に，拡張していきたいと思います。
本当にありがとうございました。
No.17257　Re: Rを使ったTable集計　　【青木繁伸】　2012/08/02(Thu) 00:30

ちっとも複雑ではないでしょう。

それに，xtabs だけではあなたの望む結果は出せません。
計算の一部として使える xtabs の機能は，この場合は table と基本的には同じ（xtabs の引数を見れば分かる）。
> xtabs(~TREATMENT+Day1, data=d)
         Day1
TREATMENT  0  1
     CAFF  4 13
     PLA   9  8
> table(d$TREATMENT, d$Day1)
      
        0  1
  CAFF  4 13
  PLA   9  8
なお，先ほどのプログラムは，もう少し短くなります。
table 関数と prop.test 関数しか使いません。
> func <- function(y) {
+ 	func2 <- function(x) {
+ 		ans <- prop.test(x[2], sum(x))
+ 		return(c(p=unname(ans$estimate), lcl=ans$conf.int[1], ucl=ans$conf.int[2]))
+ 	}
+ 	a <- table(d$TREATMENT, y)
+ 	return(c(func2(a[1,]), func2(a[2,])))
+ }
> round(sapply(d[,3:7], func), 3)
     Day1  Day2  Day3  Day4  Day5
p   0.765 0.824 0.706 0.706 0.706
lcl 0.498 0.558 0.440 0.440 0.440
ucl 0.922 0.953 0.886 0.886 0.886
p   0.471 0.647 0.529 0.471 0.471
lcl 0.239 0.386 0.285 0.239 0.239
ucl 0.715 0.847 0.761 0.715 0.715
また，このプログラムで，prop.test のかわりに binom.test を使えば，二項分布による正確な信頼限界が求められるというのも，このようなアプローチをとるメリットでしょう。
> round(sapply(d[,3:7], func), 3)
     Day1  Day2  Day3  Day4  Day5
p   0.765 0.824 0.706 0.706 0.706
lcl 0.501 0.566 0.440 0.440 0.440
ucl 0.932 0.962 0.897 0.897 0.897
p   0.471 0.647 0.529 0.471 0.471
lcl 0.230 0.383 0.278 0.230 0.230
ucl 0.722 0.858 0.770 0.722 0.722
No.17258　Re: Rを使ったTable集計　　【赤岳】　2012/08/02(Thu) 13:14

青木先生，
ありがとうございます。
理解するため，先生の下記プログラム実行した後，prop.testを実行しました。
a <- table(d$TREATMENT,d$TREATMENT, d$Day1)
0 1
CAFF 4 13
PLA 9 8

> prop.test(a)

2-sample test for equality of proportions with continuity
correction

data: a
X-squared = 1.9927, df = 1, p-value = 0.1581
alternative hypothesis: two.sided
95 percent confidence interval:
-0.66431747 0.07608218
sample estimates:
prop 1 prop 2
0.2352941 0.5294118

こういう結果が得られましたが，ここの95%CIは一つしか出てきていません。
これは一体何の値でしょうか。
また，
prop1（CAFF）, prop2（PLA）のそれぞれの95%CIはどこに出力されるのでしょうか。
何か大きく勘違いしているかも知れません。
お世話おかけし，申し訳ございません。

追記：二項分布のbinom.testの件，非常に勉強になりました。ありがとうございます。
No.17259　Re: Rを使ったTable集計　　【青木繁伸】　2012/08/02(Thu) 18:12

prop.test も指定法により何通りもの検定をするので注意。

a <- table(dTREATMENT,d$TREATMENT, d$Day1)
prop.test(a)

は，「二群の比率の差の検定」を行うことになります。
If p is NULL and there is more than one group, the null tested is that the proportions in each group are the same.
結果の出力の冒頭にも書いてあるので，注意しましょう。
"2-sample test for equality of proportions with continuity correction" とね。

しかしそれは，明らかにあやまり。
正しくは，「母比率の検定」をするのです。
If there is only one group, then the null tested is that the underlying probability of success is p, or .5 if p is not given.

あなたが挙げた「1の割合は，76.5%（13/17例）」というのは，CAFF について比率と母比率の信頼区間を求めるということです。
No. 17253 のプログラムの方が分かりやすいですが，prop.test(x, n) というのがそれです。x，n はひとつの数値です。

a <- table(dTREATMENT,dDay1) に続いてやるべきは，まず，CAFF についての分析（検定）x は a[1,2], n は a[1,1]+a[1,2] です。
prop.test(a[1,2], sum(a[1,])) は，a[1,2]=13，sum(c(a[1,1], a[1,2]))=17 なので prop.test(13, 17) とおなじです。

実行例は
> prop.test(a[1,2], sum(a[1,]))

	1-sample proportions test with continuity correction ★★ タイトルが違うことに注意

data:  a[1, 2] out of sum(a[1, ]), null probability 0.5 
X-squared = 3.7647, df = 1, p-value = 0.05235
alternative hypothesis: true p is not equal to 0.5 
95 percent confidence interval:
 0.4976163 0.9217688 
sample estimates:
        p 
0.7647059
これは，CAFF の場合です。PLA については，prop.test(a[2,2], sum(a[2,])) です。
> prop.test(a[2,2], sum(a[2,]))

	1-sample proportions test with continuity correction

data:  a[2, 2] out of sum(a[2, ]), null probability 0.5 
X-squared = 0, df = 1, p-value = 1
alternative hypothesis: true p is not equal to 0.5 
95 percent confidence interval:
 0.2385724 0.7146614 
sample estimates:
        p 
0.4705882
それぞれの分析（検定）において，prop.test が返すオブジェクト中に 比率の推定値は estimate, 信頼限界は conf.int として返されるので ansestimateやansconf.int のようにして取り出すのです。

そして，これを Day1 〜 Day5 までやるので，sapply 関数を使うのです。

これだけ面倒くさいことを，関数をつなぎ合わせるだけで（あれだけのプログラム量で）できるんですから，簡単ですよね。
No.17260　Re: Rを使ったTable集計　　【赤岳】　2012/08/02(Thu) 19:16

青木先生，

詳細に，噛み砕いてのご説明，大変恐縮いたします。
非常によく分かりました。
prop.test(a）は，「二群の比率の差の検定」
prop.test(x, n)は，一つの群のCAFF について比率と母比率の信頼区間
Day1 〜 Day5 までやるのにsapply 関数を使う
以上3点のことが分からず，悩んでおりました。
functionについては，まだぜんぜん馴染めず，悩まされています。
function文があると，ついつい引いてしまいます。慣れるしかないと思っています。

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

● 「統計学関連なんでもあり」の過去ログ--- 045 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る