R_fromOldHtml1 - RとLinuxと...

RとLinuxと...


R_fromOldHtml1

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

_ [[ ]] の使い方

例えば以下のようなデータで (Maindonald 旧版 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

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

 > 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 を参照

なお boxplot と split, rug を併用した例は Maindonald 旧版 p.225

_ リストの先頭要素を取る

  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下側確率

構文は 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 データの読み込み

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

以下は,excelファイルを読み込むだけでなく,更新なども行なう。

 library(RODBC)
 connExcel 

このあと

 data 

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

 odbcClose(connExcel)

_ 新しいグラフィックデバイス

Maindonald 旧版 , p.22,35 X11() を実行。サイズを指定するには X11(width=7.5,height=4)

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

 dev.off()

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

 dev.off(3)

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

 dev.set(2)

_ 特定の行の取り出し

例えば Maindonald 旧版 , 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]

_ 図における曲線,座標設定などの詳細

Maindonald 旧版 p.182

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

二つのカテゴリについて,変数 Y と 変数 X の対応をプロットする。Maindonald 旧版 p.134

_ aggregate の使い方

Maindonald 旧版 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 落ち

_ ニューラルネットワークその2

豊田秀樹「非線形多変量解析」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)

_ 決定木でラベルを表示する

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)

を行なうと,ラベルが付与された。

_ 多重共線性の発見

朝野 p.108より

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

_ 自己組織化マップ

Rjp Wiki を参照しました.

 #
 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) ))

_ データフレームの操作

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

_ Dalgaard の本のデータ解析

research/statistics/Intro.R.R

_ matrixの作り方

例えば自然科学の統計学 p.38 の最小二乗法の計算

 x A B X t(X) %*% X

_ 行列式を使った方程式の解

自然科学の統計学 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)

_ 松原望 蟻が黒い確率

松原望『実践としての統計学』 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)

と計算される。

_ 作業ディレクトリにあるすべてのオブジェクトの削除

 remove(objects())   # S-Plus
 rm(list = ls())     # R

_ 図のサイズ調整

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

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

これに外周を加えるのが

 oma

_ ブラウザからヘルプを見る

 help.start()

_ R のライブラリ類の場所

 /usr/local/lib/R/library/

_ データテーブルを並び替える stab

Maindonald 旧版 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 の違い

旧版 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 の使い方

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

Maindonald 旧版 p.59, p.117

_ set.seed

set.seed(x) は,R起動には,x をランダムな値に設定して生成される乱数のリストである。 したがって,同じ x の値で実行すれば,同じ乱数のセットが得られるので, 結果として,同じ結果が得られる。Maindonald 旧版 p.64

_ split

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

_ text関数の使い方

Maindonald 旧版 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 がマイナスなら上にテキストが付記される。

_ 回帰分析で,変数間の相関係数を出すオプション

Maindonald 旧版 # p.137

 allbacks0.lm summary(allbacks.lm, corr=T)

また個々のサンプルを抜いた場合の影響度を調べるには Maindonald 旧版 p.139

 round(coef(allbacks.lm0), 2)
 round(lm.influence(allbacks.lm0)$coef, 2)

_ 変数間の線型性を求める alias

Maindonald 旧版 p.166

_ 分散分析,回帰分析で base ラインを変更する relevel

Maindonald 旧版 p.90 p.178 p.210

_ option の使い方

Maindonals p.178

_ 散布図の外にプロットを打つ points(, ..., xpd =T)

Maindonald 旧版 p.202

_ nlme マルチレベルでの固定効果と,ランダム効果

Maindonald 旧版 p.226

 science.lme summary(science.lme)$tTable

_ カイ二乗検定における有意確率(両側)

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

_ 浜田知久馬 「ロジスティック回帰分析」

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 の作り方

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 を埋め尽くすと,次に次元に移り,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

_ ログリニア分析

始めにここを参照.ここここも.[[ここ>#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 表,3次の効果を含む飽和モデル,自由度0で説明力は無い

 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 表,3次の効果を除いた場合

 # K = 2
 Meyer.123<- loglin(Meyer, list(c(1,2),c(2,3),c(1,3)), 
      param=T)

#c (1,2,3) を除いたモデルの適合度は 0.926 ,つまり有意でないので,3次の交互作用を除いたモデルは適合している。 すなわち3次の適合度は省いて構わない。

 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 頁の表, 3次に加えて,さらに2次を全て省いた場合

 # 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 によるログリニア,対数線形分析

始めにここを参照.ここも.[[ここ>#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を参照しました

 
Link: Rの備忘録(1823d) R_fromOldHtml1_2(4099d)
Last-modified: 2007-09-25 (火) 17:38:06 (4099d)