myscratchpad

今までイマイチよくわからなかったデータの並べ替え

dat <- data.frame(a=factor(c("a1", "a1", "a1", "a1", "a1", "a2", "a2", "a2", "a2", "a2")), s=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), b1=c(3, 3, 1, 3, 5, 3, 5, 2, 4, 6), b2=c(4, 3, 4, 5, 7, 2, 6, 3, 6, 4), b3=c(6, 6, 6, 4, 8, 3, 2, 3, 6, 5), b4=c(5, 7, 8, 7, 9, 2, 3, 3, 4, 6))
# データフレームの中身

a	s	b1	b2	b3	b4
a1	1	3	4	6	5
a1	2	3	3	6	7
a1	3	1	4	6	8
a1	4	3	5	4	7
a1	5	5	7	8	9
a2	6	3	2	3	2
a2	7	5	6	2	3
a2	8	2	3	3	3
a2	9	4	6	6	4
a2	10	6	4	5	6

これをb1で並べ替える

class(dat$b1) # 型の確認。文字列では並べ替えできない
odr <- order(dat$b1) # odrというオブジェクトに、並べ替えの"行番号"が保存される
odr
 [1]  3  8  1  2  4  6  9  5  7 10
dat2 <- dat[odr,] # odrに保存された行のデータを取り出す。3行目のデータ、8行目のデータ…という感じ

## 型を確認しておく。数値がfactorだったら数値になおす
x <- factor(1:10)
x

class(x)
y <- as.numeric(as.character(x)) # as.numeric(x)だとfactorのlevelが適用されてしまうので
class(y)

各水準間の差得点を条件間で比較するだけ
入戸野先生のホームページの情報、pdfが大変分かりやすかった。感謝

# データの準備。データ解析テクニカルブックのp. 116より
dat <- data.frame(
s=factor(c("s1", "s2", "s3", "s4", "s5")),
a1b1=c(3, 3, 1, 3, 5), a1b2=c(4, 3, 4, 5, 7), a1b3=c(6, 6, 6, 4, 8), a1b4=c(5, 7, 8, 7, 9),
a2b1=c(3, 5, 2, 4, 6), a2b2=c(2, 6, 3, 6, 4), a2b3=c(3, 2, 3, 6, 5), a2b4=c(2, 3, 3, 4, 6)
)
dat
sapply(dat[2:9], sum)
sapply(dat[2:9], mean)

options(contrasts = c("contr.sum", "contr.sum"))
lmres <- lm(cbind(a1b1,a1b2,a1b3,a1b4,a2b1,a2b2,a2b3,a2b4)~1, dat)
library(car)
afact <- factor(c(rep("a1", 4), rep("a2", 4)))
bfact <- factor(c(rep(c("b1","b2","b3","b4"),2)))
idat <- data.frame(afact, bfact)
Anovares <- Anova(lmres, idata=idat, idesign=~afact*bfact, type=3)
summary(Anovares)

# 交互作用対比の準備
data1 <- dat[2:5]
data2 <- dat[6:9]
commat <- combn(x=1:4, 2) # combn関数で組み合わせをつくる
combnN <- ncol(commat)
c1 <- commat[1,]
c2 <- commat[2,]

# t検定結果格納用のオブジェクトを用意する
tres <- data.frame(matrix(1, 6, 9))
colnames(tres) <- c("d1mean", "d2mean", "d1sd", "d2sd", "t.value", "df", "p.value", "CIlower", "CIupper")
tresam <- data.frame(matrix(1,6,2)) # 文字列用
colnames(tresam) <- c("alternative", "method")
rwnames <- letters[1:6] # 行名用
# t検定繰り返し
for (i in 1:combnN) {
d11 <- data1[c1[i]]
d12 <- data1[c2[i]]
d21 <- data2[c1[i]]
d22 <- data2[c2[i]]
d1 <- d11[,1]-d12[,1]
d2 <- d21[,1]-d22[,1]
tobject <- t.test(d1, d2, alternative="two.sided", paired=T)
tres.vct <- c(mean(d1), mean(d2), sd(d1), sd(d2), tobject[[1]], tobject[[2]], tobject[[3]], tobject[[4]][1], tobject[[4]][2])
tres[i,] <- tres.vct
tresam[i,] <- c(tobject[[7]], tobject[[8]])
cntname <- paste(names(d11), "-", names(d12), " vs. ", names(d21), "-", names(d22), sep="")
rwnames[i] <- cntname
}

