Theoretical Sociology

太郎丸博のブログです。研究ノートや雑感などを掲載しています。(このページは太郎丸が自主的に運営しています。京都大学の公式ページではありません。)
<< October 2019 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 >>
 
RECOMMEND
後期近代と価値意識の変容: 日本人の意識 1973-2008
後期近代と価値意識の変容: 日本人の意識 1973-2008 (JUGEMレビュー »)

NHKの日本人の意識調査のデータをつっこんで分析した本です。
RECOMMEND
Labor Markets, Gender and Social Stratification in East Asia: A Global Perspective (The Intimate and the Public in Asian and Global Perspectives)
Labor Markets, Gender and Social Stratification in East Asia: A Global Perspective (The Intimate and the Public in Asian and Global Perspectives) (JUGEMレビュー »)

直下の和書の英語版です。審査を通過するためにレフェリーのコメントに従って若干修正してあります。
RECOMMEND
東アジアの労働市場と社会階層 (変容する親密圏/公共圏)
東アジアの労働市場と社会階層 (変容する親密圏/公共圏) (JUGEMレビュー »)

GCOEの成果をまとめた本です。日本を中心に韓国、台湾(中国も少し)との比較研究をしてます。
RECOMMEND
若年非正規雇用の社会学‐階層・ジェンダー・グローバル化 (大阪大学新世紀レクチャー)
若年非正規雇用の社会学‐階層・ジェンダー・グローバル化 (大阪大学新世紀レクチャー) (JUGEMレビュー »)
太郎丸 博
拙著です。非正規雇用に関する本はたくさんありますが、「なぜ正規雇用と非正規雇用では賃金格差があるのか」など当たり前と思われがちな問題を突き詰めて考えてみました。
RECOMMEND
フリーターとニートの社会学
フリーターとニートの社会学 (JUGEMレビュー »)

拙編です。オーソドックスな計量社会学の手法で、若年非正規雇用や無職にアプローチした本です。白い装丁なので、輪郭がわからないですね...
RECOMMEND
人文・社会科学のためのカテゴリカル・データ解析入門
人文・社会科学のためのカテゴリカル・データ解析入門 (JUGEMレビュー »)
太郎丸 博
拙著です。軽く読み流すのは難しいですが、まじめに一歩一歩勉強するために作りました。
ARCHIVES
RECENT COMMENT
  • 阪大を去るにあたって: 社会学の危機と希望
    charlestonblue (10/08)
  • Cohen et. al 2011 「フェミニズムの方法論的インパクト: 社会学のやっかいな問題?」
    abe daijyu (10/05)
  • アマチュア社会学の可能性
    読者 (02/20)
  • 社会システム理論の野望、あるいは全体性へのオブセッション
    宮国 (12/19)
  • 片山他 2015「図書館は格差解消に役立っているのか?」
    オカベ (12/09)
  • ランダム効果の意味、マルチレベル・モデル、全数調査データ分析
    YZ (12/07)
  • 学歴社会から「学習資本」社会へ:日本の教育と社会における階級形成の再編
    赤尾勝己 (02/11)
  • グラフィカル・モデリングとは?
    anonymous (11/30)
  • Rスクリプト覚書き:vglm関数で平行性の仮定を置かずに順序ロジット
    ほっくー (08/05)
  • 台湾の経済: 典型NIESの光と影
    おーまきちまき (07/19)
RECENT TRACKBACK
 
Rスクリプト覚え書き: latticeパッケージのxyplotで細かくグループ分けしてプロット
lattice でいろいろオプションをつけて、細かくグループわけして、好みの書式に整えてプロットする方法がわかったので、メモしておく。もっとかんたんに書けるかもしれないが、これが現在の私の知識の限界ということで。
# car パッケージの賃金に関するデータを読み込む
library(car)
head(SLID)
summary(SLID)

library(lattice)

# グループ別に年齢と賃金の散布図を作る
xyplot(wages ~ age | sex + language, data=SLID,  # 性別と言語で6つのパネルにサンプルを分割
       groups=education>12,  # プロットの点を教育年数が12年以上かどうかで色分け
       auto.key=T) # key は凡例をさすようで、凡例を作るように指定

