きわめて個人的なR言語のメモ

2005-01-22Rで対数線形モデル

近日中に全面修正します。

家事に対する意見

 放送大学の教材「社会調査法」(1992年初版)に載っている対数線形モデル?の例題を、に搭載されている対数線形モデルのための関数loglinを使って分析してみる。

 「「家事は女の仕事である」という意見に対する賛否」を「結婚時の就職状況」と「現在の就職状況」から分類したもの。

 既存の多重クロス集計表をオブジェクト"shigoto"に入力する。オブジェクトの名前は何でもいい。

shigoto <- array(data=c(7,5,14,23,10,1,16,11),dim=c(2,2,2))

 各カテゴリに名前を付ける。日本語対応版のRなら日本語のカテゴリ名を付けられます。

dimnames(shigoto) <- list(c("結婚時有職","結婚時無職"),c("現在有職","現在無職"),c("賛成","反対"))

 "shigoto"の内容を表示させる。

shigoto

 見やすい形の多重クロス集計表に整形して表示させる。

ftable(shigoto)

いよいよloglin

さまざまなモデルでログリニア分析を行い、モデルを選択する。

ここで、liglin関数の二つ目の引数でモデルを指定するのですが、この例題の場合、1=「結婚時の就職状況(有職/無職)」、2=「現在の就職状況(有職/無職)」、3=「「家事は女の仕事」に対する賛否(賛成/反対)」を示しています。また、loglin関数では階層的なモデルを想定しており、高次のモデルには下位のモデルが含まれます。つまり、 list(c(A,B)) は A*B という二次の交互作用と、AとBそれぞれの主効果の両方を含んでおり、模式的に表すと【A+B+A*B+e】という意味になります。

飽和モデル

 三変数が全て関連しているモデル。

 【A+B+C+A*B+B*C+C*A+A*B*C+e】

model0 <- loglin(shigoto, list(c(1,2,3)), param=TRUE)

対連関モデル(均一連関モデル)

 飽和モデルから三次の交互作用を取り除いたモデル。

 【A+B+C+A*B+B*C+C*A+e】

model1 <- loglin(shigoto, list(c(1, 2), c(1, 3), c(2, 3)), param=TRUE)

条件付独立モデル

 対連関モデルから特定の二変数の交互作用を取り除いたモデル。

 【A+B+C+A*B+B*C+e】

model2 <- loglin(shigoto, list(c(1, 2), c(1, 3)), param=TRUE)
model3 <- loglin(shigoto, list(c(1, 2), c(2, 3)), param=TRUE)
model4 <- loglin(shigoto, list(c(1, 3), c(2, 3)), param=TRUE)

一変数独立モデル

 変数Cが変数Aと変数Bの組み合わせから独立であるモデル。

 【A+B+C+A*B+e】

model5 <- loglin(shigoto, list(c(1, 2),3), param=TRUE)
model6 <- loglin(shigoto, list(c(1, 3),2), param=TRUE)
model7 <- loglin(shigoto, list(c(2, 3),1), param=TRUE)

三変数独立モデル

 三つの変数がそれぞれ独立しているモデル。(χ二乗検定が想定しているモデル。)

 【A+B+C+e】

model8 <- loglin(shigoto, list(1, 2, 3), param=TRUE)

分析結果を表示する。

 出力のうち、$lrtは尤度比、$pearsonはχ二乗値、$dfは自由度、$marginは出力されたモデル。

model0
model1
model2
model3
model4
model5
model6
model7
model8

 各モデルにおける観測度数と推定度数の適合度の検定の有意水準(p値)を計算する。

