myscratchpad

Rで確認的因子分析を行う。
Indiana Universityのサイトとデータを参考にする。感謝。

library(foreign)
dat <- data.frame(read.spss("http://www.indiana.edu/~statmath/stat/all/cfa/values.sav"))

library(sem)
model <- specify.model()
economic   ->  privtown,  p1,   1
economic   ->  govtresp,  p2,   NA
economic   ->  compete,   p3,   NA
moral      ->  homosex,   p4,   NA
moral      ->  abortion,  p5,   NA
moral      ->  euthanas,  p6,   NA
privtown  <->  privtown,  e1,   NA
govtresp  <->  govtresp,  e2,   NA
compete   <->  compete,   e3,   NA
homosex   <->  homosex,   e4,   NA
abortion  <->  abortion,  e5,   NA
euthanas  <->  euthanas,  e6,   NA
economic  <->  economic,  NA,   1
moral     <->  moral,     NA,   1
economic  <->  moral,     em,   NA

semres <- sem(model, cov(dat), nrow(dat))
(sumsem <- summary(semres))
std.coef(semres) # 標準化係数
# このサイトとこのサイトにある関数を利用してAIC, BIC, R二乗を出す。感謝
source("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1266106324")
sem.aic(semres) # AIC
rsquare.sem(semres) #  R二乗
sumsem$chisq - sumsem$df * log(semres$N)  # BIC


## 作図。Rgraphvizが上手くインストールできないので (libcdt-4.dllがないとかいわれる) のでGraphvizで書く
 # とりあえずインストールはできた
# 以下のコードで出力されたdigraph...のテキストをGraphvizで描画する。Graphvizを起動し、File -> Newで出たウィンドウにコピペしてF5を押す。設定ダイアログは何もせずにOKとした
path.diagram(semres, ignore.double=FALSE, edge.labels="values", digits=3, standardize=T) # file="filename"と指定するとpdfが作業ディレクトリに出力される。だけど表示がおかしい

出力されたテキスト
# 望んだ作図ではないが、今日はこの辺で勘弁してやる


元サイトと全然結果が一致しない…まいったのう

## specify.model の覚書
変数間の関係、共分散・係数の名前、初期値

変数間の関係は -> か <-> で指定。前後に変数名を入れる
 sem関数に与える相関・共分散行列にない名前は潜在変数になる
その次は共分散の名前。p1とか文字列で指定する。
その次が初期値。NAこれを推定する。
一方向 -> は回帰係数
両方向 <-> で変数名が違う場合は共分散
両方向で変数名が同じ場合はその変数の分散

係数をNAとすると固定母数として扱われる
潜在変数 (因子) の分散は固定母数1にすること
F1 <-> F1, NA, 1
観測変数の分散はNAにする
O1 <-> O1, e1, NA

同じ名前のパスが2つあると同値制約
同じ潜在変数同士　F1 <-F1 は負荷量にNA, 初期値に1を与えて固定母数にする
初心者の方へ: モデル式中にはパス係数、分散、共分散の全てを表記しなければならない。パス係数だけ指定してもモデルが推定されないので、要注意である (Rjiwikiより)

psychパッケージのstructure.sem関数でもモデルが書ける

下のspssのページをもとにした。
http://www.nyu.edu/its/statistics/Docs/intracls.html
下のページも参照。感謝
http://www.hs.hirosaki-u.ac.jp/~pteiki/research/stat/S/icc/

dat <- data.frame(
rater1=c(9,6,8,7,10,6),
rater2=c(2,1,4,1,5,2),
rater3=c(5,3,6,2,6,4),
rater4=c(8,2,8,6,9,7)
)

