goo blog サービス終了のお知らせ 

Memorandums

知覚・認知心理学の研究と教育をめぐる凡庸な日々の覚書

list() data.frame()

2008-11-29 | R
データ整理に必要なベクトルや行列の操作のうち、list()、data.frame()は基礎的なもの。以下の解説などは簡潔でわかりやすいだろう。

フリーソフトRをもちいた生態学データの加工と統計解析 I.
基礎編
宮本和樹(森林総研関西森林生態研究グループ)
http://www001.upp.so-net.ne.jp/ito-hi/stat/R1.html

R-tips
39. データフレーム事始
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/39.html
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

ANOVA 被験者内2要因 by R: 再掲載

2008-11-18 | R
 以前にもとりあげた被験者内2要因をRで実行する例だが、SPSSで「反復測定」を用いない例と比較できるように示してみた。対応のあるデータの場合でも「反復測定」に応じたデータ入力は、やはり面倒だろう。実験結果をファイルに記録する場合、測定値(反応)を一列に記録することが多いからだ。反応(従属変数)を1列に記録し(1変量)、誤差項の指定に注意すればよい。データファイルの作成、計算の実行、出力のいずれもSPSSより簡潔・容易だと思われる。
 なお、SPSS形式のデータファイル(.sav)の読み込みも示した。

cf.
被験者内要因が 2つの場合 (1変量)(星野祐司氏 ( 立命館大学))
http://www.psy.ritsumei.ac.jp/~hoshino/spss/anova02x.html

データは同じく森・吉田(1990, pp. 116-121)からで、要因a が 2水準,要因b が 4水準。

> x<-read.spss("anova02x.sav")
> x
$FACTOR_A
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

$FACTOR_B
[1] 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4

$SUBJECT
[1] 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5

$DATA
[1] 3 4 6 5 3 3 6 7 1 4 6 8 3 5 4 7 5 7 8 9 3 2 3 2 5 6 2 3 2 3 3 3 4 6 6 4 6 4 5 6

> attach(x)
> A <- factor(FACTOR_A)
> B <- factor(FACTOR_B)
> S <- factor(SUBJECT)
> aov.ex=aov(DATA~(A*B)+Error(S/(A*B)),x)
> summary(aov.ex)

Error: S
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 38.150 9.538

Error: S:A
Df Sum Sq Mean Sq F value Pr(>F)
A 1 16.9000 16.9000 8.0958 0.04662 *
Residuals 4 8.3500 2.0875
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: S:B
Df Sum Sq Mean Sq F value Pr(>F)
B 3 19.7000 6.5667 6.0383 0.00952 **
Residuals 12 13.0500 1.0875
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: S:A:B
Df Sum Sq Mean Sq F value Pr(>F)
A:B 3 30.5000 10.1667 7.0725 0.005413 **
Residuals 12 17.2500 1.4375
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

 単純主効果の分析は、以下を参考(独立な場合)に水準ごとのデータセットを用意し、上の分析を応用できる。
Rによる単純主効果の分析
 SPSSの場合は、以下を参考に。

http://www.psy.ritsumei.ac.jp/~hoshino/spss/simple02x.html
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rで多変量分散分析

2008-08-26 | R
manova()
を利用し、aovと同様に以下で結果をとりだす。
---
factor をX
従属変数をY <- cbind(y1,y2,y3) などとして結合し、
---
fit <- manova(Y ~ X)
summary.aov(fit) # univariate ANOVA tables
summary(fit, test="Wilks") # ANOVA table of Wilks' lambda
summary(fit) # same F statistics as single-df terms

Wilks のlambdaについては、以下などを参照。

http://aoki2.si.gunma-u.ac.jp/lecture/Wilks/wilks2.html

ただ、心理学では複数の従属変数をまとめて評価するよりも、それらの関係や比較が問題となることが多いように思われる。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

AIC( )とstepAIC( )

2008-07-28 | R
AIC()とstepAIC()のAIC値の比較
 cf.
Rによる多変量解析入門(下平英寿氏)
 講義資料8 「モデル選択」(2004/10/12版) pp.6-7.
http://www.is.titech.ac.jp/~shimo/class/dk2005/san08.pdf

