前回のあらすじ

どういうわけか，ANOVAコーディングによる対数線形モデルの推定値をまともに表示してくれる統計ソフトはほとんどありません。残念ながらRのloglm()やglm()も同様です。そこで，glm()の結果オブジェクトを利用して，必要な数値を自動で取得できるようにしてみます。

下準備

まず，データは前回と同じものを使用しましょう。

sample_data <- '
性別 年齢   分類      Freq
男   ~19   1.親族       92
男   ~19   2.強盗殺人  110
男   ~19   3.その他    144
男   20~29 1.親族      175
男   20~29 2.強盗殺人  191
男   20~29 3.その他    389
男   30~39 1.親族      212
男   30~39 2.強盗殺人   85
男   30~39 3.その他    217
男   40~   1.親族      263
男   40~   2.強盗殺人   34
男   40~   3.その他    188
女   ~19   1.親族       18
女   ~19   2.強盗殺人    2
女   ~19   3.その他      8
女   20~29 1.親族       68
女   20~29 2.強盗殺人    7
女   20~29 3.その他     22
女   30~39 1.親族       69
女   30~39 2.強盗殺人    4
女   30~39 3.その他     17
女   40~   1.親族       46
女   40~   2.強盗殺人    3
女   40~   3.その他      2'

df <- read.table(textConnection(sample_data),head=T)

このデータを対数線形モデルで分析しますが，モデルをできるだけ単純にするために，有意でない効果は省いておくことにします。

# ANOVAコーディング
> res <- glm( Freq~性別*年齢*分類, family='poisson', data=df,
                contrasts=list(
                    性別=contr.sum, 
                    年齢=contr.sum, 
                    分類=contr.sum))

> require(car)
> Anova(res)
Analysis of Deviance Table (Type II tests)

Response: Freq
               LR Chisq Df Pr(>Chisq)    
性別            1616.40  1  < 2.2e-16 ***
年齢             198.59  3  < 2.2e-16 ***
分類             263.06  2  < 2.2e-16 ***
性別:年齢         21.85  3  7.025e-05 ***
性別:分類        166.26  2  < 2.2e-16 ***
年齢:分類        195.51  6  < 2.2e-16 ***
性別:年齢:分類    10.40  6     0.1086    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

この結果から，性別×年齢×分類の交互作用は有意ではありませんので，これは省いてしまいます。

# 主効果と１次の交互作用までを含むようにモデルを修正
# res <- update(res, .~. -性別:年齢:分類) としても結果は同じ

> res <- glm( Freq~ (性別+年齢+分類)^2, family='poisson', data=df,
                contrasts=list(
                    性別=contr.sum, 
                    年齢=contr.sum, 
                    分類=contr.sum))

> Anova(res)
Analysis of Deviance Table (Type II tests)

Response: Freq
          LR Chisq Df Pr(>Chisq)    
性別       1616.40  1  < 2.2e-16 ***
年齢        198.59  3  < 2.2e-16 ***
分類        263.06  2  < 2.2e-16 ***
性別:年齢    21.85  3  7.025e-05 ***
性別:分類   166.26  2  < 2.2e-16 ***
年齢:分類   195.51  6  < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

考え方

基本的には，分散共分散行列と係数行列のそれぞれについて，行列全体を要因ごとに分割し，不足部分を補うという処理の繰り返しです。１次の交互作用までを含めたサンプルのモデルについて，vcov()で得られる分散共分散行列には以下の情報が含まれています。

しかし，ANOVAコーディングの場合には，これだけの分散共分散行列が必要です。

ここで，グレーになっている部分が省略されている部分です。こうしてみると，結構な部分が省略されています。

さらに行方向のみ，色をつけてみました。オレンジ色の部分は，それぞれの要因内で，その上の行の合計を0から引くことによって求められます。たとえばA2の行は0-A1，B4の行は0-(B1+B2+B3)です。交互作用についても，たとえばA2:B1は0-A1:B1，B4:C1は0-(B1:C1+B2:C1+B3:C1)です。

また，青色の部分は，その上の行をブロック単位で合計して0から引くことによって求められます。たとえば，[A1:B4 ~ A2:B4]の部分は0-( [A1:B1 ~ A2:B1]+[A1:B2 ~ A2:B2]+[A1:B3 ~ A2:B3] )，[A1:C3 ~ A2:C3]の部分は0-( [A1:C1 ~ A2:C1]+[A1:C2 ~ A2:C2] )，最後の(B1:C3 ~ B4:C3)の部分は0-( [B1:C1 ~ B4:C1]+ [B1:C2 ~ B4:C2] )となります。

また，列方向も同様の手順で求めることができ，行方向で値を埋めてから列方向で埋めると，たとえばこの表の一番右下の部分も値が入ることになります。

