デ−タ解析で出会う統計的問題
  (久保担当ぶんの補足 2 --- R プログラムたち)

とりあえずの再現手順

1. ディレクトリ (フォルダ) を作り， 必要なファイル COMMON.R , samples.R , grouping.R , estimate.R たちを保存．
2. そのディレクトリ (フォルダ) で R を起動する． あるいはいったん起動してから setwd(ディレクトリ名) で移動する．
3. 例題生成: source("samples.R")
4. その例題の図を作ってみる:

   > plot.frame()
   > plot.samples.box()
   > plot.samples()
   

   これらの作図関数は COMMON.R で定義されている．
5. ポアソンモデルの推定: source(estimate.R) して以下のように推定計算する．

   > source("estimate.R")
   > list <- estimate.poisson(samples)
   # evaluating Model(A+B+C)     AIC = 214.1 ...
   # evaluating Model(A+B)(C)    AIC = 211.1 ...
   # evaluating Model(A+C)(B)    AIC = 216.1 ...
   # evaluating Model(A)(B+C)    AIC = 210.9 ...
   # evaluating Model(A)(B)(C)   AIC = 212.0 ...
   

ファイルの解説:

---

ファイル: COMMON.R

内容:

(共用部分)

使いかた:

これ単体では特に使いみちなし． 例題生成のためのパラメーターセットなどが 定義されている．

	parameter <- data.frame(
		"treatment" = as.character(c("A", "B", "C")),
		"mean" = c(2.5, 3.5, 3.5),
		"n" = as.integer(c(20, 10, 20))
	)

---

ファイル: samples.R

内容:

Poisson 乱数の例題生成

使いかた:

R のコマンドラインで

	source("samples.R")

とすると samples.RData というファイルに乱数例が保存される (このファイル名は COMMON.R で定義されている)． このデータは load("samples.RData") で読みこめる．

---

ファイル: grouping.R

可能なグループわけを生成する generate.groups() 関数など; (注意) gregmisc package が必要なので CRAN からダウンロードしてインストールしてください．

使いかた:

R のコマンドラインで

	source("grouping.R")

とすると定義された関数たちが読みこまれる．

generate.groups(水準のvector) 関数

generate.groups(c("A", "B")) とすると

	[[1]]
	[[1]]$tag
	[1] "(A+B)"

	[[1]]$map
	      A       B 
	"(A+B)" "(A+B)" 

	[[1]]$number.groups
	[1] 1

	[[2]]
	[[2]]$tag
	[1] "(A)(B)"

	[[2]]$map
	    A     B 
	"(A)" "(B)" 

	[[2]]$number.groups
	[1] 2

という出力が得られる (データ構造は list)． これは A と B の二水準のときに (A+B) (A と B は同じ) あるいは (A)(B) (A と B は異る) といったグループわけができることを示している． 出力の list 内の各要素のうち tag はグループわけの名前， map (写像) は各水準の所属するグループをあらわす．

使用例: 4 水準のグループわけを列挙

> sapply(generate.groups(1:4), function(e) e$tag)
 [1] "(1+2+3+4)"    "(1+2+3)(4)"   "(1+2+4)(3)"   "(1+3+4)(2)"   "(1)(2+3+4)"  
 [6] "(1+2)(3+4)"   "(1+3)(2+4)"   "(1+4)(2+3)"   "(1+2)(3)(4)"  "(1+3)(2)(4)" 
[11] "(1+4)(2)(3)"  "(1)(2+3)(4)"  "(1)(2+4)(3)"  "(1)(2)(3+4)"  "(1)(2)(3)(4)"

---

ファイル: estimate.R

内容:

すべての可能なグループわけに関して， 一般化線形モデル (ポアソン回帰) を glm() 関数で評価し， その結果の list を返す

使いかた:

R のコマンドラインで

	source("grouping.R")

とすると定義された関数たちが読みこまれる． 上述の samples.R で生成された samples を読みこみ， 計算する． 使用例:

> list.results <- estimate.poisson(samples)
# evaluating Model(A+B+C)     AIC = 195.5 ...
# evaluating Model(A+B)(C)    AIC = 194.7 ...
# evaluating Model(A+C)(B)    AIC = 197.3 ...
# evaluating Model(A)(B+C)    AIC = 192.9 ...
# evaluating Model(A)(B)(C)   AIC = 194.8 ...

list.results に glm() の推定結果が入っている．

> labels(list.results) # list.results にいくつ結果が入っているか見る
[1] "1" "2" "3" "4" "5"
> list.results[[4]] # 4 番目の推定結果を表示
$tag
[1] "(A)(B+C)"

$map
      B       C       A 
      "(B+C)" "(B+C)"   "(A)" 

      $number.groups
      [1] 2

      $glm

      Call:  glm(formula = samples$observed ~ tr.mapped - 1, family = poisson(link = log)) 

      Coefficients:
        tr.mapped(A)  tr.mapped(B+C)  
	         0.875           1.243  

		 Degrees of Freedom: 50 Total (i.e. Null);  48 Residual
		 Null Deviance:      189 
		 Residual Deviance: 50.2         AIC: 193 

---

index に戻る