次は、散布図ではなく、平均値を計算して折れ線グラフにする。
# 60才未満のデータに関してグループ別に年齢と賃金の平均値をプロット
# まず平均値のデータを作る
attach(SLID[SLID$age<60, ])
(d0 <- tapply(wages, list(age, sex, language, education>12), mean, na.rm=T)) # グループ別に賃金の平均を計算
detach(SLID[SLID$age<60, ])

d1 <- as.data.frame.table(d0) # 賃金の平均のデータ (array) をデータ・フレームに変換
head(d1)
names(d1) <- c("age", "sex", "language", "education", "mean.wages") # 変数名をわかりやすく
levels(d1$education) <- c("Primary or Secondary Education", "Tertiary Education") # これもわかりやすく
xyplot(mean.wages ~ age | language + sex, data=d1,  type="b",
       groups=education,  
       auto.key=T)

さらに見やすくなるように細かく作りこんでみる。
# 凡例に線も書き加えて、年齢を5歳刻みで表示させる
xyplot(mean.wages ~ age | language + sex, data=d1,  type="b",
       groups=education,  
       xlab=list(cex=1.3), # ラベルを大きく
       ylab=list("mean of wage", cex=1.3) , # Y軸のラベル
       auto.key=list(lines=T, cex=1.2), # 凡例に線を加え、文字を少し大きくする。 詳細は ?simpleKey で調べられる
       scales=list(x=list(at=5 + 0:8*5, labels=paste(20 + 0:8*5)),  # 目盛りの位置とラベルを指定
                  cex=1.3), # 目盛りのラベルを大きく       
)

モノクロで印刷する場合はいろいろ設定が必要
# グレースケールで印刷する場合に、線の種類などを細かく変更
xyplot(mean.wages ~ age | language + sex, data=d1,  type="l",
       groups=education, 
       ylab="mean of wage",
       scales=list(x=list(at=5 + 0:8*5, labels=paste(20 + 0:8*5))),
       strip=strip.custom(which.given=1, bg=gray(.8)), # ストリップの背景を薄めのグレーに
       col="black", lty=1:2, # col で線の色を指定
       key=list(space="bottom", #凡例の位置指定。ほかに top, left, right 
                col="black", # 凡例の色
                columns = 2, # 凡例の列の数
                between = 0.8, # 凡例の記号と文字の間隔
                text=list(levels(d1$education)), # 凡例の文字
                lines=list(lty=1:2, type="l")) # 凡例の線種など 
)
Rスクリプト覚え書き:semPlotパッケージで確証的因子分析のパス図を描く
パス図を描くのがけっこう面倒なので、semPlot パッケージを使ってみた。細かい設定は、qgraph パッケージのヘルプに記載してあるものを指定すればよい。semPlot は ggraph の機能を使ったパッケージのようで、qgraph のオプションが概ね使えるようである。
# 架空のデータを人工的に作る
N <- 100

latent <- matrix(rnorm(N*8), N, 8)
latent[,2] <- latent[,1] + latent[,2]
colnames(latent) <- c("f1", "f2", paste("e", 1:6, sep=""))
#latent <- as.data.frame(latent)
(b1 <- matrix(0, 6, 2) )
b1[1:3, 1] <- c(.4, .6, .8)
b1[4:6, 2] <- c(.5, .7, .9)
x <- b1 %*% t(latent[,1:2]) + matrix(rnorm(6*N, sd=0.7), 6, N)
d0 <- as.data.frame(t(x))
head(d0)

# ここから確証的因子分析

model1 <-   'f1 =~ V1 + V2 + V3
             f2 =~ V4 + V5 + V6
              V3 ~~ V6'
library(lavaan)
cfa1 <- cfa(model1, data= d0)
summary(cfa1)
inspect(cfa1, "fit")
inspect(cfa1, "modindices")