そもそも、
AIC=-2×L(θ)+2×dim θ
 パラメタベクトル θ,その次元 dim θ
 対数尤度L(θ)
と定義される。
重回帰では、
θ = (β0 , β1 , . . . , βp , σ), dim θ = p + 1
L(θ) = -n /2 (1 + log(2π σ^2))
σ^2=Σe^2/n
であるから、
fit を線形回帰lmのobject とすると、AIC()によるAICは、
> n*(1+log(2*pi*sum(resid(fit)^2)/n)) + 2*(p+1) # AIC() の計算する AIC 値の定義
AIC = n ×(1 + log(2π σ^2 ) )+ 2 × (p + 1)
となる。

一方、stepAIC(あるいはstep() )ではextractAIC()を利用しており、これは
n log( σ^2 ) + 2p
つまり、
AIC - n × (1 + log 2π) - 2
を出力する。
> n*log(sum(resid(fit)^2)/n) + 2*p # stepAIC() の計算する AIC 値の定義




コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rで作図:物理量を横軸にとる場合

2008-07-20 | R
実験データで要因の水準別に測定値をプロットする場合、Rではデフォルトで横軸に要因の水準(名義あるいは順序尺度)をとった箱ヒゲ図を作成する。これは独立変数をfactor()あるいはordered()で要因化することによる。しかし心理学では独立変数の値(間隔尺度)を直接横軸に、水準ごとの従属変数の平均値を縦軸にとることも多い(精神物理学的測定など)。この場合、たとえば独立変数S、従属変数Rのほかに水準のラベルLを列に追加して、
y<-tapply(R,L,mean)
x<-tapply(S,L,mean)
で機械的に両変数の水準別の平均値をもとめ、
plot(x,y)
とすれば図示できる。

なお、単に水準別の平均値を図示したければ、
plot(tapply(R,L,mean))
でよい。plot() のオプション指定やtapply()の使い方は、Rヘルプか、
R-tips
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html
などを参照(pdf版あり)。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rで項目反応理論

2008-07-19 | R
4件法データに項目応答理論を適用した例。
http://ir.iwate-u.ac.jp/dspace/bitstream/10140/1865/1/arfe-v67p81-94.pdf


計算はRを利用し、IRTには以下のパッケージltm を使用している。

ltm: Latent Trait Models under IRT
by Dimitris Rizopoulos
This package provides a flexible framework for Item Response Theory analyses for dichotomous and polytomous data under a Marginal Maximum Likelihood approach.
http://wiki.r-project.org/rwiki/doku.php?id=packages:cran:ltm
多数のサンプルデータと詳細な分析例あり。
他にMASS, msm, mvtnorm, polycor が必要。

cf.
Dimitris Rizopoulos (2006).ltm: An R Package for Latent Variable Modeling and Item Response Theory Analyses. Journal of Statistical Software, Volume 17, Issue 5.

http://www.jstatsoft.org/v17/a5/paper
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rで共分散構造分析

2008-07-19 | R
パッケージ sem の使用法とAmosのようなグラフィカルモデリングの例が以下に。Windowsでは日本語処理に要工夫か。

Rで共分散構造分析・構造方程式モデル


図示には Graphvizを使用(フリー)。ファイルの拡張子は.dot
http://www.graphviz.org/
出力はたいへん美しいので、そのまま原稿に使うことが可能。
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rで重回帰

2008-07-14 | R
summary()では出力されない標準化偏回帰係数を得る方法。

http://phi.med.gunma-u.ac.jp/swtips/R.html#MISC

なおAICは、
AIC( )
で簡単に取り出すことができる。また、変数選択ためにはstep()で各モデルの比較を行うことができる。

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

R on Windows

2008-07-10 | R
R Console 内の日本語が文字化けするならば、編集->GUIプリファレンスを開いて Font を MS Gothic などに変更。
SPSSのようなGUIを使いたければ、パッケージ関数を利用できる。たとえば、
Rcmdr: R Commander
などがある。
利用法は、以下など参照。

なかなかRの世界に踏み込めないでいるあなたへ...
http://hashi8.hp.infoseek.co.jp/tips/Rintro.html

パッケージ -> CRAN ミラーサイトの設定
の後、
パッケージ ->パッケージのインストール
パッケージ->パッケージの読み込み 
でRcmdrを選ぶと、別窓で「Rコマンダー」が開く。

