FC2ブログ

 日々の研究生活での出来事,考えたことなどのメモ。あるいはお知らせ。

  お知らせのナビゲーター   トップページ > 統計  

スポンサーサイト

-- - --/-- [--] - --:--

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

多変量行列のブートストラップ

2012 - 04/21 [Sat] - 08:48


 1変量データをブートストラップリサンプリングする方法はよく見かけるのですが,多変量データをリサンプリングする方法はあまり説明されていません。なので,ここではそれをするためのRのコードをメモしておくことにします。基本は1変量の時と変わりません。ただ行全体でブートストラップをすればよいだけです。


x<- c(1:40)
(dat <- matrix(x, 10,4))
d.boot <- matrix(0,10,4)

n<-1
for(i in sample(10,replace=T)){
d.boot[n,] <- dat[i,]
n <- n+1
}
d.boot


 上記を実行すると,以下のようになります。(ただし,乱数による結果なので,実行するたびに並び順は異なります。)

リサンプリング前

[,1] [,2] [,3] [,4]
[1,] 1 11 21 31
[2,] 2 12 22 32
[3,] 3 13 23 33
[4,] 4 14 24 34
[5,] 5 15 25 35
[6,] 6 16 26 36
[7,] 7 17 27 37
[8,] 8 18 28 38
[9,] 9 19 29 39
[10,] 10 20 30 40


リサンプリング後

[,1] [,2] [,3] [,4]
[1,] 4 14 24 34
[2,] 3 13 23 33
[3,] 8 18 28 38
[4,] 10 20 30 40
[5,] 3 13 23 33
[6,] 2 12 22 32
[7,] 8 18 28 38
[8,] 10 20 30 40
[9,] 5 15 25 35
[10,] 8 18 28 38

スポンサーサイト

ネットワーク分析

2011 - 09/10 [Sat] - 08:50

ホームページの方に社会ネットワーク分析(SNA)の解説をアップしました。
もしよろしければご活用下さい。

  UCA-Works

因子分析法(データ行列の標準化)

2010 - 05/17 [Mon] - 19:45

 前回に引き続き,芝(1979)の因子分析法に掲載されているデータを使って,今度はデータ行列の標準化を行ってみます。ご存じの通り,標準化とは各変数において,平均0,分散1となるように補正することです。データを標準化することのメリットは単位に依存することなく,各値を直接比較できるようになることです。

 方法は至って簡単で,データ行列の各値から平均値を引き,標準偏差で割ってやれば値が求まります。Rを使って標準化したいなら『scale()』という関数を利用すれば一瞬で求まります。例えば以下のようにします。


# データの読み込み
X <- read.fortran('shiba1979.txt', c("9F2.0"), header=F,skip=1)
# データの標準化(関数を使う)
scale(X)

 さて,ここで,便利な関数を使わずに,データ行列をそのまま使って,行列計算で標準化する方法を見てみます。標準化されたデータ行列をZとすると,Zは以下のような計算式で求まります。
    
ここで,N は被験者数であり,I はN×Nの単位行列,1は全要素が1のN×1ベクトル,Σは各変数の標準偏差を対角成分に持つ行列です。なので,これを素直にRで書くと,以下のようになると思います。

# データの標準化(行列計算で求める)
# 被験者数を取り出す
N <- nrow(X)

# N × Nの単位行列を作る
I <- diag(1, N)

# N × 1の1ベクトルを作る
one <- matrix(1, N, 1)

# 各変数の標準偏差を対角成分に持つ行列を作る
Sigma <- diag(round(apply(X,2,sd),3))

# 標準化された行列を求める
X <- as.matrix(X)
Z <- (I-((1/N) * (one %*% t(one)))) %*% X %*% solve(Sigma)
round(Z, 3)

 すると以下のような結果が出力されました。