# psycパッケージのICC関数を使う。psychパッケージはマジ便利
library(psych)
ICC(dat)
## 勉強のため、ICCのコードの一部を貼り付けておく
x <- dat
n.obs <- dim(x)[1]
nj <- dim(x)[2]
x.s <- stack(x)
x.df <- data.frame(x.s, subs = rep(paste("S", 1:n.obs, sep = ""),
nj))
aov.x <- aov(values ~ subs + ind, data = x.df)
s.aov <- summary(aov.x)
stats <- matrix(unlist(s.aov), ncol = 3, byrow = TRUE)
MSB <- stats[3, 1]
MSW <- (stats[2, 2] + stats[2, 3])/(stats[1, 2] + stats[1,
3])
MSJ <- stats[3, 2]
MSE <- stats[3, 3]
(ICC1 <- (MSB - MSW)/(MSB + (nj - 1) * MSW)) # ICC1. Single raters absolute
(ICC2 <- (MSB - MSE)/(MSB + (nj - 1) * MSE + nj * (MSJ - MSE)/n.obs)) # ICC2. Single random raters
(ICC3 <- (MSB - MSE)/(MSB + (nj - 1) * MSE)) # ICC3. Single fixed raters
(ICC12 <- (MSB - MSW)/(MSB)) # ICC1k. Average random raters
(ICC22 <- (MSB - MSE)/(MSB + (MSJ - MSE)/n.obs)) # ICC2k. Average random raters
(ICC32 <- (MSB - MSE)/MSB) # ICC3k. Average fixed raters


head(x.df)
aov.x
s.aov
round(stats, 3)
alpha(dat) # 質問紙の信頼性係数を求めるときと同じだが、列には項目ではなく評定者がならぶ

因子変数と数値変数が混在しているデータフレームで、因子だけのと数値だけのデータフレームに分ける
# データの準備
f1 <- sample(letters[1:5], 20, replace=T)
f2 <- ordered(sample(letters[11:15], 20, replace=T))
f2["Levels"] # これは順序変数にしてみよう
n1 <- sample(500:1000, 20, replace=T)
n2 <- sample(rnorm(100), 20, replace=T)
f3 <- sample(letters[21:25], 20, replace=T)
dat <- data.frame(f1, f2, n1, n2, f3)
summary(dat)
sapply(dat, class)

# 因子変数はdatfへ、それ以外 (数値) はdatiへ
fcol <- grep("factor", sapply(dat, class))
datf <- dat[,c(fcol)]
dati <- dat[,-c(fcol)]
datf
dati

library(psych)
data(bfi) # psychパッケージ内のデータセット
dat <- bfi

# 各変数ごとに反応の種類を調べる
sapply(dat, function(x) levels(factor(x)))
  # 変数ごとに反応や単位が違う場合に使う。データはデータフレームにしておくこと
  # library(psych); data(sat.act); sapply(sat.act, function(x) levels(factor(x)))
# 全体の回答数
summary(factor(as.matrix(dat)))
# 一応プロットもしてみる
win.graph()
plot(as.vector(as.matrix(dat)))
win.graph()
hist(as.vector(as.matrix(dat)))

# 各項目ごとの反応数を集計する。
# 有効回答 (NA, Inf, NaN以外の反応) 種類をとりだす
(res <- names(table(as.matrix(dat)))) # table関数はNAとかは集計しないことを利用
resN <- length(res) # 反応の種類数
itemN <- ncol(dat) # 項目数
count.dat <- data.frame(matrix(NA, itemN, resN)) # 集計表用の行列。行が項目、列が反応
  rownames(count.dat) <- names(dat)
  colnames(count.dat) <- res
for (i in 1:resN){
count.dat[,i] <- sapply(dat, function(x) sum(na.omit(x)==res[i]))
}
count.dat

# NA, NaN, Infも調べる
na <- sapply(dat, function(x) sum(is.na(x))) # NAとNaNがカウントされる
inf <- sapply(dat, function(x) sum(is.infinite(x))) # Infをカウント

# 記述統計とかとくっつける。記述統計はpsychパッケージのdescribe関数で
descriptives <- cbind(data.frame(describe(dat)), count.dat, na, inf)
descriptives
# csvで保存
write.csv(descriptives, file="descriptives.csv", quote=FALSE)

データフレームの列ごとに数値をカウントする。
Excelでいうとcountifみたいなものか

dat <- data.frame(matrix(sample(1:5, 2000, replace=T), 200, 10)) # 10項目5件法の尺度に200人が答えた、という仮想データ
count.dat <- matrix(NA, 5, 10) # 集計表用の行列
for (i in 1:5){
count.dat[i,] <- sapply(dat, function(x) sum(x==i))
}
colnames(count.dat) <- colnames(dat)
rownames(count.dat) <- 1:5
count.dat # 集計表の表示
rowSums(count.dat) # 回答ごとの集計
colSums(count.dat) # 項目ごとの全反応数