Rコマンダーを用いたRの学習には、以下などがある。
慶応大/古谷研究室

Rの本領はこのようなパッケージの豊富さにある。
CRANパッケージリスト/RWiki

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rで因子分析

2008-07-05 | R
Rには最尤法による因子分析の関数factanal が実装されている。

factanal(x,factors,rotation,scores,...)

xはデータセット
factorsは因子の数
rotationはバリマックス(varimax)回転とプロマックス(promax)回転を指定できる
 デフォルトは直交回転 “varimax”
scoresは因子得点を求める方法で回帰方法(regression)とバートレット法 (Bartlett)から選択
 デフォルトは “none”で因子得点を出力しない

出力のうち、
カイ2乗検定統計量は元のデータの分散と指定した共通因子のモデルに基づいて求めたデータの分散との間の有意差に関する検定

factanalは因子数を指定する点が不便なので、これをもとめるプログラムとして以下のものがある。
青木研究室(群馬大)
http://aoki2.si.gunma-u.ac.jp/R/factanal2.html

渡邊研究室(慶応大)
http://web.sfc.keio.ac.jp/~watanabe/adstat10.htm#1

分析を実行するためのデータ例として、職業適性能力検査データが以下で利用できる。
http://www1.doshisha.ac.jp/~mjin/data/

cf.
http://www1.doshisha.ac.jp/~mjin/R/0508_25.pdf
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rと非線形回帰分析

2008-06-30 | R
自由に関数式を指定することができる非線形回帰分析の関数nlsの使い方。

nls(formula,data,start,trace)

具体例は下記に。
http://www1.doshisha.ac.jp/~mjin/R/0411_16.pdf
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rでベイズ統計学

2008-06-30 | R
パッケージの一括インストール
install.packages("ctv")
library(ctv)
install.views("Bayesian")


cf.
RjpWiki

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rによる心理統計再入門10:分散分析 単純主効果の検定

2008-06-25 | R
2要因(被験者間)で交互作用が有意な場合の単純主効果の検定をRで行った例。誤差項の選択(各水準別)がSPSS(全体)の場合と異なる。

dataは下記から。
小塩研究室(中部大)
http://psy.isc.chubu.ac.jp/%7Eoshiolab/teaching_folder/datakaiseki_folder/05_folder/da05_02.html

failure(f1,f2,f3)、perfection(p1,p2)の2要因被験者間計画。従属変数はdepression。

1) データの読み込み
> data.ex=read.csv("data5_1.csv",header=T)
> data.ex
sub failure perfection depression
1 1 f1 p1 10
2 2 f1 p1 13
3 3 f1 p1 21
4 4 f1 p1 16
5 5 f1 p2 16
6 6 f1 p2 19
7 7 f1 p2 13
8 8 f1 p2 8
9 9 f2 p1 15
10 0 f2 p1 16
11 11 f2 p1 12
12 12 f2 p1 15
13 13 f2 p2 21
14 14 f2 p2 23
15 15 f2 p2 16
16 16 f2 p2 19
17 17 f3 p1 21
18 18 f3 p1 14
19 19 f3 p1 24
20 20 f3 p1 20
21 21 f3 p2 31
22 22 f3 p2 36
23 23 f3 p2 24
24 24 f3 p2 34

2) 2要因の分散分析
> aov.ex=aov(depression ~ failure*perfection, data.ex)
> summary(aov.ex)
Df Sum Sq Mean Sq F value Pr(>F)
failure 2 528.08 264.04 15.6727 0.0001143 ***
perfection 1 165.37 165.37 9.8162 0.0057477 **
failure:perfection 2 156.25 78.13 4.6373 0.0237486 *
Residuals 18 303.25 16.85
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

交互作用が有意なので、まずperfectionの水準別にfailureの効果(単純主効果)を検定する。

3) 下位のデータセットの準備と単純主効果の検定
> attach(data.ex)
> failure <- factor(failure)
> perfection <- factor(perfection)