> round(Z, 3)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1.742 1.827 1.705 0.121 -0.285 0.158 0.460 0.583 1.592
[2,] 1.533 2.184 1.262 -0.351 0.482 0.505 -0.554 0.189 -1.190
[3,] 1.481 0.696 1.135 -0.148 -0.558 0.678 1.007 -0.075 1.049
[4,] 1.272 1.113 0.249 0.189 0.700 0.736 1.241 0.847 0.914
[5,] 1.272 1.648 0.629 -1.093 -1.434 -0.651 0.070 1.110 -0.036
[6,] 0.854 -0.434 1.072 -0.216 0.865 1.488 0.148 1.768 1.389
[7,] 0.593 0.161 0.439 -0.688 -0.011 -0.536 -0.788 -0.798 -1.122
[8,] 0.228 0.161 -0.511 -0.283 -1.817 -1.866 0.695 1.373 0.303
[9,] 0.176 -0.315 -0.068 -1.700 -0.832 -1.172 -1.959 -1.522 -1.258
[10,] 0.071 -0.732 -0.891 -1.767 -1.270 -1.750 -1.413 -1.588 -1.393
[11,] -0.294 0.101 0.122 -0.283 -1.434 -1.519 0.382 0.123 -0.986
[12,] -0.398 -0.196 -0.131 -1.498 -0.285 -1.172 0.226 0.583 0.100
[13,] 0.854 -0.137 1.705 -0.621 -0.285 0.100 -0.632 0.978 0.642
[14,] -1.599 -1.327 -1.207 -1.565 -0.832 0.042 -0.008 -1.259 -0.579
[15,] 1.011 0.518 1.135 1.268 1.138 0.563 -0.554 -0.601 -0.986
[16,] 0.959 0.696 1.389 0.998 1.686 0.216 -0.008 1.176 0.100
[17,] 0.437 0.815 0.945 1.808 0.372 1.893 -0.398 -1.390 -1.461
[18,] 0.124 0.875 -0.574 1.673 2.178 1.257 0.851 1.110 -0.104
[19,] -0.085 0.577 0.629 0.863 1.357 0.447 2.099 0.912 2.067
[20,] -0.242 -0.613 -1.144 0.324 0.919 0.505 -1.647 -1.785 -0.579
[21,] -0.503 0.458 -0.511 1.673 -0.558 0.736 0.226 -0.206 0.642
[22,] -0.607 -0.732 -1.207 0.796 0.263 0.794 -1.022 -1.325 -1.732
[23,] -0.659 -1.327 -0.891 0.256 0.263 0.389 1.943 0.057 1.389
[24,] -0.868 0.042 -0.194 0.324 0.700 -0.189 -0.398 -0.206 0.235
[25,] -0.868 -0.256 0.122 0.594 0.591 1.257 -0.476 -0.733 0.642
[26,] -0.973 -0.851 -1.903 -0.216 0.044 0.447 -0.476 0.189 -0.579
[27,] -1.129 -1.506 -1.524 0.796 0.865 -0.073 1.007 1.241 0.235
[28,] -1.286 -0.434 -0.257 -0.890 -1.270 -0.941 1.241 0.057 0.371
[29,] -1.390 -0.970 -0.574 -0.823 -0.777 -0.709 -1.413 -0.206 0.642
[30,] -1.703 -2.041 -0.954 0.459 -0.777 -1.634 0.148 -0.601 -0.308

 『scale()』関数で求めた結果と一致していることが確認できると思います。

因子分析法(基礎統計量)

2010 - 05/16 [Sun] - 11:35

 因子分析の歴史的名著の1つに『芝 祐順 (著)  1979 因子分析法,東京大学出版会』という書籍があります。この本は今ではすでに絶版で,いまAmazonで古本の情報を調べたところ『62000円』というプレミア価格がついていました。

 なのですが,最近,訳あって,新品同然のキレイさで,しかも定価で手に入れることができました。そういうわけで,最初からゆっくりと読み始めています。この本は因子分析を数理的理論から理解するにはかなり良いですね。

 ところで,この本の最初の方のページに以下のような人工データがあります。


V1V2V3V4V5V6V7V8V9
808475475055536074
769068406461405433
756566434564605066
717252486865636464
718158292941486850
634665427178497871
585655355543373934
515640412220567255
504847204032222832
484134193222292730
415550412926525336
395046235032506052
635175365054396660
163129224053473242
666266647662404236
656570608656476952
556763726285423029
496839709574586849
456358588060746581
424330507261262442
376140704565504860
354129576066343125
343134496059725271
305445506849424854
304950546674414060
283918425660415442
252824577151607054
224644323236635256
203739334140294860
141933524124494246

 これをテキストで保存してRで読み込み,基礎統計量を求めてみます。

# データの読み込み
X <- read.fortran('shiba1979.txt', c("9F2.0"), header=F,skip=1)

# 基本統計量を求める
bstat<-function(X){
i<-ncol(X)
mat<-matrix(0,6,i)
rownames(mat)<-c("平均","中央値","分散","標準偏差","最小値","最大値")
colnames(mat)<-colnames(X)
mat[1,]<-round(apply(X,2,mean),3)
mat[2,]<-round(apply(X,2,median),3)
mat[3,]<-round(apply(X,2,var),3)
mat[4,]<-round(apply(X,2,sd),3)
mat[5,]<-round(apply(X,2,min),3)
mat[6,]<-round(apply(X,2,max),3)
mat<-t(mat)
data.frame(mat)
}