# ここからパス図の描画
library(semPlot)
semPaths(cfa1, "par") # ほとんどデフォルトの図(下の最初の図)
semPaths(cfa1, "par", style="lisrel", # LISREL型の図(下の二番目の図)
 mar=c(6,1,3,1), # 余白の指定、下、左、上、右の順
 edge.label.cex=.8, # 矢印の係数の文字の大きさ
 fade=F, # デフォルトでは係数の大きさに比例して、係数や矢印の色が薄くなるので、色の濃さを均一に指定
 gray=T) # モノクロにしたいときに指定

Rスクリプト覚書き:選好度の写像分析法

現実逃避のために選好度の写像分析法のプログラムを書いたので、掲載しておく。選好度の写像分析法とは、多次元空間上に外部データの情報を重ねるために用いられる。私の場合、職業威信のデータを多次元尺度法 (MDS) にかけて、二次元空間上に職業をプロットしたあと、職業の特性(在職者の平均収入や学歴など)をプロットするという作業をするために使ったことがある。このように MDS の結果を解釈するときに役に立つが、別の利用法もあるのかもしれない。下の図は米国の各州を犯罪の発生率のデータを使って犯罪発生率が類似しているほど近くになるように各州をプロットしたものである。

例えば左上の方に Georgia と Louisiana が近くに布置されているが、これらの州は犯罪発生率のパターンが似ているということである。矢印と赤字でしめしてあるのが、選好度の写像分析法の結果で Murder の矢印が左上に伸びているが、こちらの方向にあるほど殺人の発生率が高いことを示す。正確にはこの矢印上の座標の値が殺人の発生率とできるだけ高く相関するように矢印の角度を決めている。例えば、North Carolina では殺人が多く、Hawaii では逆に少ない、ということになる。矢印の長さに意味はないが、相関の大きさに比例するように作ることもできよう。

データは必ずしも選好度でなくても良いのだが、岡太・今泉 (1994) では選好度の分析が想定されていたので、このように呼ばれたと思われる。以下は R のスクリプト。我流なのでちょっと恥ずかしい。 prefmap という関数が選好度の写像分析のための関数で、s が選好度行列(例では各州の犯罪発生率と都市人口比率)で x は各州の座標である。結果は、選好度をベクトルで表す場合と点で表す場合の二通りで、適合度の指標として調整決定係数を示してある。適合度とは、点やベクトルがどの程度うまく実際の犯罪発生率や都市度を示しているかを示す指標である。この例の場合、都市度の調整決定係数が低く、この空間上に点やベクトルで都市度を示すのが難しいことがわかる。

自作スクリプトの出力
$ベクトル
             x[, 1]     x[, 2]
Murder   -0.7695444  0.6385933
Assault  -0.9386242  0.3449415
UrbanPop -0.3401471 -0.9403722
Rape     -0.5672143 -0.8235703

$点の座標
              x[, 1]     x[, 2]
Murder   -23.1943939  18.916051
Assault   36.1361960 -13.511878
UrbanPop  -0.8182509  -1.820102
Rape     -29.1844961 -42.215030

$調整決定係数
                  adj.r.squared.vector adj.r.squared.points
Response Murder              0.9176126            0.9166173
Response Assault             0.8881425            0.8860540
Response UrbanPop            0.1492228            0.2326966
Response Rape                0.9863168            0.9864470

上記の計算のためのスクリプト(2次元空間にしか対応していないので注意)
prefmap <- function(s, x){ # s が選好度行列(行が対象)、x が布置行列(2行)
  x.squared <- rowSums(x^2)
  s <- as.matrix(s)
  lm.vector <- lm(s ~ x[,1] + x[,2])
  lm.points <- lm(s ~ x.squared + x[,1] + x[,2]) # まとめて回帰分析
  ivector0 <- coef(lm.vector)[2:3, ]
  ivector <- t(ivector0)/sqrt(colSums(ivector0^2))
  ipoints <- 0.5*t(coef(lm.points)[3:4, ])/coef(lm.points)[2, ]
  getR2 <- function(x) x$adj.r.squared # lmサマリー・オブジェクトから調整決定係数をとってくる関数 
  adj.r.squared.vector <- sapply(summary(lm.vector), getR2)
  adj.r.squared.points <- sapply(summary(lm.points), getR2)
  return(list("ベクトル"=ivector, "点の座標"=ipoints,
              "調整決定係数"=cbind(adj.r.squared.vector, adj.r.squared.points)))
}
# 分析例
head(USArrests) # 米国各州の各種犯罪発生度と都市化の程度のデータ
d0 <- scale(USArrests[,-3]) #データを標準化
dist0 <- dist(d0) # 距離行列の作成
library(MASS)
mds0 <- isoMDS(dist0, k=2) # 非計量的多次元尺度法、2次元解