要因perfection の水準p1のみ取り出す。
> data.ex[data.ex$perfection=="p1",]
sub failure perfection depression
1 1 f1 p1 10
2 2 f1 p1 13
3 3 f1 p1 21
4 4 f1 p1 16
9 9 f2 p1 15
10 0 f2 p1 16
11 11 f2 p1 12
12 12 f2 p1 15
17 17 f3 p1 21
18 18 f3 p1 14
19 19 f3 p1 24
20 20 f3 p1 20

> perf1 <- data.ex[data.ex$perfection=="p1",]
> summary(aov(perf1$depression~factor(perf1$failure)))
Df Sum Sq Mean Sq F value Pr(>F)
factor(perf1$failure) 2 67.167 33.583 2.3659 0.1494
Residuals 9 127.750 14.194

次に水準p2での検定。
> perf2 <- data.ex[data.ex$perfection=="p2",]
> summary(aov(perf2$depression~factor(perf2$failure)))
Df Sum Sq Mean Sq F value Pr(>F)
factor(perf2$failure) 2 617.17 308.58 15.825 0.001131 **
Residuals 9 175.50 19.50
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p2ではfailureの効果が有意なので多重比較をおこなう。
> pairwise.t.test(perf2$depression, factor(perf2$failure), p.adj = "bonf")

Pairwise comparisons using t tests with pooled SD

data: perf2$depression and factor(perf2$failure)

f1 f2
f2 0.2960 -
f3 0.0011 0.0152

P value adjustment method: bonferroni


4) さらにfailure の水準別にperfectionの効果を検定。
f1の場合。
> data.ex[data.ex$failure=="f1",]
sub failure perfection depression
1 1 f1 p1 10
2 2 f1 p1 13
3 3 f1 p1 21
4 4 f1 p1 16
5 5 f1 p2 16
6 6 f1 p2 19
7 7 f1 p2 13
8 8 f1 p2 8
> fail1 <- data.ex[data.ex$failure=="f1",]
> t.test(fail1$depression~factor(fail1$perfection), var.equal=T)

Two Sample t-test

data: fail1$depression by factor(fail1$perfection)
t = 0.3015, df = 6, p-value = 0.7732
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-7.115489 9.115489
sample estimates:
mean in group p1 mean in group p2
15 14

f2の検定。
> fail2 <- data.ex[data.ex$failure=="f2",]
> t.test(fail2$depression~factor(fail2$perfection), var.equal=T)

Two Sample t-test

data: fail2$depression by factor(fail2$perfection)
t = -3.0417, df = 6, p-value = 0.02275
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-9.473434 -1.026566
sample estimates:
mean in group p1 mean in group p2
14.50 19.75

最後にf3での検定。
> fail3 <- data.ex[data.ex$failure=="f3",]
> t.test(fail3$depression~factor(fail3$perfection), var.equal=T)

Two Sample t-test

data: fail3$depression by factor(fail3$perfection)
t = -3.4223, df = 6, p-value = 0.01410
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-19.722376 -3.277624
sample estimates:
mean in group p1 mean in group p2
19.75 31.25

>





cf.
http://www.educ.kyoto-u.ac.jp/cogpsy/personal/Kusumi/datasem07/hand.txt


コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

多重比較

2008-06-02 | R
被験者内計画の際の多重比較について。

> pairwaise.t.test

Usage

pairwise.t.test(x, g, p.adjust.method = p.adjust.methods,
pool.sd = TRUE, ...)


> help(p.adjust)
p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
# "fdr", "none")

:

Details

The adjustment methods include the Bonferroni correction ("bonferroni") in which the p-values are multiplied by the number of comparisons. Less conservative corrections are also included by Holm (1979) ("holm"), Hochberg (1988) ("hochberg"), Hommel (1988) ("hommel"), Benjamini & Hochberg (1995) ("BH"), and Benjamini & Yekutieli (2001) ("BY"), respectively. A pass-through option ("none") is also included. The set of methods are contained in the p.adjust.methods vector for the benefit of methods that need to have the method as an option and pass it on to p.adjust.

The first four methods are designed to give strong control of the family wise error rate. There seems no reason to use the unmodified Bonferroni correction because it is dominated by Holm's method, which is also valid under arbitrary assumptions.

:

References

Holm, S. (1979). A simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, 6, 65–70.

Hommel, G. (1988). A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika, 75, 383–386.

Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of significance. Biometrika, 75, 800–803.