rownames(tres) <- rwnames
tres <- data.frame(tres, tresam)
# t検定結果
tres
# holm法でp値の調整
(p.holm <- p.adjust(tres["p.value"][,1], "holm"))
(tresp <- data.frame(tres, p.holm=p.holm))

## 入戸野先生の解説は有意水準を=.05として、これをp値の順位で割って有意水準自体を調整している。
hprank <- c(3,5,6,1,4,2) # 差の小ささ。最も差が小さいのは3番目ということ
hp <- rep(0.05, 6)
hp/hprank # これとtres["p.value"][,1]を比較する
(tresp2 <- data.frame(tres, p.holm=p.holm, p.ntn=hp/hprank))
## なので、有意水準を.05とし、tres["p.value"][,1]を順位倍して.05を超えるかどうか見ても同じ。p.adjustでやってるのはこれ。ただし、第n+1ステップが第nステップの有意水準より小さくならないように調整する。小さくなった場合はn+1の水準を使う


## 別のやり方
dat2 <- dat[c(2,5,6,9)]
options()$contrasts

lmres2 <- lm(cbind(a1b1, a1b4, a2b1, a2b4)~1, dat2)
afact2 <- factor(c("a1", "a1", "a2", "a2"))
bfact2 <- factor(c("b1", "b4", "b1", "b4"))
idat2 <- data.frame(afact2, bfact2)
idat2
library(car)
Anovares2 <- Anova(lmres2, idata=idat2, idesign=~afact2*bfact2, type=3)
Ares2 <- Astat(Anovares2)
Ares2$unv
tres[3,]

## グラフにしてみる
pmat <- matrix(sapply(dat2, mean), 2)
dimnames(pmat) <- list(c("b1", "b4"), c("a1", "a2"))

palette(grey.colors(2))
barplot(pmat, beside=T)
legend("top", rownames(pmat), fill = 1:2)
win.graph()
barplot(t(pmat), beside=T)
legend("top", colnames(pmat), fill = 1:2)

固有値とスクリープロット
library(psych)
faev <- fa.parallel(dat3)
faev

因子推定
# 最尤法
factMLres <- fa(r=dat3, nfactors=3 ,rotate="promax", fm="ml", scores=T)
# 主因子法
factPAres <- fa(r=dat3, nfactors=3 ,rotate="promax", fm="pa", scores=T)
# 一般化最小二乗法
factGLSres <- fa(r=dat3, nfactors=3 ,rotate="promax", fm="gls", scores=T)
# 重みつき最小二乗法
factWLSres <- fa(r=dat3, nfactors=3 ,rotate="promax", fm="wls", scores=T)
# 最小残差法 (重みづけない最小二乗法)
factOLSres <- fa(r=dat3, nfactors=3 ,rotate="promax", fm="minres", scores=T)

# 結果の出力
print(factMLres, cut=0, sort=T, digits=5)

回転 (rotate="") のオプション
# 回転なし。いわゆる初期解ってやつ
"none"
# 直交解
"varimax", "quartimax", "bentlerT", and "geominT"
# 斜交解
"promax", "oblimin", "simplimax", "bentlerQ, and "geominQ" or "cluster"

因子得点はscores=T。推定は (たぶん) 回帰法

library(psych)
? fa
のDetailも勉強になる

いつの間にか、psychパッケージの因子分析のオプションが激増していた。William Revelle先生に感謝。


