R に関して －目次－
【分析】
・ 正準相関分析
/// LifeCycleSavings
・ クラスター分析
/// USArrests
・ 因子分析
/// USArrests
・ 主成分分析
/// USArrests
・ 数量化Ⅲ類
/// PC 内データ
・ 対応分析
/// caith（髪色と瞳の色）
/// PC 内データ
【テクニック】
・ ラベルの付け方
・ 色の付け方
【正準相関分析】
#データの読み込み(ここでは R 内蔵の LifeCycleSavings を用いる)
x<-LifeCycleSavings
#変数を 2 つの群に分ける
v1<-c(2,3)
v2<-c(1,4,5)
#各群の相関係数行列 S,T と両群の共分散行列 R を計算する
S<-cor(x)[v1,v1]
T<-cor(x)[v2,v2]
R<-cor(x)[v1,v2]
#固有値問題を解くことになる
z<-eigen(solve(S) %*% R %*% solve(T) %*% t(R))
Z<-sqrt(z$values[1])
zz<-1/Z*solve(T) %*% t(R) %*% z$vectors[,1]
#ここから微妙
#u<-z$vectors[1,1]*x[,2]+z$vectors[2,1]*x[,3]
#v<-zz[1,1]*x[,1]+ zz[2,1]*x[,4]+ zz[3,1]*x[,5]
#cor(u,v)
#実は R の stats ライブラリに正準相関を求める関数が内蔵されている
x<-LifeCycleSavings
library(stats)
pop<-x[,2:3]
oec<-x[,-c(2,3)]
#生データの場合
cancor(pop,oec)
#基準化した場合
cancor(scale(pop),scale(oec))
【クラスター分析】
#データの読み込み(ここでは R 内蔵の USArrests を用いる)
#read.table("C:/Users/ogi/Documents/kekkon.csv",sep=",",row.names=1,header=T)
x<-USArrests
y<-hclust(dist(scale(x)))
plot(y)
p<-identify(y, function(k) apply(x[k,], 2, mean ) )
#各クラスターに含まれる個体の行番号ベクトル k を引数とする関数を実行し、その結果を
出力せよ、という意味になる。関数の中身は、データ x(=USArrests)から、当該クラスター
に含まれる行だけを抜き出し、列方向（変数方向）に、平均を計算せよ、という意味
エラー
以下にエラー cutree(x, k = 2:MAXCLUSTER) :
elements of 'k' must be between 1 and 19
原因：行数（サンプル数が 20 個未満だと identify ができない）
解決案：identify(z,MAXCLUSTER=19)
※上記の 19 は一例であって、その時の行数（サンプル数）を入れる。
よって…identify(z,MAXCLUSTER=nrow(x))となり上記のでいうと
p<-identify(y, function(k) apply(x[k,], 2, mean ), MAXCLUSTER=nrow(x) )
#ウォード法
hclust(x,method="ward")
#K-Means 法
kmeans(x,n)
【因子分析】
#データの読み込み(ここでは R 内蔵の USArrests を用いる)
#read.table("ank1.csv",sep=",",row.names=1,header=T)
x<-USArrests
#固有値・固有ベクトルを計算
y<-eigen(cor(x))
e<-y$values
e
#相関行列の対角成分を SMC に置き換えた行列の固有値を計算（正の数が基準）
corx<-cor(x)
diag(corx)<-diag((1-(1/solve(corx)))[1:length(x),1:length(x)])
eigen(corx)
#累積寄与率を計算
sum(e[1:2])/sum(e)
#因子負荷量の計算（ここではまだ因子の解釈は行わない）
z<-t(sqrt(y$values)*t(y$vectors))[,1:2]
#バリマックス回転＆因子の解釈
varimax(z)
#これは勝手に基準化を行っている
#基準化を行わない場合には varimax(z,normalize=FALSE)
#因子得点を計算する
f<-varimax(z)$loadings
s<-scale(x)%*%solve(cor(x))%*%f
s
#solve(x)とは行列 x の逆行列である
#因子得点のプロット（打点なし・行列名使用）
plot(s,type="n")
text(s,rownames(x))
#最後に biplot をする
biplot(s,f,xlabs=rownames(x),ylabs=colnames(x))
【主成分分析】
#データの読み込み(ここでは R 内蔵の USArrests を用いる)
#read.table("ank1.csv",sep=",",row.names=1,header=T)
x<-USArrests
#固有値・固有ベクトルを計算
y<-eigen(cor(x))
e<-y$values
e
#累積寄与率を計算
sum(e[1:2])/sum(e)
#因子負荷量の計算＆主成分の解釈
z<-t(sqrt(y$values)*t(y$vectors))[,1:2]
z
#主成分得点を計算する
s<-scale(x)%*%y$vectors[,1:2]
s
#主成分得点のプロット（
（打点なし・行列名使用）
plot(s,type="n")
text(s,rownames(x))
#実は R に主成分分析の関数は用意されている
#主成分分析型のデータ出力
#分散共分散行列から計算する場合
y<-princomp(x)
#相関係数行列から計算する場合
y<-princomp(x,cor=TRUE)
#結果の概要
#標準偏差・寄与率・累積寄与率
summary(y)
#上記＋因子負荷量
summary(y,loadings=TRUE)
#プロット
biplot(y)
【数量化Ⅲ類】
#データの読み込み(ここでは PC から読み込み)
x<- read.table("3rui.csv",sep=",",row.names=1,header=T)
#このままでは行列の形になっていないので行列にする
x<-as.matrix(x)
u<-apply(x,1,sum)
v<-apply(x,2,sum)
X<-diag(1/sqrt(u)) %*% x %*% diag(1/sqrt(v))
z<-eigen(X%*%t(X))
z
xaxis<-diag(sqrt(sum(u)/u))%*%z$vectors
#1 番目の固有値が 1 で対応する座標がすべて一点に集中する「無意味」な解
#よって 2 番目以降の固有ベクトルを考えれば良い
drink<- xaxis[,2:3]
plot(drink,type="n")
text(drink,rownames(x))
#個人マップの作成
zz<-eigen(t(X)%*%X)
member<-diag(sqrt(sum(v)/v))%*%zz$vectors[,2:3]
plot(member,type="n")
text(member,colnames(x))
#同時布置
rowd<-length(drink[,1])
rowm<-length(member[,1])
toge<-matrix(,rowd+rowm,2)
toge[1:rowd,1:2]<-drink
toge[-c(1:rowd),1:2]<-member
#布置する際には片方の色を変えるとわかりやすい
plot(toge,type="n")
text(drink,rownames(x))
text(member,colnames(x), col="red")
#biplot をする必要ならば最終行の前に library(mva)
yaxis<-cbind(diag(1/(z$values[2]*v))%*%t(x)%*%xaxis[,2],diag(1/(z$values[3]*v))%*%t(
x)%*%xaxis[,3])
biplot(xaxis[,2:3],yaxis,xlabs=rownames(x),ylabs=colnames(x))
【対応分析】
#データの読み込み(ここでは R 内蔵の caith を用いる)
#caith データを出力する際には MASS パッケージを読み込む必要がある
library(MASS)
x<-caith
m<-nrow(x)
n<-ncol(x)
z<-corresp(x,nf=min(m,n)-1)
koyu<-z$cor^2
round(koyu,3)
round(100*koyu/sum(koyu),2)
#ここで寄与率を見て、次元をどこまで見ればいいのか解釈する
#ここでは 2 次元とする
s<-2
r<-matrix(0,m+n,s)
for(i in 1:(m+n)){
for(j in 1:s){
if(i<=m){
r[i,j]<-z$rscore[i,j]}
if(i>m){
r[i,j]<-z$cscore[i-m,j]}}}
r
plot(r[,1],r[,2],main="Result_ank",type="n")
text(r[1:m,1],r[1:m,2],label=row.names(x))
text(r[(m+1):(m+n),1],r[(m+1):(m+n),2],label=colnames(x), col="red")
#髪にも瞳にも brown がありわかりづらいので色を変えた
abline(0,0,0,0)
【対応分析】～PC 内データ～
x<-read.table("ank1.csv",sep=",",row.names=1,header=T)
#m は行数、n は列数
m<-nrow(x)
n<-ncol(x)
require(MASS)
z<-corresp(x,nf=min(m,n)-1)
koyu<-z$cor^2
round(koyu,3)
round(100*koyu/sum(koyu),2)
#ここで寄与率を見て、次元をどこまで見ればいいのか解釈する
#ここでは 2 次元とする
s<-2
r<-matrix(0,m+n,s)
for(i in 1:(m+n)){
for(j in 1:s){
if(i<=m){
r[i,j]<-z$rscore[i,j]}
if(i>m){
r[i,j]<-z$cscore[i-m,j]}}}
r
plot(r[,1],r[,2],main="Result_ank",type="n")
text(r[1:m,1],r[1:m,2],label=row.names(x))
text(r[(m+1):(m+n),1],r[(m+1):(m+n),2],label=colnames(x))
abline(0,0,0,0)
【ラベルの付け方】
まず plot と text とは？
plot とは座標を入れることでグラフを作成し、打点すること。
通常○で打点されてしまうが、後から text を入れる場合は邪魔なので
type=”n”と入れることで目に見えない打点をしてくれる
text とは打点したポイントに名前を表示すること
※しかし通常、名前は入っておらず、番号が表示される
ラベルの付け方の流れ
×：各行に名前をつける→text で表示
○：text で表示する際に表示方法の名前を指定する
通常：text(r), text(r[,1],r[,2])
指定：text(r,label=○○), text(r[,1],r[,2],label=○○)
※ r.label なら【r の label 要素】だが
r,label なので【r を表示、ただし label は…】というニュアンス
label の指定の仕方（どれでも良い）
スムーズなやり方
データの読み込みの際に行名と列名を指定して読み込む
x<-read.table(“rsample.csv”,sep=”,”,row.names=1,header=T)
行名の取り出し
rownames(x) もしくは row.names(x)
列名の取り出し
colnames(x) ※列名の取り出しは col.names(x)ではできない
plot(r[,1],r[,2],type=”n”)
片方 text(r,label=rownames(x))
## text(r[,1],r[,2],label=rownames(x))
同時 text(r,label=c(rownames(x),colnames(x)))
ゴリ押し（事前用意）
rlabel<-c(”nameA”,”nameB”,”nameC”,…,”nameZ”)
plot(r[,1],r[,2],type=”n”)
text(r,label=rlabel)
## text(r[,1],r[,2],label=rlabel)
ゴリ押し（直接ぶちこむ）
plot(r[,1],r[,2],type=”n”)
text(r,label= c(”nameA”,”nameB”,”nameC”,…,”nameZ”))
## text(r[,1],r[,2],label= c(”nameA”,”nameB”,”nameC”,…,”nameZ”))
【色の付け方】
プロット後に名前をつける際に用いるテクニック
特に同時布置の場合には色を変えた方が見やすい
使用例
※因子分析などで plot 後に名前をつけるときに便利
plot(sample,type="n")
text(sample,rownames(x),col="red")
※クラスター分析の樹形図の色も変更可能（単色のみ）
plot(hclust(dist(scale(USArrests))),col="red")
※パレットを見ることもできる（以下の例ではこれを用いる）
pie(rep(1, 16),col=rainbow(16))
具体的な色指定
pie(rep(1, 16),col="red")
この red を変えれば良い
R には 657 色の名前が登録されているので適当にやればできるはず
一覧の表示コマンド
colors()
グラデーション
登録されているグラデーション名＋色数（最大 16）でそれに沿った色が付いてくれる
※クラスター分析の樹形図では最初の一色しか表示されない
因子分析・主成分分析などの text ではグラデーション確認済み
pie(rep(1, 16),col=cm.colors(16))
登録されているグラデーションには以下のようなものがある
rainbow(n)
虹色
heat.colors(n)
熱（白－黄－赤）
cm.colors(n)
シアン－マゼンタ
topo.colors(n)
topography：地質（黄－緑－青－紫）
terrain.colors(n) 地形（白－薄茶－黄－緑）
grey(0:15/16)
グレースケール（これだけ数値は 0~1 の少数）
rainbow の開始点となる色を指定（rainbow 限定）
pie(rep(1, 16),col=rainbow(10,start=.7,end=.2))
© Copyright 2024 ExpyDoc