Shaffer, J. P. (1995). Multiple hypothesis testing. Annual Review of Psychology, 46, 561–576. (An excellent review of the area.)



コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Rによる心理統計再入門9:分散分析 被験者内2要因

2008-06-02 | R
前述に引き続き、心理学で最も多く使われる計画のひとつとして2要因被験者内計画の分析を取り上げる。
前と同様、岩原(1965)p.308の例(被験者要因をいれて3要因計画とみなす)を実行してみる。実験要因はA(2水準)とB(3水準)。Cが被験者(4人)で、セルの大きさは2。

> data.ex=read.csv("Iwahara24-4.csv",header=T)
> data.ex
A B C REP SCORE
1 a1 b1 c1 1 2
2 a1 b1 c1 2 1
3 a1 b1 c2 1 0
4 a1 b1 c2 2 2
5 a1 b1 c3 1 2
6 a1 b1 c3 2 0
7 a1 b1 c4 1 1
8 a1 b1 c4 2 1
9 a1 b2 c1 1 4
10 a1 b2 c1 2 5
11 a1 b2 c2 1 3
12 a1 b2 c2 2 4
13 a1 b2 c3 1 5
14 a1 b2 c3 2 4
15 a1 b2 c4 1 4
16 a1 b2 c4 2 3
17 a1 b3 c1 1 0
18 a1 b3 c1 2 1
19 a1 b3 c2 1 1
20 a1 b3 c2 2 0
21 a1 b3 c3 1 0
22 a1 b3 c3 2 0
23 a1 b3 c4 1 1
24 a1 b3 c4 2 1
25 a2 b1 c1 1 0
26 a2 b1 c1 2 0
27 a2 b1 c2 1 0
28 a2 b1 c2 2 2
29 a2 b1 c3 1 2
30 a2 b1 c3 2 0
31 a2 b1 c4 1 2
32 a2 b1 c4 2 2
33 a2 b2 c1 1 4
34 a2 b2 c1 2 6
35 a2 b2 c2 1 4
36 a2 b2 c2 2 5
37 a2 b2 c3 1 6
38 a2 b2 c3 2 4
39 a2 b2 c4 1 2
40 a2 b2 c4 2 3
41 a2 b3 c1 1 0
42 a2 b3 c1 2 2
43 a2 b3 c2 1 1
44 a2 b3 c2 2 0
45 a2 b3 c3 1 0
46 a2 b3 c3 2 2
47 a2 b3 c4 1 0
48 a2 b3 c4 2 1
> aov.ex=aov(SCORE~A*B+Error(C/(A*B)),data.ex)
> summary(aov.ex)

Error: C
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 3 1.06250 0.35417

Error: C:A
Df Sum Sq Mean Sq F value Pr(>F)
A 1 0.18750 0.18750 0.5294 0.5195
Residuals 3 1.06250 0.35417

Error: C:B
Df Sum Sq Mean Sq F value Pr(>F)
B 2 116.375 58.188 40.478 0.0003285 ***
Residuals 6 8.625 1.438
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: C:A:B
Df Sum Sq Mean Sq F value Pr(>F)
A:B 2 0.3750 0.1875 0.2 0.824
Residuals 6 5.6250 0.9375

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 24 21.5000 0.8958
>


Error 内でError(C/(A*B)) のように(A*B) の()を忘れないように。
岩原(1965) p.309 の表24.4F はA,B のF値をもとめる際の分母をまちがえているようだ。 直前の表(24.4E)にあるようにMSA/MSAC でAのF値をもとめるので、上に計算のようにAは0.5294であろう。同様にBのF値は40.478と思われる。
C(被験者)は変量モデル(ランダム効果)であるから、検定の必要はない。

Bについて多重比較をBonferroni法でおこなうと、以下のようにb1とb2、b2とb3に有意差が認められる。

> pairwise.t.test(data.ex$SCORE,data.ex$B,p.adj="bonferroni")

Pairwise comparisons using t tests with pooled SD

data: data.ex$SCORE and data.ex$B

b1 b2
b2 1.1e-11 -
b3 0.56 1.7e-13

P value adjustment method: bonferroni
>

References
岩原信九郎 (1965). 教育と心理のための推計学 新訂版 日本文化科学社
コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする