myscratchpad

各種統計ソフトの解析例。UCLA
http://www.ats.ucla.edu/stat/

spssでは標準で出力される (らしい) Kaiser-Meyer-Olkinの標本妥当性と、Bartlettの球面性検定の関数を見つけたのでメモしておく

# Kaiser-Meyer-Olkinの標本妥当性
青木先生のページから。感謝。
kmo(dat) # datは数値のみのデータフレームか行列

全く同じものがこのページにもあった
kmo.test <- function(df)
{
        cor.sq = cor(df)^2
    cor.sumsq = (sum(cor.sq)-dim(cor.sq)[1])/2
    library(corpcor)
    pcor.sq = cor2pcor(cor(df))^2
    pcor.sumsq = (sum(pcor.sq)-dim(pcor.sq)[1])/2
    kmo = cor.sumsq/(cor.sumsq+pcor.sumsq)
    return(kmo)
}


# Bartlettの球面性
# psychパッケージは偉大すぎる
library(psych)
cortest.bartlettdat)
#cortest.bartlett(cor(dat)) # 相関行列を与えること、とマニュアルにはあるが、spssの結果と一致するのはローデータを与えたときのようだ

主因子法 (principal factor) と主軸法 (principal axis) って同じなのかな？

dat <- read.table("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1261662305", header = T)
# 小塩先生のページを参考にさせていただいた。感謝
# http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/09_folder/da09_01.html

# 変数名が長いので別にしておこう
ornames <- names(dat)
dat2 <- dat
names(dat2) <- c(paste("A_0", 1:9, sep = ""), paste("A_", 10:23, sep = ""))

# まずは記述統計
mres <- mean(dat2)
sdres <- sd(dat2)
resdes <- data.frame(name=names(dat2), mean= mres, sd = sdres, mm1sd=mres-sdres, mp1sd=mres+sdres)
resdes
## 以下でもいい
#summary(dat2)
#sd(dat2)
#library(psych)
#describe(dat2)

# 固有値を調べプロットする
evres <- eigen(cor(dat2))
evres$value # 固有値を表示
plot(evres$value, type="b") # プロット。いわゆるスクリープロット
# psychパッケージのfa.parallelが便利
fa.parallel(dat2)

# psychパッケージのfactor.pa関数により主因子解の因子分析
factpares <- factor.pa(dat2, 3 ,rotate="promax")
print(factpares, cut=0, sort=T, digits=5)

# 小塩先生のページどおり、一部の項目を除外して因子分析
dat3 <- dat2[c(-19, -8, -23)] # A19, A08, A23を除外
fa.parallel(dat3)
factpares2 <- factor.pa(dat3, 3, rotate="promax")
print(factpares2, cut=0, sort=T, digits=5)

微妙に小塩先生のspssの結果と一致しないがまあいいだろう
因子分析なんてそんなものだ、ということにしよう

## 以下と同じ分析のはずだが、なぜか第2主軸の符号が逆になってる。たぶん上のでいいんだろう
fares <- fa(dat3, 3, rotate="promax", fm="pa")
print(fares, cut=0, sort=T, digits=5)

データファイルで個人ごとの代表値などをまとめる
tapply, by等も試したが、たぶんaggregate関数が一番便利

dat <- read.table("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1261414199", header = T)
dat # データファイル。
summary(dat)

ag.dat <- aggregate(x=dat[5:6], by=list(dat[,1],dat[,2], dat[,3], dat[,4]), FUN=mean)
ag.dat # 個人別に条件ごとのrtの平均値をだす
## tapplyはINDEXの3番目からリスト要素になる。tapplyはデータに2列指定できない
tapply(X=dat[,5], INDEX=list(dat[,1], dat[,2], dat[,3],dat[,4]), FUN=mean)

# 分析しやすいようにreshapeパッケージのcast関数で横長にする
library(reshape)
pn <- substr(ag.dat[,2], 1, 3) # 横長データセットの変数名をつくるため、substrで一部をとりだす
hl <- substr(ag.dat[,3], 1, 1)
on <- substr(ag.dat[,4], 1, 1)
variable <- paste(pn, hl, on, sep = "_") # variableという変数名じゃないとcast関数ではエラー
value <- ag.dat$rt # valueという変数名じゃないとcast関数ではエラー
pid <- ag.dat$Group.1
ag.w <- cast(data.frame(pid, value, variable))
ag.w

普通はy軸の目盛が左に倒れている
plot(1:10)

lasというパラメータで右に90°回転させる
win.graph()
plot(1:10, las = 1)
## 軸と水平、という指定
## axisの中でも指定できる
plot(1:10, yaxt = "n")
axis(2, las = 1)


# 以下参考
win.graph()
plot(1:10, las = 0)
win.graph()
plot(1:10, las = 1)
win.graph()
plot(1:10, las = 2)
win.graph()
plot(1:10, las = 3)