p0 <- 1-pchisq(model0$lrt, model0$df)
p1 <- 1-pchisq(model1$lrt, model1$df)
p2 <- 1-pchisq(model2$lrt, model2$df)
p3 <- 1-pchisq(model3$lrt, model3$df)
p4 <- 1-pchisq(model4$lrt, model4$df)
p5 <- 1-pchisq(model5$lrt, model5$df)
p6 <- 1-pchisq(model6$lrt, model6$df)
p7 <- 1-pchisq(model7$lrt, model7$df)
p8 <- 1-pchisq(model8$lrt, model8$df)

 計算したp値を表示させます。

 0.05以下のモデルは当てはまりが悪いということになる(んだろうと思う)。

p0
p1
p2
p3
p4
p5
p6
p7
p8

 堀先生のページを読むと、AIC=χ2-2*df ということなので、これも算出してみます。AICが最も小さいモデルを選択することになります。ただし、「2.0以内の違いは意味のある違いと解釈しない方がいい」ということです。

AIC0 <- model0$pearson-2*model0$df
AIC1 <- model1$pearson-2*model1$df
AIC2 <- model2$pearson-2*model2$df
AIC3 <- model3$pearson-2*model3$df
AIC4 <- model4$pearson-2*model4$df
AIC5 <- model5$pearson-2*model5$df
AIC6 <- model6$pearson-2*model6$df
AIC7 <- model7$pearson-2*model7$df
AIC8 <- model8$pearson-2*model8$df

 計算したAICを表示させます。

AIC0
AIC1
AIC2
AIC3
AIC4
AIC5
AIC6
AIC7
AIC8

 ついでにモザイク図を作成し、指定したモデルでの残差で色を付ける。

 コピペすると一気にグラフが表示されますが、「PageUp」と「PageDown」キーで前後のグラフを見ることが出来ます。見られない場合は、グラフのウィンドウのメニューからHistory→Recordingと選択して、もう一度コマンドを入れてください。

mosaicplot(shigoto, main="model0飽和モデル", margin=list(c(1,2,3)), shade=TRUE)
mosaicplot(shigoto, main="model1対連関モデル", margin=list(c(1, 2), c(2, 3), c(1, 3)), shade=TRUE)
mosaicplot(shigoto, main="model2条件付独立モデル", margin=list(c(1, 2), c(1, 3)), shade=TRUE)
mosaicplot(shigoto, main="model3条件付独立モデル", margin=list(c(1, 2), c(2, 3)), shade=TRUE)
mosaicplot(shigoto, main="model4条件付独立モデル", margin=list(c(1, 3), c(2, 3)), shade=TRUE)
mosaicplot(shigoto, main="model5一変数独立モデル", margin=list(c(1, 2),3), shade=TRUE)
mosaicplot(shigoto, main="model6一変数独立モデル", margin=list(c(1, 3),2), shade=TRUE)
mosaicplot(shigoto, main="model7一変数独立モデル", margin=list(c(2, 3),1), shade=TRUE)
mosaicplot(shigoto, main="model8三変数独立モデル", margin=list(1, 2, 3), shade=TRUE)

最も当てはまりが良さそうなmodel2のモザイク図がこれ。残差が小さいので白いです。

f:id:bob3:20050123222855:image

#ここまでは出来ましたが、結果をどう解釈すればよいのか解りません(^^;;;

トラックバック - http://r-user.g.hatena.ne.jp/bob3/20050122

2004-09-29Rで因子分析その6 最尤法で因子分析

今度は最尤法の斜交回転(プロマックス)で因子分析をしてみましょう。

主因子法のときと同じようにデータを読み込みます。

baseball <- read.csv("C:/R/data/data.csv")
attach (baseball)
baseball[[1]] <- NULL

次に、青木先生作られた最尤法での因子分析用関数“factanal2”を使います。

ここでは因子の数を3個に決めて分析します。また、因子得点も回帰法で求めます。

factanal2(baseball,factors=3,scores="regression")

アウトプットは整形されて出てきます。

[1] H0: 3 factors are sufficient.
   Chi sq.       d.f.    P value 
35.0058400 25.0000000  0.0880976 
[1] Factor loadings(rotation:promax)
                Factor1     Factor2    Factor3 Communality