plot(mds0$points, type="n") # 結果をプロット
text(mds0$points, labels=labels(dist0))

(p0 <- prefmap(USArrests, mds0$points)) # 理想点と理想ベクトルの計算 
arrows(0, 0, p0$ベクトル[,1], p0$ベクトル[,2]) # 点は書きにくいのでベクトルを矢印で
text( p0$ベクトル[,1], p0$ベクトル[,2], names(USArrests), col="red")
Rスクリプト覚え書き: deviance 関数が lm に対して返すのは deviance ではない?

今日はじめて気づいたのだが、R で lm の結果から deviance を取り出そうとすると、"-2 * 対数尤度" ではなく、Residual Sum of Square (RSS) を返すのである。このあたりの関係についての知識があいまいだったので、以下では簡単にまとめておく。 逸脱度 (deviance) とは、一般的には入れ子の関係にある複数のモデルの当てはまりのよさを比較するための統計量であるが、ふつうは、あるモデルと飽和モデルのあてはまりのよさを比較したものであり、

逸脱度 = −2 × 対数尤度  (1)
で定義される。glm を使ってロジスティック回帰分析をした結果から対数尤度と逸脱度を計算すると、以下のように確かに (1)式の関係が成り立っているのがわかる。
> library(AER) # AER パッケージの SwissLabor というデータを使う
> data(SwissLabor) 
> head(SwissLabor)
  participation   income age education youngkids oldkids foreign
1            no 10.78750 3.0         8         1       1      no
2           yes 10.52425 4.5         8         0       1      no
3            no 10.96858 4.6         9         0       0      no
4            no 11.10500 3.1        11         2       0      no
5            no 11.10847 4.4        12         0       2      no
6           yes 11.02825 4.2        12         0       1      no
> 
> glm1 <- glm(participation~income, data=SwissLabor, family=binomial)
> -2 * logLik(glm1) # logLik は対数尤度を返す関数
[1] 1175.972
attr(,"nobs")
[1] 872
attr(,"df")
[1] 2
attr(,"class")
[1] "logLik"
> deviance(glm1) # deviance は逸脱度を返す関数
[1] 1175.972
ところが、OLSで線形回帰分析を行うと (1)式の通りにはならない。
> data(CPSSW04)
> head(CPSSW04)
  earnings     degree gender age
1 34.61538   bachelor   male  30
2 19.23077   bachelor female  30
3 13.73626 highschool female  30
4 19.23077   bachelor female  30
5 19.23077   bachelor   male  25
6 38.46154   bachelor female  32
> lm1 <- lm(earnings ~ degree, data=CPSSW04)
> -2 * logLik(lm1)
[1] 56150.05
attr(,"nall")
[1] 7986
attr(,"nobs")
[1] 7986
attr(,"df")
[1] 3
attr(,"class")
[1] "logLik"
> deviance(lm1)
[1] 528939.3
このような違いが生じるのは、deviance は OLS の結果に対しては、"−2 × 対数尤度" ではなく、残差平方和 (Residual Sum of Square: RSS) を返すからである。RSS を計算すると、
> sum(residuals(lm1)^2) # residuals はすべてのケースの残差のベクトルを返す
[1] 528939.3
となり、一致していることがわかる。