# データの準備
dat <- read.table("https://blog.cnobi.jp/v1/blog/user/89d80905c7038b4121822249e9062fba/1261838567", header = T)
# UCLAのページから拝借。感謝。
## 0をNAで置換する。この記事のコメントで教えていただいた。感謝
dat2 <- dat
dat2[dat2==0] <- NA
dat3 <- na.omit(dat2) # NAを除外する


# 記述統計
## 変数名を保存しておく
orname <- data.frame(itnumber=names(dat), itname=c("INSTRUC WELL PREPARED","INSTRUC SCHOLARLY GRASP","INSTRUCTOR CONFIDENCE","INSTRUCTOR FOCUS LECTURES","INSTRUCTOR USES CLEAR RELEVANT EXAMPLES","INSTRUCTOR SENSITIVE TO STUDENTS","INSTRUCTOR ALLOWS ME TO ASK QUESTIONS","INSTRUCTOR IS ACCESSIBLE TO STUDENTS OUTSIDE CLASS","INSTRUCTOR AWARE OF STUDENTS UNDERSTANDING","I AM SATISFIED WITH STUDENT PERFORMANCE EVALUATION","COMPARED TO OTHER INSTRUCTORS, THIS INSTRUCTOR IS","COMPARED TO OTHER COURSES THIS COURSE WAS"))
dscrb <- data.frame(ITEM=orname$itname, MEAN=mean(dat3), SD=sd(dat3), N=sapply(dat3, length))
dscrb

# スクリープロットと平行分析
library(psych)
fa.parallel(dat3)

# psychパッケージ、fa関数で最尤解の因子分析。プロマックス回転
library(psych)
factMLres <- fa(dat3, 3 ,rotate="promax", fm="ml", scores=T) # 因子得点を出すためにはscores=Tとしておく
print(factMLres, cut=0, sort=T, digits=5) # factanal関数と違って、fa関数はきちんと並べ替えされる

## factanal関数と比較してみる
fctanlres <- factanal(data=dat3, factors=3, rotation="promax", scores="regression")
print(fctanlres, cut=0.00001, digits=5)

## いずれにせよ自由度33のカイ二乗値=135.35, p < .001なので、この因子モデルは適合しない。因子数を増やして再度分析する
factMLres4 <- fa(dat3, 4 ,rotate="promax", fm="ml", scores=T)
factMLres5 <- fa(dat3, 5 ,rotate="promax", fm="ml", scores=T)
factMLres6 <- fa(dat3, 6 ,rotate="promax", fm="ml", scores=T)
## カイ二乗値とp値だけとりだす
factMLres$STATISTIC
factMLres$PVAL
factMLres4$STATISTIC
factMLres4$PVAL
factMLres5$STATISTIC
factMLres5$PVAL
factMLres6$STATISTIC
factMLres6$PVAL
## 5因子解が適合した (ただし、カイ二乗値による結果が絶対ではない。おそらく3因子解が適切だろう)
print(factMLres5, cut=0, sort=T, digits=5)
## 因子得点
factMLres5$scores



# psychパッケージのfa関数はfactanal関数で出力されなかったものも出てくるので便利。負荷量は一致するが、どうも平方和の計算がなんか違う。因子の順番を変えて計算しているようだ？

# データフレームの変数ごとに欠損値があるかどうか
sapply(dat, function(x) sum(is.na(x)))

# 特定の値を欠損値とする
dat[dat=="setosa"] <- NA
dat[dat==9999] <- NA # 数値もできる

# 変数を指定して欠損値とする。データ内に既に欠損値があるとエラーになる
dat[dat$Species=="virginica", "Species"] <- NA

# 欠損値の除去
# 関数内では大体na.rm=Tとする
mean(x, na.rm=T)

# データフレームから行 (ケース) ごと除去
na.omit(dat)


Amelia II, Mice, mitools といったパッケージがあるらしい

# おまけ
dat2 <- dat
for(i in 1:ncol(dat)){
dat2[,i] <- replace(dat[,i], which(dat[,i]==0), NA)
}