試合数       0.08692896  0.61345987  0.1618269  0.41007759
安打         0.13100079  0.79228560  0.3577115  0.77283520
二塁打      -0.10197822  0.29937034  0.8626870  0.84425105
三塁打      -0.42225483  0.26920935 -0.1534471  0.27431883
本塁打       0.92100110  0.03529310 -0.1251296  0.86514604
三振         0.59384076  0.24216411 -0.4478380  0.61184920
四球         0.47323091  0.04773030 -0.1017807  0.23658499
死球         0.31730815  0.05930770  0.2077376  0.14735678
犠打        -0.66492790 -0.01359067 -0.1065619  0.45366926
犠飛        -0.01606036  0.06812931  0.2237036  0.05494284
盗塁        -0.44385786  0.58826249 -0.1466628  0.56457252
SS.loadings  2.37833522  1.58270947  1.2745596  5.23560430
Proportion  21.62122926 14.38826790 11.5869056 47.59640275
Cum.Prop.   21.62122926 36.00949715 47.5964028          NA
[1] Factor scores(regression)
    
         Factor1     Factor2     Factor3
  1   0.82284220  1.19656161 -0.20089422
  2   0.89770074  0.37227000  0.61553463
  3  -0.02846898 -0.57487912 -0.98405791
(以下省略)

んー、p値が0.088と0.05を越えています。モデルとしていまいちのようですね。

(あ、因子間の相関係数は出力されないのか……)

特定の因子に対する負荷量や共通性から見て「犠飛」は分析から除外した方が良さそうです。

baseball[[10]] <- NULL

として、「犠飛」を削除しました。

続けて再度分析をしましょう。

factanal2(baseball,factors=3,scores="regression")

結果はこうなります。

[1] H0: 3 factors are sufficient.
    Chi sq.        d.f.     P value 
31.48900855 18.00000000  0.02525308 
[1] Factor loadings(rotation:promax)
               Factor1     Factor2      Factor3 Communality
試合数       0.1735365  0.59466957  0.126250862   0.3996861
安打         0.1252477  0.87842045  0.035327756   0.7885575
二塁打      -0.1964831  0.53849450 -0.543264388   0.6237181
三塁打       0.5879066  0.11238270  0.009059499   0.3583461
本塁打      -0.6823296  0.08378079  0.535535746   0.7593914
三振        -0.1207402  0.09961408  0.706445446   0.5235663
四球        -0.3057383  0.02804047  0.337839583   0.2083978
死球        -0.3496064  0.15590027  0.005603665   0.1465609
犠打         0.6102822 -0.13434296 -0.234661543   0.4455584
盗塁         0.7703390  0.39772686  0.090372922   0.7597761
SS.loadings  2.0917607  1.64612862  1.275669593   5.0135589
Proportion  20.9176069 16.46128620 12.756695933  50.1355891
Cum.Prop.   20.9176069 37.37889313 50.135589060          NA
[1] Factor scores(regression)
    
         Factor1     Factor2      Factor3
  1   0.03886573  1.28682950  0.784845204
  2  -0.89193362  0.59099874  0.050250995
  3   0.20394105 -0.70609811  0.418533208
(以下省略)

p値が0.025と0.05を切っています。悪くないモデルのようです。

因子負荷量を書き直すと

因子1 因子2 因子3

盗塁 0.77 0.40 0.09

本塁打 -0.68 0.08 0.54

犠打 0.61 -0.13 -0.23

三塁打 0.59 0.11 0.01

死球 -0.35 0.16 0.01

安打 0.13 0.88 0.04

試合数 0.17 0.59 0.13

三振 -0.12 0.10 0.71

二塁打 -0.20 0.54 -0.54

四球 -0.31 0.03 0.34

この因子の名前を付けるとすると、因子1は盗塁が多くホームランは少ない「ランナー因子」、因子2はヒットを多く打つ「アベレージヒッター因子」、因子3はホームランが多い反面三振は多い「パワーヒッター因子」といったところでしょうか。

