myscratchpad

# 小塩先生のホームページからデータを拝借、出力を参照。感謝。
# http://psy.isc.chubu.ac.jp/~oshiolab/teaching_folder/datakaiseki_folder/add_folder/daad_01.html

ex <- c(3,6,6,6,5,4,6,5,6,6,5,5,6,5,5,6,4,6,5,4)
so <- c(4,6,5,7,7,5,6,5,6,5,4,5,6,5,6,6,4,6,3,6)
ac <- c(4,7,7,5,6,5,7,4,6,6,4,6,5,4,4,6,3,7,4,6)
il <- c(5,8,5,4,5,5,6,5,7,6,5,5,5,4,5,4,6,4,3,3)
re <- c(4,7,5,6,5,6,4,5,7,5,5,4,6,5,6,4,5,5,5,5)
fl <- c(4,7,6,5,5,6,4,6,6,5,5,5,5,3,6,5,6,5,4,4)
dat <- data.frame(ex, so, ac, il, re, fl)

library(psych)

fa.parallel(dat) # 固有値の推移。PC Actual Dataが主成分解 (だと思う) 。

fit <- principal(dat, nfactors=2, rotate=F, scores=T) # 固有値1を超えるのは2つなのでnfactorsで指定する。

print(fit, digits=3)
fit$scores
cor(fit$scores)

plot(fit) # 主成分得点からみる散布図

# 結論：psychパッケージは偉大である。

# 足立浩平 (2006). 多変量データ解析法――心理・教育・社会系のための入門――　ナカニシヤ出版の例より
 # 因子分析と主成分分析は似て非なる方法であることに注意
# 以下のサイトを参考。というかRコマンダーでの例をそのままコードにしただけ。感謝
http://www.nanzan-u.ac.jp/~kamiya/r/content/rcmdrprin.html
http://www.statmethods.net/advstats/factor.html

# 表　3.1　p24
ppt <- c(88,52,77,35,60,97,48,66,78)
wrt <- c(70,78,87,40,43,95,62,66,50)
int <- c(65,88,89,43,40,91,83,65,48)
dat <- data.frame(ppt, wrt, int)
mean(dat)
sapply(dat, var)

fit <- princomp(~ppt+wrt+int, cor=F, data=dat) # 相関行列をもとにした分析をするにはcor=Tとする
fit2 <- prcomp(dat, cor=F)
 # 主成分得点 表3.1　p22
 # 単にprincomp(dat) でも同じ
fit$scores
fit2$x

 # 各主成分の係数　表3.5　p32
loadings(fit)
 unclass(loadings(fit)) # 行列にしただけ
fit2$rotation
 # 固有ベクトルなので以下と同じ
eigen(cov(dat))$vectors

 # 各主成分得点どうしの共分散　表3.6　p32
cov(fit$scores)
 round(cov(fit$scores),1) # みやすく四捨五入
 round(cov(fit2$x),1)

 # 成分行列 (いわゆる負荷量)
t(fit$sd * t(fit$loadings ) )[, drop = FALSE]
t(fit2$sdev * t(fit2$rotation ) )[, drop = FALSE]
 # 結果が異なる。おそらくprcompは Unlike princomp, variances are computed with the usual divisor N - 1. であるため、標準偏差の値が異なる
 fit$sd
 fit2$sdev
 sqrt(fit2$sdev^2 *(8/9)) # 不偏分散を標本に戻すと一致

 # 説明率、累積説明率　図3.6　p27
summary(fit)
print(summary(fit2), digits=10) #
基本的には同じ結果だが、prcompのsummaryは小数点第5桁までで区切っているのがわかる

 # 固有値
fit$sd^2
fit2$sdev^2
 # 固有値なので以下と同じ
eigen(cov(dat))$value


## おまけ。psychパッケージでもやってみる

fit2 <- princomp(~ppt+wrt+int, cor=T, data=dat) # 相関行列を使うよう指定してprincomp関数で実行
library(psych)
fit3 <- principal(dat, nfactors=3, rotate=F, scores=T)# psychパッケージでは標準化した相関行列がデフォルトで用いられる。回転もデフォルトだとバリマックス回転が行われる

summary(fit2) # princomp関数
print(fit3, digits=5) # psychパッケージのprincipal関数

fit2$scores
fit3$scores # psychパッケージでは標準化した主成分得点が算出される