この表の太線で示してある部分は，それぞれ主効果や交互作用の境目の部分で，この表は上から2, 4, 3, 8, 6, 12行に分割できます。また，交互作用の部分（10行目以降）はそれぞれ，8 =2*4, 6 =2*3, 12 =4*3に分割できます。

これをvcov()の分散共分散行列でいうと，主効果・交互作用の境目は上から1, 3, 2, 3, 2, 6行，交互作用の部分は3 =1*3, 2 =1*2, 6 = 3*2で，それぞれ各要因の水準数が1ずつ少ないだけですので，この部分を計算によって埋めてやることになります。

一番複雑なB*C（年齢:分類）の交互作用のところだけを取り出してみると，より具体的には次のようにC1とC2の境目の部分で分割し，分割によってできた２つのブロックを合算して0から引けばC3のブロックを作成できます。

そして，出来上がったそれぞれのブロックで，すでに値の埋まっている３行を合成してB4に関する行を作成します。

もっと要因が多くなれば，分割して省略部分を補充したものを，さらに要因の水準に応じて分割し，またそこで省略部分を補充し，そしてそれをさらに分割し・・・と進めていき，必要な部分がすべて補充されたらそれらを合体させることになります。

このような，分割して補充したものをさらに分割し・・・という処理は，再帰処理を使えばうまくいきそうできそうです。また，分割によって行ブロックが増えた場合，それぞれのブロックに対して繰り返し同じ処理を行うことになりますが，この手の繰り返し作業を行う関数として，Rにはlapply()やsapply()などのapply()系関数が備わっています。

もっと効率的な方法もあるかもしれませんが，今回はこうした考え方に基づいてスクリプトを作成してみました。

係数算出スクリプト

当初の想定ではもっとシンプルになるはずだったのですが，出力の見た目を整えたりしていたら，思ったより長くなってしまいました。作成したスクリプトは以下の通りです。自分の勉強も兼ねてできるだけforを使わずapply系を多用したため，ちょっとわかりづらくなってしまったかもしれません。また，処理に無駄が多いのでしょう，飽和モデルのように２次の交互作用があったりすると，（数秒程度ですが）やや計算に時間がかかります。