通行人通行人2005/07/01 15:16因子数のP値は,帰無仮説が「その因子数で十分」ということであることに注意。
試合数を一緒に分析することは問題があるのでは?また,試合数で割った値を分析する方がよいのではないでしょうか。

トラックバック - http://r-user.g.hatena.ne.jp/bob3/20040929

2004-09-28Rで因子分析その5 変数選択をして主因子法で

変数選択

では、「二塁打」「死球」「犠飛」の変数をデータから削除してしまいます。

baseball[[3]] <- NULL
baseball[[7]] <- NULL
baseball[[8]] <- NULL

これで、「二塁打」「犠飛」「死球」のデータが取り除かれました。

(ここはもっとスマートなやり方がありそうに思うので、ご存知の方がいらっしゃったら教えてください。)

再度分析

では、もう一度分析します。

result <- pfa(baseball)
fl.matrix.out(result)

アウトプットは

After Varimax rotation
        Fac001 Fac002 Communality
 Var001 -0.064 -0.642  0.416
 Var002 -0.088 -0.672  0.460
 Var003  0.520 -0.230  0.324
 Var004 -0.873 -0.105  0.773
 Var005 -0.466 -0.317  0.317
 Var006 -0.434 -0.100  0.198
 Var007  0.672  0.016  0.452
 Var008  0.547 -0.555  0.607
 Sq.Sum  2.200  1.347
  Cont.   27.5   16.8
   Cum.   27.5   44.3

書き直すと、

因子1 因子2 共通性

------------------------------

試合数 -0.06 -0.64 0.42

安打 -0.09 -0.67 0.46

三塁打 0.52 -0.23 0.32

本塁打 -0.87 -0.11 0.77

三振 -0.47 -0.32 0.32

四球 -0.43 -0.10 0.20

犠打 0.67 0.02 0.45

盗塁 0.55 -0.56 0.61

-----------------------------

固有値 2.20 1.35

寄与率 27.5 16.8

累積 27.5 44.3

うーん、マイナスの因子負荷量が多いので解釈が難しいですが……

バッターのタイプには3種類存在するなどを参考にして考えると、因子1はホームランよりも三振を少なくしボールにバットを上手く当てることを重視する「アベレージヒッター因子」、因子2はホームランを多く打ち、フルスイングをするので三振は多くなる「パワーヒッター因子」といったところでしょうか。

なお、因子得点は

result$scores

とすれば出てきます。

もとデータの選手の名前に対応させて散布図を描けば、選手のタイプ分けも出来るでしょう。もちろん、クラスタ分析にかけることも出来ます。

(本当はRでもそこまで出来るはずですが、私の勉強不足でやり方がよくわかんないです。すいません)

トラックバック - http://r-user.g.hatena.ne.jp/bob3/20040928

2004-09-27Rで因子分析その4 主因子法で因子分析

主因子法の直交回転

まずは青木先生の関数“pfa”を使って、主因子法の直交回転(バリマックス)で因子分析をしてみる。

result <- pfa(baseball)
fl.matrix.out(result)

これで分析結果が整形されて表示されます。

After Varimax rotation
        Fac001 Fac002 Fac003 Communality
 Var001  0.087 -0.610 -0.195  0.418
 Var002  0.098 -0.686 -0.392  0.633
 Var003 -0.025 -0.129 -0.855  0.748
 Var004 -0.493 -0.272  0.158  0.342
 Var005  0.898 -0.064  0.022  0.811
 Var006  0.540 -0.377  0.383  0.581
 Var007  0.427 -0.084  0.036  0.191
 Var008  0.280 -0.018 -0.254  0.143
 Var009 -0.663 -0.029  0.162  0.467
 Var010  0.029 -0.057 -0.206  0.047
 Var011 -0.516 -0.602  0.151  0.652
 Sq.Sum  2.327  1.454  1.252
  Cont.   21.2   13.2   11.4
   Cum.   21.2   34.4   45.8