bstat(X)

 すると以下のような結果が得られます。

> bstat(X)
平均 中央値 分散 標準偏差 最小値 最大値
V1 46.633 46.5 366.999 19.157 14 80
V2 53.300 52.5 282.424 16.805 19 90
V3 48.067 46.5 249.513 15.796 18 75
V4 45.200 45.0 219.752 14.824 19 72
V5 55.200 55.5 333.890 18.273 22 95
V6 52.267 55.5 299.099 17.294 20 85
V7 47.100 47.5 164.162 12.813 22 74
V8 51.133 52.0 231.016 15.199 24 78
V9 50.533 52.0 217.223 14.738 25 81

 値を見てみると,標準偏差の結果が教科書と若干違っています。これは,Rのsdという関数が偏差平方和をN-1で割っているのに対して,教科書ではNで割っているからです。被験者数であるNの値がもっと大きければさほど影響はないのですが,今回はN=30なのでこのような大きな差異がでているのですね。

今日はゼミの日

2009 - 04/16 [Thu] - 00:35

4月15日 18:00~

 今日はゼミの日。前半は検定力分析について。後半は線形代数についてです。

 まず検定力分析についてですが,今日やったことは検定力分析を行えるRの関数の説明が中心です。テーマは『相関係数の検定』『カイ二乗検定と適合度検定』『一元配置の分散分析』『比率の差の検定』です。

 このうち僕は『比率の差』に関することで発表しました。しかしこの内容,googleで検索をかけても日本語で説明したサイトは皆無。と,言うわけで,ここでちょっとだけまとめておこうかと・・・

 比率の差で検定力分析を行う場合,必要となる指標は有意水準α,標本数n,そして効果量hの3つです。

 まず,2つの異なる母集団からそれぞれランダムに標本が得られたとき,そのサンプルサイズが等しかったとします。また,母集団Aにおける標本比率をP1,母集団Bにおける標本比率をP2とします。このときの効果量hは以下の式によって構成されます。



 また,このときに利用するRの関数はpwrパッケージの中にある『pwr.2p.test()』です。例えば,以下のような問題を考えてみましょう。

問題 
 今,2つの集団(東京と大阪)からそれぞれ80人ずつ標本が得られたとします。このとき花粉症だった人の比率を計算するとP(東京)=0.67,p(大阪)=0.50でした。有意水準α=0.05として比率の差における検定力を求めてみよう。

解答
 まずはじめに,効果量の値を求める必要があります。そのためRで以下のようなスクリプトを書き,実行します。ここで『asin()』は『arcsin』のことを意味しています。


> (phi東京 <- 2*asin(sqrt(.67)))
[1] 1.917713
> (phi大阪 <- 2*asin(sqrt(.50)))
[1] 1.570796
> (h <- phi東京 - phi大阪)
[1] 0.3469169


次に,検定力は以下のようにして求めます.


> pwr.2p.test(h=h,n=80,sig.level=0.05)

Difference of proportion power calculation for binomial distribution (arcsine transformation)

h = 0.3469169
n = 80
sig.level = 0.05
power = 0.5925747
alternative = two.sided

NOTE: same sample sizes


 このように検定力は0.59と求まっています。簡単ですね。では2つの母集団で標本数が異なる場合はどうすればよいかというと,『pwr.2p2n.test()』を使えば良いだけです。

 そんなこんなで,前半のゼミが終わると,すぐに後半の線形代数のゼミに突入です。今日は翻訳初日というわけで,3節訳した時点で終了。しかし,さすが先輩,見事な翻訳をしておられます。来週は僕が当番の日なので,先輩方に負けないようにがんばろうと思います。

 ゼミ終了後,O先輩と駅近くに新しくできたラーメン屋に行きました。食べたのは北海道味噌ラーメン。・・・味は普通でした。味噌ラーメンと言えば,やはり『純連』が思い浮かびますが,これを超えるような味噌ラーメンにはなかなか出会えないですねぇ。

 | HOME |  »

プロフィール

管理者

 専門は心理統計学で構造方程式モデリングを中心に探索的分析法の研究をしています。

 ・著書:5冊
 ・翻訳:2冊
 ・査読論文:4本
 ・紀要論文:2本
 ・報告書等:2本
 ・国際発表:5回
 ・国内発表:12回

最近の記事+コメント

カテゴリ

全記事一覧表示

過去の記事

カウンタ

現在の閲覧者数

カレンダー

07 | 2018/08 | 09
- - - 1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30 31 -

ブログ内検索

リンク

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。