## 覚書
 prcomp関数は各個体の平均から偏差を求めるQモード法。対してprincomp, principal関数はRモード法。Qモード法は変数の数 > サンプルサイズ　でも実行できる。

##
# 足立浩平 (2006). 多変量データ解析法――心理・教育・社会系のための入門――　ナカニシヤ出版


# 表2.1　p16
mat <- matrix(c(3.2,3.4,3,3.2,4.2,4,4,3.7,3.6,3.6,3.7,3.4,3.2,3.2,2.7,3.5,3.2,3.2,4.6,4,4.8,4.6,4.3,3.6,3.2,3.7,3.2,3.8,3.7,3.4,3.5,3.5,4.5,3.8,3.9,4.1,3.7,3.5,3.7,3.5,3.6,3.8,2.8,2.5,2.2,2.6,3.1,3.4,3.5,3.4,2.9,3.5,3.9,3.1,2.9,2.8,2.6,2.2,2.1,2.5,3,3.2,3.8,4,3.5,4.2,4.7,2.7,2.2,2.3,2.6,2.6,2.2,2.6,3.2,3.1,3.7,4.1,3.6,4.2,4.7,2.4,2.5,2.4,2.2,3.2,3.3,3.6,2.8,2.4,3.2,4.3,4.7,3.5,4.9,2.3,3.3,3.9,1.4,2.1,3.4,2.9,3.3,1.5,2.1,3.4,4.2,3.5,3.5,1.8,3.3,2.5,1.7,3.6,4.1,4.2,4.1,1.6,2.6,3.5,4.1,3.7,4.2,2.3,3.4,4.7,3.3,4.1,3.4,3.2,4.5,3.7,3.7,4.2,3.9,3.5,3.7,3.3,2.8,3.9,3.8,4.7,1.3,1.5,2.3,3.9,3.6,4.4,3.7,2.5,2.8,2.9,1.8,2.3,1.8,4.2,4.3,4,4.9,3,4.5,4,5,3.5,4.1,3.3,4.3,4.3),
nr=14,
dimnames=list(c("僧侶", "銀行員", "漫画家", "デザイナー", "保母", "大学教授", "医師", "警察官", "新聞記者", "船乗り", "プロスポーツ選手", "作家", "俳優", "スチュワーデス"),
c("立派な", "役立つ", "よい", "大きい", "力がある", "強い", "速い", "騒がしい", "若い", "誠実な", "かたい", "忙しい"))
)

# 表　3.1　p22
ppt <- c(88,52,77,35,60,97,48,66,78)
wrt <- c(70,78,87,40,43,95,62,66,50)
int <- c(65,88,89,43,40,91,83,65,48)
dat <- data.frame(ppt, wrt, int)
mean(dat)

# 表　3.2　p28
ppt <- c(77,60,97,35,66,48,78,52,88)
wrt <- c(18,8,20,8,14,12,10,16,14)
int <- c(45,20,45,22,32,41,24,44,32)
dat <- data.frame(ppt, wrt, int)
mean(dat)