なぜ deviance 関数がこのような不規則な挙動をするのか、よくは知らないのだが、それほど不思議なことではない。実際、逸脱度はしばしば RSS に対応するものとしてテキストで記述されている(私も似たようなことを書いた気がする)。なぜなら、どちらもモデルのフィッティングが悪いほど大きくなるからであり、これらを最小化するようにパラメータが推定されるからである。また逸脱度の代わりに N * log(RSS / N) が使われることもある(extractAIC ではそうである)。さらに逸脱度と RSS は単調増加の関係にある(RSSを最小化するのが OLS だが、これは残差が正規分布していれば最尤推定値を導く)ので、どちらも似たようなものだと考えることはできるのである。社会学者にとって、こういう話はどうでもいいのだが、教えるときに数値が一致しないと授業としてはすっきりしないので、少し調べてしまいました。

Rスクリプト覚え書き: 暦日の連番を作る
日時の形式のデータを作りたいことがまれにあるが、そのための関数。
# 暦日の連番を作る ?seq.POSIXt でヘルプが見られる
# ISOdate は数値から日時形式 (POSIXct) のオブジェクトを作る
x <- seq(ISOdate(2008,1,1), ISOdate(2013,3,1), "months")
#, 2008/1/1, 2008/2/1, ... 2013/3/1 というベクトルができる
Rスクリプト覚え書き: lattice での文字の大きさの調整
lattice の使い方を少し勉強したのでメモ。文字の大きさを指定する方法。下記のスクリプトで作った図は、このページにあります。
country <- c("西独", "東独", "英", "米", "ハンガリ", "スロベニア", "日") # 国名のラベルの用意
a <- c(6.98, 7.06, 7.21, 7.87, 7.07, 7.41, 7.28) # 各国の感情的コミットメントの平均値
c <- c(5.90, 5.62, 5.54, 5.48, 5.52, 5.24, 6.55) # 各国の継続コミットメントの平均値
names(a) <- names(c) <- country # ラベルをつける

library(lattice)
barchart(~c(a,c) | gl(2,7, labels=c("感情的", "継続")), 
         xlab=list(label="コミットメントの強さの平均", cex=2), # xlab をリスト形式で指定し、cex= でサイズ指定
         scales=list(cex=1.5),           # scales は軸の目盛りラベルの書式設定で、その中で cex= でサイズ指定
         par.strip.text=list(cex=1.5), ) # par.strip.text はストライプ中の文字書式。同じく cex=でサイズ指定
Rスクリプト覚書き:混合効果(マルチレベル)・順序ロジット
先日の研究会で順序ロジットモデルで、交差ランダム効果をつけたものをRで推定できるのか、という問題が話題になっていたので、調べてみた。現在では ordinal というパッケージがあって、このパッケージの clmm 関数を使うとできる。また clmm2 という関数もあり、こちらのほうがより柔軟にモデリングできるらしいが、実装 (implementation) が古いらしい。それゆえ、単純な順序ロジットで十分という人は clmm を使い、clmm では思ったようなモデルが組めないという人は clmm2 を使ったらいいように感じた。以下は clmm のスクリプトの例。

> library(ordinal)
> 
> ## サンプル・データの取得と中身の確認
> data(soup)
> head(soup)
  RESP PROD PRODID SURENESS DAY SOUPTYPE SOUPFREQ COLD EASY GENDER AGEGROUP
1    1  Ref      1        6   1   Canned  >1/week  Yes    7 Female    51-65
2    1 Test      2        5   1   Canned  >1/week  Yes    7 Female    51-65
3    1  Ref      1        5   1   Canned  >1/week  Yes    7 Female    51-65
4    1 Test      3        6   1   Canned  >1/week  Yes    7 Female    51-65
5    1  Ref      1        5   2   Canned  >1/week  Yes    7 Female    51-65
6    1 Test      6        5   2   Canned  >1/week  Yes    7 Female    51-65
  LOCATION