show.all.coef <- function(lmobj){ 
	####################################
	# 補助的な処理関数
	
	#### 行列をn分割
	divide.row<-function(mat, n_split){
		if(!is.data.frame(mat)) mat<- as.data.frame(mat)
		n_row <- nrow(mat)
		split(mat, rep(1:n_split, each= n_row/n_split))	
	}
	
	#### 行列を要因ごとに分割
	split.row<-function(vcovmat, lvvec){
		x.vcov<-vcovmat
		res.list<-list()
	
		sapply(lvvec, function(X) {
			l <- Reduce('*', X)
			res.list <<- c(res.list, list(c(list(X), list(as.data.frame(x.vcov[1:l,,drop=F])))))
			x.vcov <<- x.vcov[-(1:l),]
		})
		return(res.list)
	}
	
	#### 各要因の水準を処理
	get.lv.vec<- function(lmobj){
		x.terms <- attr(terms(lmobj$terms), "term.labels")
		x.levels <- lmobj$xlevels
		lv.vec <- list()	
	
		sapply(x.terms, function(X) {
			tms <- unlist(strsplit(X,':'))
			lvs <- numeric(0)
			sapply(tms, function(Y){
				lvs <<- c(lvs, length(x.levels[[Y]])-1)
			})
			lv.vec <<- c(lv.vec, list(lvs))
		})
		names(lv.vec) <- x.terms
		return(lv.vec)
	}
	
	#### 完全な行名を作成
	get.rownames<- function(lmobj){
		x.terms <- attr(terms(lmobj$terms), "term.labels")
		x.levels <- lmobj$xlevels

		sapply(1: length(x.levels), function(X) {
			x.levels[[X]] <<- paste(names(x.levels[X]), x.levels[[X]])
		})
		
			
		# 各termを要因に分解
		t.list <- sapply(x.terms, function(X) {
			# 多次元の処理
			x <- strsplit(X,':')
			x <- lapply(x, function(Y) x.levels[Y])
		})
		# termごとに組み合わせを作成
		x.comb <- sapply(t.list, function(X) {
			n <- length(X)
			expand.grid(X)
		})
		# 組み合わせリストを結合
		res <- sapply(x.comb, function(X){
			if(ncol(X)==1){
				as.character(unlist(X))
			} else {
				sapply(1:nrow(X), function(Y){
					Reduce(function(a,b){paste(a,b,sep=':')},X[Y,])
				})
			}
		})
		return (unlist(res))
	}
	
	#### 結果清書処理
	res.show.1<-function(mat, n, l){
		cat('\n',n[[1]],'\n','========================\n',sep='')
		n.row <- Reduce('*', (l[[ n[[1]] ]]+1))
		print(format(mat[1:n.row, ], justify='left'))
		cat('---\n',"Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",sep='')
		if(nrow(mat[-(1:n.row),])>1){ 
			res.show.1(mat[-(1:n.row),], n[-1], l[-1])
		}
	}
	
	res.show<-function(mat, lmobj){
		x.terms <- (get.lv.vec(lmobj)) 
		x.names<-names(x.terms)
		
		print(format(mat[1, ], justify='left'))
		cat('---\n',"Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",sep='')
		mat<-mat[-1, ]
		res.show.1(mat, x.names, x.terms)	
	}
	
	
	####################################
	# 主要処理の関数
	
	#### 省略部分の算出（多次元用）
	calc.missing.1 <- function(mat, n){
		if(length(n[-1]) == 0) mat
		else {
			n <- n[-1]
			x <- divide.row(mat, n[1])
			missing <- Reduce('+', x) * (-1)
			x <- c(x, list(missing))
			
			x <- lapply(x,function(X){
				calc.missing.1(X, n)
			})
			x <- Reduce(rbind, x)
			return (x)
		}
	}
	
	#### 省略部分の算出
	calc.missing <- function(mat, level){
		level=rev(level)
		
		# 行を分割
		x <- divide.row(mat, level[1])
	
		# 省略されている行を計算して追加
		missing <- Reduce('+', x) * (-1)
		x <- c(x, list(missing))
	
		# ２要因以上の場合は分割されたブロックごとに再分割して省略部分を計算
		x <- lapply(x, function(X) {
			calc.missing.1(X, level)
		})
		
		# 結果を１つにまとめる
		x <- Reduce(rbind, x)
	}
	
	####  標準誤差の算出
	calc.stderr<-function(lmobj){
		
		# 切片を除去
		x <- vcov(lmobj)[-1,-1]
	
		# 分散共分散行列の行を効果ごとに分割
		x <- split.row(x, get.lv.vec(lmobj))
	
		# 分散共分散行列の省略部分を算出して加える
		x <- lapply(x, function(X) calc.missing(X[2],X[[1]]))
	
		# 分散共分散行列を１つにまとめる
		x <- Reduce(rbind, x)
	
		# 行と列の入れ替え
		x<- t(x)
	
		# 入れ替え後の共分散行列で再度省略部分を算出して結合
		x <- split.row(x, get.lv.vec(lmobj))
		x <- lapply(x, function(X) calc.missing(X[2],X[[1]]))
		x <- Reduce(rbind, x)
	
		# 分散だけを抽出
		x <- diag(as.matrix(x))
	
		# 切片の情報を戻す
		x <-rbind(vcov(lmobj)[1,1],as.data.frame(x))
	
		# 分散の正の平方根が標準誤差
		return(sqrt(x))
	}
	
	####  係数の算出
	calc.coef <- function(lmobj){
		# 係数を抽出
		x <-as.data.frame(lmobj$coeff[-1])
		x$dummy <- 0 # １列だけだとエラーになるのでダミーで列を追加
	
		# 分散共分散行列と同様に省略部分を算出
		x <- split.row(x, get.lv.vec(lmobj))
		x <- lapply(x, function(X) calc.missing(X[2],X[[1]]))
		x <- Reduce(rbind,x)
	
		# 切片の情報を戻す
		x <- rbind(lmobj$coeff[1],x[1])
	
		return(x)
	}
	
	####################################
	# メイン処理
	
	#### 標準誤差の算出
	full.stderr <- calc.stderr(lmobj)
	
	#### 係数の算出
	full.coef <- calc.coef(lmobj)
	
	#### z値の算出
	full.z.value <- full.coef/full.stderr
	
	#### 両側検定でのzの確率の算出
	full.Pr <- sapply(full.z.value, function(X) pnorm(abs(X),0,1,lower.tail=F)*2)
	
	#### 有意記号
	full.sig <- sapply(full.Pr, function(X) {
		if(X < .001) '***'
		else if(X < .01) '**'
		else if(X < .05) '*'
		else if(X < .10) '.'
		else ' '
	})
	
	#### 全部をまとめる
	full.table <- cbind(full.coef,full.stderr,full.z.value,full.Pr,full.sig)
	
	#### 行名と列名
	rownames(full.table)<-c('(Intercept)',get.rownames(lmobj))
	colnames(full.table)<-c('Estimate', 'Std. Error', 'z value', 'Pr(>|z|)','')
	
	#### 多少見やすく表示
	res.show(full.table, lmobj)
}