# 表8-1　p76
se <- c(9,2,5,4,6,4,6,6,7,4,5,6,7,4,3,5,4,7,5,5,6,4,3,7,4,5,3,4,6,5,4,3,5,5,2,4,7,5,6,5,3,5,6,5,5,2,5,5,7,6,6,6,7,7,4,7,7,6,6,4,4,6,4,2,5,4,4,5,6,4,6,6,6,5,6,4,7,8,8,4,4,5,6,8,4,5,4,3,4,5,3,3,6,4,5,4,3,4,5,4)
yo <- c(7,3,6,6,5,5,7,6,6,4,6,4,5,5,6,6,5,6,7,5,7,6,6,7,6,6,5,6,5,6,3,3,5,6,5,6,7,7,5,5,4,7,4,7,7,5,9,5,6,5,6,6,6,6,4,6,5,6,5,4,6,4,5,6,6,5,6,4,6,6,5,7,5,6,6,4,5,6,5,5,6,6,6,6,6,7,4,5,6,6,5,3,6,5,7,4,7,7,4,5)
sa <- c(9,5,7,6,7,5,6,7,8,6,6,6,8,6,5,7,8,8,7,7,7,6,5,8,7,5,4,6,9,8,7,6,6,7,5,7,9,5,7,4,6,7,7,5,9,4,6,7,8,6,7,8,8,7,6,8,8,6,7,7,6,5,6,5,6,4,6,7,7,5,7,8,8,6,8,5,7,9,8,3,5,7,7,6,6,7,5,6,6,7,6,6,8,5,6,6,6,5,5,7)
bu <- c(2,8,6,3,6,5,5,4,5,8,4,5,5,7,6,3,4,4,4,5,4,4,7,5,6,4,6,5,5,5,7,9,7,5,5,7,5,3,4,7,9,6,7,4,3,7,4,8,5,7,3,4,4,5,8,3,6,3,6,6,5,7,6,5,6,7,4,6,4,7,5,5,5,6,5,6,6,6,4,6,6,6,6,7,6,4,9,6,5,6,7,7,5,6,4,7,4,6,9,7)
ha <- c(9,1,8,8,6,6,8,8,6,4,6,5,5,5,4,9,8,6,9,4,8,5,4,7,5,8,5,5,8,6,4,3,5,4,5,5,6,8,7,5,4,6,5,9,8,3,9,4,7,4,6,8,8,7,2,7,7,6,5,4,6,4,6,5,7,5,7,6,7,5,5,7,6,6,6,5,5,6,6,5,5,5,4,6,6,8,2,4,6,5,5,3,5,6,6,3,7,7,4,3)
ya <- c(8,3,4,5,5,3,3,7,5,3,4,6,6,4,2,5,5,6,5,4,6,3,3,8,4,5,4,4,6,6,5,4,6,4,3,5,8,5,6,3,3,4,5,6,7,4,4,5,7,6,5,8,6,6,3,8,6,4,5,5,5,3,4,2,6,3,4,8,4,4,6,7,6,2,6,4,7,8,9,3,4,4,4,8,5,6,2,3,4,4,4,3,7,4,6,5,3,3,5,5)
ch <- c(3,7,6,7,6,5,6,6,3,6,5,4,7,7,6,7,5,4,4,6,6,7,8,5,5,4,6,7,5,6,5,7,6,8,6,6,3,5,2,8,6,6,5,5,3,6,4,3,2,5,6,3,4,5,7,2,2,3,3,5,9,5,7,9,5,7,6,4,5,6,4,4,2,8,3,5,3,1,2,6,6,6,3,2,6,5,8,7,9,6,9,7,3,8,3,6,7,7,6,4)
ni <- c(8,3,6,7,6,5,8,7,4,3,7,6,6,5,5,7,6,6,7,8,6,5,2,6,6,7,5,5,6,6,5,3,5,8,5,5,5,5,7,4,4,7,4,6,7,5,8,5,6,4,7,9,7,7,3,7,6,4,7,4,6,4,5,5,6,5,6,4,8,6,6,7,5,5,6,3,6,7,7,5,5,5,5,7,6,9,4,4,6,6,5,4,7,5,5,5,5,4,7,3)
dat <- data.frame(se, yo, sa, bu, ha, ya, ch, ni)


# 表15.1　p150
grp <- rep(c("a", "b"), c(15, 12))
sha  <- c(15,11,16,19,18,15,17,12,13,14,16,11,20,15,13,11,10,11,10,10,13,11,15,12,10,12,10)
kyo  <- c(14,13,14,21,26,28,19,15,22,26,20,15,21,20,13,15,13,14,10,14,19,10,20,22,11,10,14)
kin  <- c(15,17,17,18,21,18,12,18,16,18,18,20,17,19,17,18,16,24,13,22,23,20,20,23,18,19,21)
shi  <- c(14,17,26,15,15,12,10,12,10,6,18,15,20,12,16,17,9,16,12,18,24,28,16,13,10,27,19) 
dat <- data.frame(grp, sha, kyo, kin, shi)


# 表15.4　p153
grp   <- factor(c(1,1,1,1,2,2,2,3,3,3,3,4,4,4,4,4,5,5,5,5))
sak  <- c(9,6,6,4,7,1,4,5,5,3,6,8,4,10,5,5,4,5,1,6)
ens   <- c(8,8,6,8,10,9,9,4,5,5,5,6,6,5,6,3,6,4,4,4)
tmp   <- c(2,5,4,8,5,9,7,9,10,8,7,4,7,6,6,8,3,2,1,4)
mss  <- c(6,7,8,8,6,7,5,6,7,5,5,4,4,3,3,5,5,6,5,4) 
arg   <- c(10,6,5,7,4,7,6,0,7,3,5,5,1,10,2,3,7,5,3,6)
org   <- c(7,9,5,6,3,4,6,8,5,5,6,6,5,6,4,5,6,9,7,10)
dat <- data.frame(grp, sak, ens, tmp, mss, arg, org)