1 Region 1
2 Region 1
3 Region 1
4 Region 1
5 Region 1
6 Region 1
> summary(soup)
      RESP        PROD      PRODID  SURENESS DAY          SOUPTYPE  
 1      :  10   Ref : 739   1:739   1:228    1:921   Self-made:958  
 2      :  10   Test:1108   2:369   2:260    2:926   Canned   :589  
 3      :  10               3:184   3:115            Dry-mix  :300  
 4      :  10               4:185   4: 98                           
 5      :  10               5:185   5:277                           
 6      :  10               6:185   6:869                           
 (Other):1787                                                       
      SOUPFREQ    COLD           EASY        GENDER      AGEGROUP  
 >1/week  :799   No :1608   7      :470   Male  : 609   18-30:349  
 1-4/month:948   Yes: 239   6      :388   Female:1238   31-40:460  
 <1/month :100              8      :270                 41-50:400  
                            5      :210                 51-65:638  
                            4      :190                            
                            9      :110                            
                            (Other):209                            
     LOCATION  
 Region 1:539  
 Region 2:588  
 Region 3:720  
               
               
               
               
> 
> ## 交差ランダム効果・順序ロジットモデルの推定
> mm1 <- clmm(SURENESS ~ PROD + (1|RESP) + (1|PRODID), data = soup)
> summary(mm1)
Cumulative Link Mixed Model fitted with the Laplace approximation

formula: SURENESS ~ PROD + (1 | RESP) + (1 | PRODID)
data:    soup

 link  threshold nobs logLik   AIC     niter    max.grad cond.H 
 logit flexible  1847 -2667.58 5351.16 24(1261) 9.67e-07 1.7e+03

Random effects:
           Var Std.Dev
RESP   0.32628  0.5712
PRODID 0.06658  0.2580
Number of groups:  RESP 185,  PRODID 6 

Coefficients:
         Estimate Std. Error z value Pr(>|z|)    
PRODTest   1.2708     0.2983    4.26 2.04e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Threshold coefficients:
    Estimate Std. Error z value
1|2  -1.4867     0.2750  -5.406
2|3  -0.4601     0.2714  -1.696
3|4  -0.1202     0.2711  -0.443
4|5   0.1473     0.2711   0.543
5|6   0.8607     0.2720   3.165
> 
> 
> ## lme4 や ordinal ではランダム効果の分散の検定ができてないが、
> ## ランダム効果無しのモデルとフィッティングを比較する尤度比検定を行えばよい。
> 
> mm2 <- clmm(SURENESS ~ PROD + (1|RESP), data = soup)
> mm3 <- clmm(SURENESS ~ PROD + (1|PRODID), data = soup)
> 
> 
> anova(mm2, mm1)
Likelihood ratio tests of cumulative link models:
 
    formula:                                    link: threshold:
mm2 SURENESS ~ PROD + (1 | RESP)                logit flexible  
mm1 SURENESS ~ PROD + (1 | RESP) + (1 | PRODID) logit flexible  

    no.par    AIC  logLik LR.stat df Pr(>Chisq)    
mm2      7 5361.4 -2673.7                          
mm1      8 5351.2 -2667.6  12.259  1  0.0004631 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> 
> anova(mm3, mm1)
Likelihood ratio tests of cumulative link models:
 
    formula:                                    link: threshold:
mm3 SURENESS ~ PROD + (1 | PRODID)              logit flexible  
mm1 SURENESS ~ PROD + (1 | RESP) + (1 | PRODID) logit flexible  

    no.par    AIC  logLik LR.stat df Pr(>Chisq)    
mm3      7 5384.4 -2685.2                          
mm1      8 5351.2 -2667.6   35.28  1  2.855e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
> 
> 
> ## ランダム効果を除いたふつうの順序ロジット・モデルの推定 clm 関数を使う
> mm0 <- clm(SURENESS ~ PROD, data = soup)
> 
> anova(mm0, mm1)
Likelihood ratio tests of cumulative link models:
 
    formula:                                    link: threshold:
mm0 SURENESS ~ PROD                             logit flexible  
mm1 SURENESS ~ PROD + (1 | RESP) + (1 | PRODID) logit flexible  

    no.par    AIC  logLik LR.stat df Pr(>Chisq)    
mm0      6 5392.7 -2690.3                          
mm1      8 5351.2 -2667.6  45.506  2  1.314e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
> 
Rスクリプト覚書き:vglm関数で平行性の仮定を置かずに順序ロジット
先日の数理社会学会で平行性の仮定を置かない順序ロジットについて阪大の院生の人たちが報告していたので、R でできないのか気になっていたのだが、vglm 関数を使えばそんなに難しくないことが分かったので、スクリプトをメモしておきます。
library(VGAM)
head(ugss) # 学部生のライフスタイル調査 
?ugss