実行結果

このスクリプトをサンプルデータの結果オブジェクトに対して実行してみます。この結果と，coef(loglm(Freq~(性別+年齢+分類)^2,data=df))を実行した場合の結果で係数を比較してみてください。glmとloglmでは係数の計算方法が違うのか，まったく同じ数字にはなりませんが，小数点５桁目程度までは一致した値が得られたはずです。

> show.all.coef(res)

            Estimate Std. Error  z value Pr(>|z|)    
(Intercept) 3.702562 0.05527773 66.98108        0 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

性別
========================
         Estimate Std. Error   z value      Pr(>|z|)    
性別 女 -1.308874 0.05403883 -24.22098 1.337453e-129 ***
性別 男  1.308874 0.05403883  24.22098 1.337453e-129 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

年齢
========================
             Estimate Std. Error   z value     Pr(>|z|)    
年齢 ~19   -0.3773068 0.08139305 -4.635615 3.558780e-06 ***
年齢 20~29  0.6119044 0.05662153 10.806921 3.192040e-27 ***
年齢 30~39  0.1962872 0.06074367  3.231402 1.231848e-03 ** 
年齢 40~   -0.4308847 0.07615437 -5.658043 1.531089e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

分類
========================
                  Estimate Std. Error   z value     Pr(>|z|)    
分類 1.親族      0.7756114 0.05882659 13.184709 1.074682e-39 ***
分類 2.強盗殺人 -0.9202835 0.09299232 -9.896339 4.317930e-23 ***
分類 3.その他    0.1446721 0.06815129  2.122808 3.376994e-02 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

性別:年齢
========================
                     Estimate Std. Error   z value    Pr(>|z|)   
性別 女:年齢 ~19   -0.0938277 0.08132109 -1.153793 0.248585040   
性別 男:年齢 ~19    0.0938277 0.08132109  1.153793 0.248585040   
性別 女:年齢 20~29  0.1604780 0.05569295  2.881478 0.003958144 **
性別 男:年齢 20~29 -0.1604780 0.05569295 -2.881478 0.003958144 **
性別 女:年齢 30~39  0.1476254 0.05667241  2.604891 0.009190349 **
性別 男:年齢 30~39 -0.1476254 0.05667241 -2.604891 0.009190349 **
性別 女:年齢 40~   -0.2142758 0.06604856 -3.244215 0.001177746 **
性別 男:年齢 40~    0.2142758 0.06604856  3.244215 0.001177746 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

性別:分類
========================
                          Estimate Std. Error    z value     Pr(>|z|)    
性別 女:分類 1.親族      0.6349615 0.05728629  11.084005 1.500097e-28 ***
性別 男:分類 1.親族     -0.6349615 0.05728629 -11.084005 1.500097e-28 ***
性別 女:分類 2.強盗殺人 -0.4017288 0.09006429  -4.460468 8.178097e-06 ***
性別 男:分類 2.強盗殺人  0.4017288 0.09006429   4.460468 8.178097e-06 ***
性別 女:分類 3.その他   -0.2332327 0.06629483  -3.518113 4.346274e-04 ***
性別 男:分類 3.その他    0.2332327 0.06629483   3.518113 4.346274e-04 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

年齢:分類
========================
                              Estimate Std. Error    z value     Pr(>|z|)    
年齢 ~19:分類 1.親族       -0.36316013 0.06473808 -5.6096835 2.026970e-08 ***
年齢 20~29:分類 1.親族     -0.41655689 0.05094462 -8.1766614 2.918166e-16 ***
年齢 30~39:分類 1.親族      0.13821529 0.05389808  2.5643826 1.033595e-02 *  
年齢 40~:分類 1.親族        0.64150174 0.06102525 10.5120710 7.600585e-26 ***
年齢 ~19:分類 2.強盗殺人    0.48237457 0.06822017  7.0708496 1.539880e-12 ***
年齢 20~29:分類 2.強盗殺人  0.29990646 0.05722116  5.2411812 1.595520e-07 ***
年齢 30~39:分類 2.強盗殺人 -0.09586531 0.06867398 -1.3959481 1.627301e-01    
年齢 40~:分類 2.強盗殺人   -0.68641571 0.09083564 -7.5566786 4.134921e-14 ***
年齢 ~19:分類 3.その他     -0.11921443 0.05944967 -2.0053003 4.493095e-02 *  
年齢 20~29:分類 3.その他    0.11665043 0.04596720  2.5376883 1.115873e-02 *  
年齢 30~39:分類 3.その他   -0.04234997 0.05321073 -0.7958916 4.260951e-01    
年齢 40~:分類 3.その他      0.04491397 0.06248561  0.7187891 4.722709e-01    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1