st <- "a,b,c,d,e"
st2 <- "f,g,h"
st3 <- "i,j,k,l"
lst <- list(st, st2, st3)
# これを3行6列の表にしたい。空白はNAでうめる
(sts <- lapply(lst, function(x) unlist(strsplit(x, ","))))
lens <- lapply(sts, length)

reslist <- list()
for (i in 1:length(sts)) {
lens <- length(sts[[i]])
ko <- 6
sta <- c(sts[[i]], rep(NA, 6-lens))
reslist[[i]] <- sta
}

res2 <- t(data.frame(reslist))
dimnames(res2) <- NULL
data.frame(res2)

# 行ごとにx内の要素に一致するものの数を数える
apply(res2, 1, function(a) sum(complete.cases(match(a,x))))
# 列ごとにx内の要素に一致するものの数を数える
apply(res2, 2, function(a) sum(complete.cases(match(a,x))))

# 棒グラフの中に数値を表示させる

x <- matrix(c(10,10,20, 15, 5,20, 25,10,5, 5,5,10), nr=3,
dimnames=list(paste("r", 1:3, sep=""), paste("c", 1:4, sep="")))
x
cprp(x) #　この関数はこれ

## 行と列の順番を変える
 # x <- t(x) # 行列の転置
 # x <- x[nrow(x):1, ncol(x):1] # 順序を変える

par(mfrow=c(2,3))
## 棒グラフ縦型に数値をいれる
bp <- barplot(x, beside=T, main="棒縦")
for (i in 1:ncol(x)) {
xl <- bp[,i]
yl <- x[,i]+2
 # yl <- x[,i]/2 # 中心にいれるとき
txt <- round(x[,i], 2)
txt[txt=="0"] <- NA
text(xl, yl, txt)
}
 # 凡例のコードは共通
 legend("topright", legend=rownames(x), fill=grey.colors(nrow(x)))

## 積み上げ棒グラフ縦型に数値を入れる
 # x <- t(x) # 行と列を逆にしたいなら。
bp <- barplot(x, main="積み縦")
for (i in 1:ncol(x)) {
xl <- bp[i]
yl <- cumsum(x[,i]) - x[,i]/2
txt <- round(x[,i], 2)
txt[txt=="0"] <- NA
text(xl, yl, txt)
}

## 帯グラフ縦型に数値を入れる
x
cx <- cprp(x)$cp
bp <- barplot(cx, main="帯縦\n 数値の色も変えてみた")
for (i in 1:ncol(x)) {
xl <- bp[i]
yl <- cumsum(cx[,i]) - cx[,i]/2
txt <- round(x[,i], 2)
txt[txt=="0"] <- NA
 # パーセントをいれるとき
 # txt <- round(cx[,i], 2); txt[txt=="0"] <- NA
text(xl, yl, txt, col=c("white", "black", "blue"))
}

## 棒グラフ横型に数値をいれる
bp <- barplot(x, beside=T, horiz=T, main="棒横")
for (i in 1:ncol(x)) {
xl <- bp[,i]
yl <- x[,i]+2
 # yl <- x[,i]/2 # 中心にいれるとき
txt <- round(x[,i], 2)
txt[txt=="0"] <- NA
text(yl, xl, txt)
}

## 積み上げ棒グラフ横に数値を入れる
bp <- barplot(x, horiz=T, main="積み横")
for (i in 1:ncol(x)) {
xl <- bp[i]
yl <- cumsum(x[,i]) - x[,i]/2
txt <- round(x[,i], 2)
txt[txt=="0"] <- NA
text(yl, xl, txt) # 縦型のときとxl, ylを逆にする
}

## 帯グラフ横型に数値を入れる
x
cx <- cprp(x)$cp

bp <- barplot(cx, horiz=T, main="帯横\n軸を%にする", xaxt="n")
axis(1, at=seq(0,1,0.2), labels=paste(seq(0,100,20), "%"))

for (i in 1:ncol(x)) {
xl <- bp[i]
yl <- cumsum(cx[,i]) - cx[,i]/2
txt <- round(x[,i], 2)
txt[txt=="0"] <- NA
 # パーセントをいれるとき
 # txt <- round(cx[,i], 2); txt[txt=="0"] <- NA
text(yl, xl, txt) # 縦型のときとxl, ylを逆にする
}