# 目的変数の準備:ピアスの穴の数(なし、1〜2個、3個以上の3カテゴリの順序変数を作る)
pierce <- rep("None", nrow(ugss))
pierce[ugss$piercings==1 | ugss$piercings==2] <- "One or Two"
pierce[ugss$piercings > 2] <- "Three or More"
pierce <- factor(pierce)
xtabs(~pierce)

# 平行性を仮定しない
vglm1 <- vglm(pierce ~ sex + status, data=ugss, family=cumulative(reverse=T))
summary(vglm1)

# 平行性を仮定する
vglm2 <- vglm(pierce ~ sex + status, data=ugss, family=cumulative(reverse=T, parallel=T))
summary(vglm2)

# sex については平行性を仮定するが、status については平行性を仮定しない
vglm3 <- vglm(pierce ~ sex + status, data=ugss, family=cumulative(reverse=T, parallel=F~status))
summary(vglm3)

#モデルの比較(推定パラメータ数の昇順でモデルの順番を並べ替えてある)
dev1 <- c(deviance(vglm2), deviance(vglm3), deviance(vglm1)) # 逸脱度
df1 <- c(vglm2@df.residual, vglm3@df.residual, vglm1@df.residual) # 残差の自由度
np1 <- 1608 - df1 # 推定パラメータ数
rbind(dev1, df1, np1)

# AICの計算
dev1 + 2 * np1

# 上記3モデルの比較:尤度比検定の有意確率の計算
1 - pchisq(dev1[1] - dev1[2], df1[1] - df1[2]) 
1 - pchisq(dev1[1] - dev1[3], df1[1] - df1[3])
1 - pchisq(dev1[2] - dev1[3], df1[2] - df1[3])

Rスクリプト覚書: 多次元正規分布する乱数の生成 rmvnorm
多次元正規分布する架空のデータを作りたいなーと思うことがたまにあるので覚書。以下は3変数からなる10ケースの架空のサンプルを作るスクリプト。
# まず変数間の共分散行列を作る。
(covariance1 <- matrix(c(2, 1, 2,
                        1, 3, 0,
                        2, 0, 4), 3,3))

# サンプルサイズ、各変数の平均値、そして共分散行列を指定する
data1 <- rmvnorm(10, mean=1:3, sigma=covariance1)

# どんな分布になったか確認
plot(as.data.frame(data1))
Rスクリプト覚書: リストをベクトルに変換する関数 unlist
リストをベクトルに変換したいときがある。例えば3学年分の試験の点数のデータがあり、学年ごとに偏差値を計算したい。偏差値は scale(成績) * 10 + 50 で得られる(scale はz得点を返す関数)。グループ別に計算するには、by を使えばよい。ところが by の返り値はリストであるため、そのままでは分析に使えない。そこで unlist を使ってベクトルに変換する。例えば以下のように。
> # 架空のサンプルデータ
> d1
   学年 成績
1     1   61
2     1   70
3     1   68
4     1   81
5     1   46
6     1   55
7     2   54
8     2   63
9     2   70
10    2   62
11    3   77
12    3   53
> 
> # 偏差値を計算する関数を作る
> hensa <- function(x){ scale(x) * 10 + 50}
> 
> # byの返り値はリスト形式なのでunlistでベクトルに変換する
> d1$偏差値 <- unlist( by(d1$成績, d1$学年, hensa) )
> d1
   学年 成績   偏差値
1     1   61 47.96350
2     1   70 55.29489
3     1   68 53.66569
4     1   81 64.25547
5     1   46 35.74453
6     1   55 43.07592
7     2   54 37.40666
8     2   63 51.14485
9     2   70 61.83011
10    2   62 49.61838
11    3   77 57.07107
12    3   53 42.92893

Copyright (C) 2004 paperboy&co. All Rights Reserved.

Powered by "JUGEM"