分かりやすく書き直すとこんな感じでしょうか。

【回転後の因子負荷量】

項目 因子1 因子2 因子3 共通性
本塁打 0.90 -0.06 0.02 0.81
犠打 -0.66 -0.03 0.16 0.47
三振 0.54 -0.38 0.38 0.58
三塁打 -0.49 -0.27 0.16 0.34
四球 0.43 -0.08 0.04 0.19
死球 0.28 -0.02 -0.25 0.14
安打 0.10 -0.69 -0.39 0.63
試合数 0.09 -0.61 -0.20 0.42
盗塁 -0.52 -0.60 0.15 0.65
二塁打 -0.03 -0.13 -0.86 0.75
犠飛 0.03 -0.06 -0.21 0.05
固有値 2.33 1.45 1.25
寄与率 21.2 13.2 11.4
累積 21.2 34.4 45.8

特定の因子への負荷量や共通性から見て「犠飛」「死球」は除外した方が良さそうです。また、第3因子に強い負荷量を持つのは「二塁打」だけですので、これも除外して、因子を二つに決めて分析し直した方が良さそうです。

トラックバック - http://r-user.g.hatena.ne.jp/bob3/20040927

2004-09-26Rで因子分析その3 二変量の分析

回帰直線入り散布図行列

二編量の関係を吟味してみる。

散布図行列を描く

普通に散布図行列を描くには

pairs(baseball)

で良いですが、ここではRjpWikiで公開されている

pairs用 回帰直線つきpanelを使わせて頂いて、回帰直線も表示したいと思います。

panel.lsfit <- function(x,y,...)
{
f <- lsfit(x,y)$coef
xx<- c ( min (x) , max(x) )
yy <- f["X"] * xx+ f["Intercept"]
lines(xx,yy,col="red")
points(x,y)
}
pairs(baseball, panel = panel.lsfit, pch = ".") 

相関行列を算出

ここでは青木先生の関数“my.cor”を使います。

my.cor(baseball)

で、上三角行列が「単相関係数」、下三角行列が「偏相関係数」、対角要素が「重相関係数」である相関行列が得られる。

            Var 1       Var 2       Var 3       Var 4       Var 5        Var 6
Var 1  0.49835624  0.40778275  0.04897438  0.12344932  0.10101434  0.013764419
Var 2  0.46190776  0.60631408 -0.07673645 -0.04679785 -0.30636509 -0.008015155
Var 3  0.15587627 -0.12310150  0.42688692 -0.33025436 -0.14101041 -0.209701240
Var 4  0.07122171 -0.08357778 -0.11436855  0.78567585  0.48191262  0.474239488
Var 5  0.21238691 -0.34848652 -0.02442644  0.35389118  0.61150159  0.259878898
Var 6 -0.05639442  0.04290625 -0.06561933  0.41223813  0.02322316  0.538118607
Var 7 -0.01568597  0.13425732  0.01626393  0.39812289 -0.18950813 -0.174814423
Var 8 -0.02568533 -0.13160480  0.21565671 -0.42459077 -0.05901787  0.070749878
Var 9 -0.04652111  0.19419433  0.03019263 -0.10770892  0.03879625  0.205311341
            Var 7      Var 8        Var 9
Var 1  0.09508995 -0.1386073  0.043771291
Var 2  0.19210838 -0.1167896  0.217348290
Var 3 -0.12014682  0.3644515 -0.009790575
Var 4  0.36857085 -0.5891611 -0.016409118
Var 5 -0.03273129 -0.2842107 -0.045642082
Var 6  0.03870287 -0.2441358  0.167593265
Var 7  0.48812646 -0.2525798  0.054702262
Var 8 -0.01260427  0.6341374 -0.032818085
Var 9  0.05793332 -0.0293333  0.303570602

トラックバック - http://r-user.g.hatena.ne.jp/bob3/20040926