アールメカブ

• 追加された行はこの色です。
• 削除された行はこの色です。
• R_fromOldHtml1 へ行く。
• R_fromOldHtml1 の差分を削除

個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-mこの部分を"@"に変更下さいias.tokushima-u.ac.jp
 までご連絡ください．
個人メモを整理中です．不適切な記述が多々あるかと思います．お気づきの際は，ishida-m(この部分を"@"に変更下さい)ias.tokushima-u.ac.jp  までご連絡ください．

#contents

*   [[ ]] の使い方 [#z6109ecf]


例えば以下のようなデータで ([[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] p.200)

  > z
    0.8 1 1.2 1.4 1.6 2.5
    0 6   4 2   2   0   0  
    1 1   1 4   4   4   2  
  > z[[1]] // 一列目の一行目
  [1] 6
  > z[[2]] // 一列目の二行目
  [1] 1
  > z[[3]] // 二列目の一行目
  [1] 4
  > z[[4]] // 二列目の二行目
  [1] 1
  > z[[5]] // 三列目の一行目
  [1] 2
  > z[[6]] // 三列目の二行目
  [1] 4
  > 
  > z

このとき以下は,行名を

  > dimnames(z)[[1]]
  [1] "0" "1"

次は列名を表示する.

  > dimnames(z)[[2]]
  [1] "0.8" "1"   "1.2" "1.4" "1.6" "2.5"

今は二次元のデータだから,三次元を指定するとエラーとなる.
  > dimnames(z)[[3]]
  Error in dimnames(z)[[3]] : subscript out of bounds





*  boxplot [#qb7738a5]

例えばカッコーの巣であれば，

  > source("/source/statistics/R/book/DAAG/data/cuckoos.R")
  > cuckoos
      length breadth       species  id
  1     21.7    16.1  meadow.pipit  21
  ...
  120   21.0    16.0          wren 238
  > boxplot(length ~ species, data=cuckoos)
  > 

ちなみに，外れ値はレンジで指定するが，デフォルトでは range = 1.5 となっている。詳細は『データサイエンス入門』p.14 

またrug() で，y軸右(side=4)ないし左(side=2)に，プロットに対応したバーコードを表示できる。

ラグプロットについては[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?%A5%E9%A5%B0%A5%D7%A5%ED%A5%C3%A5%C8]] を参照

なお boxplot と split, rug を併用した例は [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]   p.225






*  リストの先頭要素を取る [#r796e7ac]


   x[1:2]
 
  $string
   [1] "a" 
 
  $integer 
  [1] 1 2 
 > x["string"] 
 
  $string
   [1] "a"
 
   > x[-1] 
 
  $integer 
  [1] 1 2 
 
  $float 
  [1] 3.14 









*  qt下側確率 [#uf69c824]

構文は qt(下側確率，自由度) 坪田 p.29

例えば qt(0.05,10) とすると，qt()は分位点関数なので，正規曲線の総面積の累積が5%となる x の座標点 
  > qt(0.05,10)
  [1] -1.812461
を与える。上側は 
  > qt(0.95,10)
  [1] 1.812461








*  Excel データの読み込み [#va762f1e]

RODBCパッケージを利用する。

以下は，excelファイルを読み込むだけでなく，更新なども行なう。
  library(RODBC)
  connExcel 
このあと

  data 

などとクエリーを行なえる。使い終われば
  odbcClose(connExcel)









*  新しいグラフィックデバイス [#t51cebd3]


[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] , p.22,35
X11() を実行。サイズを指定するには X11(width=7.5,height=4)

現在アクティヴなグラフィックデバイスを閉じるには

  dev.off()

特定のデバイス, 例えば ３ 番を閉じるには

  dev.off(3)

デバイスの変更, 例えば ２ 番に切替えるなら

  dev.set(2)











*  特定の行の取り出し [#l372d554]

例えば [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]] , p.179 の砂糖の平均を，コントロール群だけ取り出すなら

  > sugar
     weight     trt
  1    82.0 Control
  2    97.8 Control
  3    69.9 Control
  4    58.3       A
  ...
  12   48.9       C
 
  > mean(sugar$weight[sugar$trt == "Control"])
  [1] 83.23333
 
あるいは
  > splitplot[1,]
    Block  Insect Mollusc Rabbit    Lime Competition 
  1     A Sprayed   Slugs Fenced Unlimed     Control  
  ...

のようなデータフレームがあるとして，Block が A の項目は
  splitplot[splitplot$Block=="A",]
として得られる。

~/mysplus/crawley.s , p.352

次の例は，term.cos で，行名が Berg である行の各列から，0.999 以上の項目があれば取り出す

  [berg >=  0.999]


# あるマトリックスの，ある名前の行の (1列目の) 数値を見る
  stif.keller.without.svd.Doc3
    [rownames(stif.keller.without.svd.Doc3)=="Gouverneur", 1]











*  図における曲線，座標設定などの詳細 [#j0fca1af]

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.182

なお plot で引数に type="n"を指定すると，点は刻印されない。これに続いて
points() を実行する。

二つのカテゴリについて，変数 Y と 変数 X の対応をプロットする。[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.134







*  aggregate の使い方 [#cbfa9ae1]


[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]   p.43 
例えば
  yield block shade plot
  101   north none  north.none
  108   north none  noth.none

のようなデータがあった時，ブロックとシェードをまとめた平均を作成する場合

  aggregate(yield, by=list(block,shade), mean)
  block shade meanOFyield
  east  none  99

他に  Maindonals   p.225






*  rank 落ち [#s5b1c832]


//  例えば一つの水準に三つの項目があり，どれか一つだけを選択する場合，
//  任意に選んだ二つの水準のそれぞれについて，どちらかが選択されているかどうかが分かれば，
//  残り一つの水準について選択されているかどうかは自動的に分かる。
//
//  つまり一つの水準の値は，残り２水準の値によって計算される。
//  すなわちある水準のベクトルは，残り２水準のベクトルの和で表され，
//  それらに従属している。







*  ニューラルネットワークその２ [#l4142b36]


豊田秀樹「非線形多変量解析」87ページを実行するには

  data(iris)
  library(nnet)

iris.2 # このままだと元のデータフレムの情報として，Species 変数のレベルが三つになっている。
  levels(iris.2$Species)
  #[1] "setosa"     "versicolor" "virginica"

そこで三つ目のレベルを欠損値とする。これによりレベルは二つとして扱われる。
  levels(iris.2$Species)[3]   
  iris.nnet iris.nnet.pre table(iris.2$Species, iris.nnet.pre)
  summary(iris.nnet)








*  決定木でラベルを表示する [#hdf800ef]

text 関数が有効な時と，有効でない時がある。これは，データフレームの作成方法に問題がある。

有効

  senkei senkei.tree plot(senkei.tree,type="u")
  text(senkei.tree)
  summary(senkei.tree)

無効

  base base1 <- base.tree  + man + mit + nicht
     +noch+nun+nur+oder+schon+sich+so+und+zu, data=base1)
  plot(base.tree,type="u")
  text(base.tree)
  summary(base.tree)

ところが，一度
  text(base.tree,labels=col.names(base.tree))
 
  Error in 
 text(xy.coords(x, y, recycle = TRUE), labels, adj, pos, 
   offset,  : 
  	invalid pos value
  In addition: Warning message: 
  NAs introduced by coercion 

を実行 (エラーで無表示となる)後，再度
 text(base.tree)
を行なうと，ラベルが付与された。










*  多重共線性の発見 [#p762e5ab]

朝野 p.108より


説明変数の相関行列 R について R の固有値を求める。このとき固有値の積は，R の行列式に一致する。偏回帰係数の分散は，固有値の逆数の和と一致する。









*  自己組織化マップ [#pffef85e]
[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=R%A4%C7SOM%28%BC%AB%B8%CA%C1%C8%BF%A5%B2%BD%A5%DE%A5%C3%A5%D7%29&word=som]] を参照しました．
  #
  library(som)
  library(class)
  ex gr#ex[1,]
  testplot(test)
# 六角形配置にしたがって円を配置する (下記の 12*8は xdim, ydimに対応)
  symbols(test$grid$pts[,1],test$grid$pts[,2],
    circles=rep(0.5,8*12),inches=F,add=T)
# knn1を使って ex 中のデータが SOMで割り当てた test$code (8*12ある)のどれに近いのかを調査
  bins # exから抜粋した値がどこに割り当てられるのかを示す。 
# 全てをプロットすると収集がつかなくなるので 
  seq(101, 1000, by=25)だけをプロットしてみる
  text(test$grid$pts[bins[seq(101,1000,by=25)],]
    + rnorm(36,0,0.2),col="blue",
     as.character(seq(101,1000,by=25) ))










*  データフレームの操作 [#mf1618ef]

  stifstom4[stifstom4 
  stifstom4[,stifstom4[,c(4:16)]
  stifstom4[,c(4:16)][stifstom4[,c(4:16)]  
  #stifstom4[,c(4:16)][stifstom4 
  stifstom4[stifstom4 > 1] 





*  Dalgaard の本のデータ解析 [#c1ad686c]

research/statistics/Intro.R.R 

*  matrixの作り方 [#n6b7dd23]

例えば自然科学の統計学 p.38 の最小二乗法の計算
  x A B X t(X) %*% X








*  行列式を使った方程式の解 [#l7fd71d4]

自然科学の統計学  p.67
  3x - 2y = 7
  x + 3y = -5 
は R で以下のように解く

  x1 X #X = matrix(c(3,-2,1,3),ncol=2,byrow=T) でももちろん良い。
 
  Y = c(7,-5);Y
  solve(X) %*% Y
  # あるいは
  solve(X,Y)










*  松原望 蟻が黒い確率 [#g145d398]

松原望『実践としての統計学』 p.16, p.163, p.222


B1 = 100% 蟻は黒い, とする。
B2 = 99 % の蟻は黒い, とする。そしてそれぞれの確率を
p(B1) = p(B2) = 0.5 

すると

p(D|B1) = 1, 蟻は100%黒いという仮定が成立して，その時に黒い蟻を見つける確率 。
p(D|B2) = 0.99, 蟻は99%黒いという仮定が成立して，その時に黒い蟻を見つける確率。

このとき p(B1|D) , 黒い蟻を見つけた時に，B1 が正しい確率は，

p(B1|D) = p(B1 * D) / p(D) = p(B1) * p(D|B1) / p(B1) * p(D|B1) + p(B2) * p(D|B2) = 0.5 * 1 / ( 0.5 * 1 + 0.5 * 0.99)

と計算される。








*  作業ディレクトリにあるすべてのオブジェクトの削除 [#u8b09aa4]
  remove(objects())   # S-Plus
  rm(list = ls())     # R






*  図のサイズ調整 [#ve59356c]


  # dev.off()
  # plot.new() 場合によっては描画前に必要となる
  par.old par.old plot(x,y)
  par(par.old)

Dalgaard p. 28, 30, 71
Venables p. 84

これに外周を加えるのが
  oma






*  ブラウザからヘルプを見る [#je7d86ba]

  help.start()








*  R のライブラリ類の場所 [#f343fc64]

  /usr/local/lib/R/library/







*   データテーブルを並び替える stab [#u850d32b]


[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.37 以下のようなデータがあった場合
  > jobs
       BC Alberta Prairies Ontario Quebec Atlantic     Date
  1  1752    1366      982    5239   3196      947 95.00000
  2  1737    1369      981    5233   3205      946 95.08333
  3  1765    1380      984    5212   3191      954 95.16667
 
  ...
 
  Jobs 
によって，次のようなデータに変更できる。
  > Jobs
      values      ind
  1     1752       BC
  2     1737       BC
  3     1765       BC
  ...
 
  142    944 Atlantic
  143    943 Atlantic
  144    942 Atlantic






*  abline と spline, lines の違い [#p65d45d0]

[[旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  Chapter3, p.53, p.54, p.55

abline(x,y) は切片が x 傾きが y の直線を引く

splineは lines と組み合わせて，曲線

  lines(spline(x,y))

他に

  lines(lowess(x,y),) p.35

も使われる。








*   pretty, axis の使い方 [#n74c2833]

pretty は，第一引数をもとに，第二引数で指定された数のベクトルを作成する。

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.59, p.117








*  set.seed [#peb6a616]

set.seed(x) は，R起動には，x をランダムな値に設定して生成される乱数のリストである。
したがって，同じ x の値で実行すれば，同じ乱数のセットが得られるので，
結果として，同じ結果が得られる。[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.64





*  split [#s1d4876e]


第一引数で指定したベクトルを，第二引数で指定したカテゴリーで，リストに分ける。
Maindonald p. 134, p.213










*  text関数の使い方 [#a0c82837]

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.134 より

  plot(weight.split$hb ~ volume.split$hb,
            pch=16, xlim=range(volume), ylim = range(weight),
            ylab="Weight (g)", xlab="Volume (cc)")

  text(weight.split$hb ~ volume.split$hb, labels=c(1:7), 
     adj=c(1.5,-1.5))

ここで adj は(x,y)点に対して x がマイナスなら，左に，また y がマイナスなら上にテキストが付記される。









*  回帰分析で，変数間の相関係数を出すオプション [#f7ea1c56]


[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  # p.137
  allbacks0.lm summary(allbacks.lm, corr=T)

また個々のサンプルを抜いた場合の影響度を調べるには [[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]   p.139
  round(coef(allbacks.lm0), 2)
  round(lm.influence(allbacks.lm0)$coef, 2)








*  変数間の線型性を求める alias [#j5f121ed]

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.166









*  分散分析，回帰分析で base ラインを変更する relevel [#sc3be186]

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]  p.90  p.178  p.210








*   option の使い方 [#e595e701]

Maindonals  p.178






*  散布図の外にプロットを打つ points(, ..., xpd =T) [#cf96f320]

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]   p.202





*   nlme マルチレベルでの固定効果と，ランダム効果 [#w367daae]

[[Maindonald 旧版:http://wwwmaths.anu.edu.au/~johnm/r-book.html]]   p.226
  science.lme summary(science.lme)$tTable







*  カイ二乗検定における有意確率(両側) [#teb7b8ff]

青木先生のサイト
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/one-two.html




*  浜田知久馬 「ロジスティック回帰分析」 [#rebda9b4]

hamada.R を参照  p.141

  x y z # z
  #      x  y
  # [1,] 95  5
  # [2,] 90 10
  #
  fac # fac.z fac.z gl(2,1,2,fac); fac.z
  # fac.z
  # [1] A B
  # Levels: A B
  test.glm glm(z ~ fac.z, family=binomial)
  summary(test.glm)
  summary(test.glm)$coef[1,3]^2 # χ二乗
  summary(test.glm)$coef[2,3]^2 # χ二乗
  #exp(abs(summary(test.glm)$coef[2,1]))
  exp(summary(test.glm)$coef[2,1])
  anova(test.glm, test="Chisq")
  fisher.test(z)
  # p.143
  x y z #> z 
  #     x  y
  #[1,] 5 95
  #[2,] 2 48
  #[3,] 8 42
 fac fac.z #> fac.z
  #[1] A B C
  #Levels: A B C
  test.glm summary(test.glm)
  summary(test.glm)$coef[1,3]^2 # χ二乗
  summary(test.glm)$coef[2,3]^2 # χ二乗
  summary(test.glm)$coef[3,3]^2 # χ二乗
  exp(summary(test.glm)$coef[2,1])
  exp(summary(test.glm)$coef[3,1])
  anova(test.glm, test="Chisq")
  # p.144
  x y z 
  fac 
  test.glm summary(test.glm)
  summary(test.glm)$coef[1,3]^2 # χ二乗
  summary(test.glm)$coef[2,3]^2 # χ二乗
  anova(test.glm, test="Chisq")







*   array の作り方 [#p6acf9ec]

Meyer English Corpus Linguistics p.134

    Type    Form 1-4 Words  5 or more Words 
    PT    Simple NP  216   23 
    PT    Gen. NP  0   0 
    PT    Post.Mod.  1   7 
    APOS    Simple NP  52   10 
    APOS    Gen. NP  18   9 
    APOS    Post.Mod.  12   97 
 
  Meyer array(c(
          216, 0, 1, 
          23, 0, 7,
          52, 18, 12,
          10,  9, 97),
        dim = c(3,2, 2),
        dimnames =
        list(
             Form = c("SimpleNP", "GenNP", "PostMed"), 
                   # = c(1)# = c(1)         
             Length = c("lessWords","moreWords"),      
                   # = c(2)
             Type= c("PT", "AP")     
                   # = c(3)
             ))
  class(Meyer) ftable(Meyer)

c(3, 2, 2) で 3 行 2 列 2 層，入力順に，層(行→列)→層(行→列)と埋める


  > Meyer
  , , Type = PT
 
            Length
  Form       lessWords moreWords
    SimpleNP 216        23      
    GenNP      0         0      
    PostMed    1         7      
 
  , , Type = AP
 
            Length
  Form       lessWords moreWords
    SimpleNP  52        10      
    GenNP     18         9      
    PostMed   12        97      

  > ftable(Meyer)
                     Type  PT  AP
  Form     Length                
  SimpleNP lessWords      216  52
           moreWords       23  10
  GenNP    lessWords        0  18
           moreWords        0   9
  PostMed  lessWords        1  12
           moreWords        7  97
  > 
 
  ####################
  Meyer array(c(
          216, 23,
          0,0,
          1,7,
          52,10,
          18,9,
          12,97),
        dim = c(2,3, 2),
        dimnames =
        list(
             Length = c("Less", "More"),           
             Form = c("SNP", "GNP","PMD"),
             Style= c("PT", "AP")
             ))
c(2, 3, 2) で２行３列，入力順に行→列と埋める　
2*3 を埋め尽くすと，次に次元に移り，2 * 3 を再度埋めていく
  > Meyer
  , , Style = PT
 
        Form
  Length SNP GNP PMD
    Less 216   0   1
    More  23   0   7
 
  , , Style = AP
 
        Form
  Length SNP GNP PMD
    Less  52  18  12
    More  10   9  97
 
  class(Meyer) # "array"
  class(Meyer)<- "table"
  class(Meyer) # "table" 
  アレイをテーブルに変更する
  > ftable(Meyer)
              Style  PT  AP
  Length Form              
  Less   SNP        216  52
         GNP          0  18
         PMD          1  12
  More   SNP         23  10
         GNP          0   9
         PMD          7  97


*  ログリニア分析 [#q5ce84c6]

始めに[[ここ>#matsuda121]]を参照．[[ここ>#loglin2]]と[[ここ>#loglm]]も．[[ここ>#social.mobility2]]も．さらにMeyer.Rを参照。

Meyer.Howa 
ここで list(c(1,2,3)) は list(1:2:3)としても同じで，三次交互作用までを含んだモデル。
list(c(1,3), 2) ならば，1,3 相互の交互作用を含めるが，
2 については，1,3 とは独立の因子とみなす。
list(1,2,3) とすれば，すべてが独立のモデルである。

なお階層モデルで三次の交互作用までを含むということは，2 次，また単独の効果を含む。柳井/緒方『統計学-基礎と応用』,p.150

飽和モデルの結果をみる
# meyer p.135 表，３次の効果を含む飽和モデル，自由度０で説明力は無い

  Meyer.Howa$lrt # 尤度比
  1-pchisq(Meyer.Howa$lrt,Meyer.Howa$df) # 有意度
  Meyer.Howa$pearson# χ二乗値
  1-pchisq(Meyer.Howa$pearson,Meyer.Howa$df) # 有意度
  Meyer.Howa$df # 自由度
  Meyer.Howa$margin # 選択されたモデル。この場合は飽和モデル
 
  #[[1]]
  #[1] "Form"   "Length" "Type"  

# meyer p.135 表，３次の効果を除いた場合
  # K = 2
  Meyer.123<- loglin(Meyer, list(c(1,2),c(2,3),c(1,3)), 
       param=T)
#c (1,2,3) を除いたモデルの適合度は 0.926 ，つまり有意でないので，３次の交互作用を除いたモデルは適合している。
すなわち３次の適合度は省いて構わない。
  Meyer.123$lrt #0.1530465
  1-pchisq(Meyer.123$lrt,Meyer.123$df) #0.9263314
  1-pchisq(Meyer.123$pearson,Meyer.123$df)
  Meyer.123$pearson
  Meyer.123$df # 2
  Meyer.123$margin
  #[[1]]
  #[1] "Form"   "Length"
 
  #[[2]]
  #[1] "Length" "Type"  
 
  #[[3]]
  #[1] "Form" "Type"

# meyer p.135 頁の表, ３次に加えて，さらに２次を全て省いた場合
  # K = 2
  Meyer.1.2.3<- loglin(Meyer, list(1,2,3), param=T)
  Meyer.1.2.3$lrt #  488.0099 
  1-pchisq(Meyer.1.2.3$lrt,Meyer.1.2.3$df) # 0
  1-pchisq(Meyer.1.2.3$pearson,Meyer.1.2.3$df)
  Meyer.1.2.3$pearson
  Meyer.1.2.3$df # 7
  Meyer.1.2.3$margin
  #[[1]]
  #[1] "Form"
 
  #[[2]]
  #[1] "Length"
 
  #[[3]]
  #[1] "Type"





*  loglm によるログリニア，対数線形分析 [#ge2b6b48]
&aname(loglm);
始めに[[ここ>#matsuda121]]を参照．[[ここ>#loglin2]]も．[[ここ>#social.mobility2]]も

loglm 使うためにMASSパッケージを取り込み
  library(MASS)

初めに飽和モデルを分析する。
  Meyer.loglm<- loglm(~Length*Form*Type, data=Meyer)

次に，AIC が最小になるモデルを選び出す。
  # 飽和モデルから，効果を減らしてく
  Meyer.loglm<- loglm(~Length*Form*Type, data=Meyer)
  anova(Meyer.loglm)
  Meyer.loglm2<- update(Meyer.loglm, .~. -Form:Type)
  anova(Meyer.loglm2)
  stepAIC(Meyer.loglm, ~.^2, data=Meyer)
 
  #######
  Start:  AIC= 24 
   ~Length * Form * Type 

                     Df    AIC
  - Length:Form:Type  2 20.153
 #  Length:Form:Type を省いた方が小さい
  <none>                24.000
 
  Step:  AIC= 20.15 
   ~Length + Form + Type + Length:Form + Length:Type + Form:Type 
 
                     Df     AIC
  - Length:Type       1  19.978
 #  Length:Type を省いた方が小さい
  <none>                 20.153
  + Length:Form:Type  2  24.000
  - Length:Form       2 145.147
  - Form:Type         2 153.044
 
  Step:  AIC= 19.98 
   ~Length + Form + Type + Length:Form + Form:Type 
 
                Df     AIC
  <none>            19.978 # これ以上 AIC は小さくならない
  + Length:Type  1  20.153
  - Length:Form  2 255.045
  - Form:Type    2 262.943
  Call:
  loglm(formula = ~Length + Form + Type + 
      Length:Form + Form:Type, 
      data = Meyer, evaluate = FALSE)
 
  Statistics:
                        X^2 df P(> X^2)
  Likelihood Ratio 1.977911  3 0.577004
  Pearson               NaN  3      NaN
  > 
  ##### よって最適なモデルは
  # ~Length + Form + Type + Length:Form + Form:Type
  #####

別データによる分析

それぞれの効果の推定値には，ダミーコーディングによる場合と，ANOVAコーディングによる場合がある。
  rocket<- array(data=c(5,7,8,9,3,21,7,9,6), dim=c(3,3))
  dimnames(rocket)<- list(c("A1","A2","A3"),c("B1","B2","B3"))
  rocket

   # モデル選択
   model0<- loglin(rocket, list(c(1, 2)), param=TRUE)
   model1<- loglin(rocket, list(1, 2), param=TRUE)
   model0
   model1
   p0<- 1-pchisq(model0$lrt, model0$df)
   p1<- 1-pchisq(model1$lrt, model1$df)
   p0
   p1
   AIC0<- model0$pearson-2*model0$df
   AIC1<- model1$pearson-2*model1$df
   AIC0
   AIC1
 
# ここでは飽和モデルを採用し分析してみる
   rocket.df<- as.data.frame.table(rocket)
   colnames(rocket.df)<- c("kozoku","yoko","dosu")
  > loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
  2 iterations: deviation 0 
  $lrt
  [1] 0
 
  $pearson
  [1] 0
 
  $df
  [1] 0
 
  $margin
  $margin[[1]]
  [1] 1 2
 
  $margin[[2]]
  [1] 1
 
  $margin[[3]]
  [1] 2
  $param
  $param$"(Intercept)"
  [1] 1.990005
 
  $param$"1"
           A1          A2          A3 
  -0.07248058 -0.24275578  0.31523636 
 
  $param$"2"
           B1          B2          B3 
  -0.11174159  0.12344831 -0.01170672 
 
  $param$"1.2"
             B1         B2          B3
  A1 -0.1963447  0.1562521  0.04009266
  A2  0.3104027 -0.7720850  0.46168230
  A3 -0.1140580  0.6158330 -0.50177496

# ダミーコーディングによる効果の推定値
   glm0<- glm(dosu ~ kozoku * yoko, data = rocket.df, 
    family = poisson)
   anova(glm0, test = "F")
   summary(glm0)
  #########
  >  summary(glm0)
 
  Call:
  glm(formula = dosu ~ kozoku * yoko, family = poisson, 
     data = rocket.df)
 
  Deviance Residuals: 
  [1]  0  0  0  0  0  0  0  0  0
 
  Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
  (Intercept)      1.60944    0.44721   3.599 0.000320 ***
  kozokuA2         0.33647    0.58554   0.575 0.565538    
  kozokuA3         0.47000    0.57009   0.824 0.409689    
  yokoB2           0.58779    0.55777   1.054 0.291970    
  yokoB3           0.33647    0.58554   0.575 0.565538    
  kozokuA2:yokoB2 -1.43508    0.88730  -1.617 0.105800    
  kozokuA3:yokoB2  0.37729    0.69551   0.542 0.587492    
  kozokuA2:yokoB3 -0.08516    0.77254  -0.110 0.912227    
  kozokuA3:yokoB3 -0.62415    0.79657  -0.784 0.433303    
  ---
  Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
 
  (Dispersion parameter for poisson family taken to be 1)
 
      Null deviance:  2.0874e+01  on 8  degrees of freedom
  Residual deviance: -4.4409e-16  on 0  degrees of freedom
  AIC: 52.681
 
  Number of Fisher Scoring iterations: 3
  #######
   # ANOVA コーディング効果の推定値
  loglm(dosu ~ kozoku * yoko, data = rocket.df, 
       family = poisson)$param
  loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
 
  #######
  > loglm(dosu ~ kozoku * yoko, data = rocket.df, 
       family = poisson)$param
  $"(Intercept)"
  [1] 1.990005
 
  $kozoku
           A1          A2          A3 
  -0.07248058 -0.24275578  0.31523636 
 
  $yoko
           B1          B2          B3 
  -0.11174159  0.12344831 -0.01170672 
 
  $kozoku.yoko
        yoko
  kozoku         B1         B2          B3
      A1 -0.1963447  0.1562521  0.04009266
      A2  0.3104027 -0.7720850  0.46168230
      A3 -0.1140580  0.6158330 -0.50177496
  ########
  > loglin(rocket,margin=list(c(1,2),c(1),c(2)),param=TRUE)
  2 iterations: deviation 0 
  $lrt
  [1] 0
 
  $pearson
  [1] 0
 
  $df
  [1] 0
 
  $margin
  $margin[[1]]
  [1] 1 2
 
  $margin[[2]]
  [1] 1
 
  $margin[[3]]
  [1] 2
 
  $param
  $param$"(Intercept)"
  [1] 1.990005
 
  $param$"1"
           A1          A2          A3 
  -0.07248058 -0.24275578  0.31523636 
 
  $param$"2"
           B1          B2          B3 
  -0.11174159  0.12344831 -0.01170672 
 
  $param$"1.2"
             B1         B2          B3
  A1 -0.1963447  0.1562521  0.04009266
  A2  0.3104027 -0.7720850  0.46168230
  A3 -0.1140580  0.6158330 -0.50177496


なおこのデータ，またコーディングについては 
http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03-01.pdf
の表3.2「ロケットの発射試験」
この例題の飽和モデルによる推定値は 
http://www.sci.kagoshima-u.ac.jp/~ebsa/matsuda01/pdf/ch03-04.pdf
の表3.11にダミーコーディング、表3.12にANOVAコーディング



なお SPSS で出力される Partial association  は「標準化係数（標準効果）」のこと。

以上は[[Rjp Wiki:http://www.okada.jp.org/RWiki/index.php?cmd=read&page=%A3%D1%A1%F5%A3%C1%20%28%BD%E9%B5%E9%BC%D4%A5%B3%A1%BC%A5%B9%29&word=loglin]]を参照しました