R GUI Projects
複数(被験者別)のデータファイル(.txt スペース区切り)を結合し,さらにその一部(観測値)を抽出する。ここでの例は,9変数(V1,,,V9)で,正反応(V9=1)のみを取り出し,かつ反応時間(V8)が200ms以下と3000ms以上を除外する。
data1 = read.table("subject1.txt",sep=" ",header=F)
data2 = read.table("subject2.txt",sep=" ",header=F)
data12= rbind(data1,data2)
newdata12 = subset(data12, V9=="1" & V8 < 3000 & V8 > 200)
3つ以上のデータフレームの結合も同様
data123 = rbind(data1,data2,data3)
さらに,保存
write.table ( data123, "data123.txt" )
なお,データフレームに代入する = は <ー と同じ。
以下のサイトのデータを使い、Rのsemパッケージで検証的因子分析をおこなった結果を、GraphVizを使って図示してみた。
> path.diagram(MODEL1.sem, "ex1", ignore.double=FALSE, edge.labels="values", digits=2, standardize=TRUE)
Running dot -Tpdf -o ex1.pdf ex1.dot
ファイル名ex1で出力→ ex1.dot → ex1.pdf を自動生成
cf.
GraphViz Download
irisデータを使って、目的とする変数の基本統計量をグループ別にもとめる。
by(目的変数,グループ化変数,summary)
> iris
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
:
148 6.5 3.0 5.2 2.0 virginica
149 6.2 3.4 5.4 2.3 virginica
150 5.9 3.0 5.1 1.8 virginica
> by ( iris[1], iris[5], summary )
Species: setosa
Sepal.Length
Min. :4.300
1st Qu.:4.800
Median :5.000
Mean :5.006
3rd Qu.:5.200
Max. :5.800
------------------------------------------------------------
Species: versicolor
Sepal.Length
Min. :4.900
1st Qu.:5.600
Median :5.900
Mean :5.936
3rd Qu.:6.300
Max. :7.000
------------------------------------------------------------
Species: virginica
Sepal.Length
Min. :4.900
1st Qu.:6.225
Median :6.500
Mean :6.588
3rd Qu.:6.900
Max. :7.900
SDならば,
> by ( iris[1], iris[5], sd)
Species: setosa
[1] 0.3524897
------------------------------------------------------------
Species: versicolor
[1] 0.5161711
------------------------------------------------------------
Species: virginica
[1] 0.6358796
一般に,複数の要因の水準の組み合わせ(たとえば:データxの2から4列目)ごとに求めるならば,グループ化変数をx[2:4]のように指定する。
cf. tapply()
Rで数量化類(双対尺度法、コレスポンデンス分析)を行う。
cf. 群馬大学 青木研究室
http://aoki2.si.gunma-u.ac.jp/R/qt3.html
アイテムデータ形式のデータを用意。
> dat
+Cat.A=c("a", "c", "b", "a", "b", "b"),
+ Cat.B=c("b", "b", "c", "b", "a", "c"),
+ Cat.C=c("c", "a", "a", "c", "b", "b")
+ )
>
> dat
Cat.A Cat.B Cat.C
1 a b c
2 c b a
3 b c a
4 a b c
5 b a b
6 b c b
>
ダミー変数を使って、カテゴリーデータに変換。
> source("http://aoki2.si.gunma-u.ac.jp/R/src/make.dummy.R", encoding="euc-jp")
> cat.dat
> cat.dat
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1 0 0 0 1 0 0 0 1
[2,] 0 0 1 0 1 0 1 0 0
[3,] 0 1 0 0 0 1 1 0 0
[4,] 1 0 0 0 1 0 0 0 1
[5,] 0 1 0 1 0 0 0 1 0
[6,] 0 1 0 0 0 1 0 1 0
数量化類を実施。
> source("http://aoki2.si.gunma-u.ac.jp/R/src/qt3.R", encoding="euc-jp")
> a
> summary(a)
$Eigen.value
解1 解2 解3 解4
0.93253 0.59861 0.36192 0.10694
$Correlation.coefficient
解1 解2 解3 解4
0.96568 0.77370 0.60160 0.32701
$Category.score
解1 解2 解3 解4
Cat.A.a -1.24668 0.90967 -0.37059 0.22888
Cat.A.b 0.99434 0.18573 -0.33881 0.38433
Cat.A.c -0.48967 -2.37654 1.75761 -1.61076
Cat.B.a 1.18320 1.23348 2.56331 1.94775
Cat.B.b -0.99434 -0.18573 0.33881 -0.38433
Cat.B.c 0.89991 -0.33814 -1.78987 -0.39738
Cat.C.a 0.11411 -1.70558 -0.18806 1.47834
Cat.C.b 1.13257 0.79591 0.55865 -1.70722
Cat.C.c -1.24668 0.90967 -0.37059 0.22888
$Sample.score
解1 解2 解3 解4
#1 -1.20389 0.70381 -0.22295 0.074847
#2 -0.47287 -1.83872 1.05738 -0.526739
#3 0.69325 -0.80048 -1.28365 1.493610
#4 -1.20389 0.70381 -0.22295 0.074847
#5 1.14259 0.95434 1.54209 0.636940
#6 1.04480 0.27724 -0.86991 -1.753505
attr(,"class")
[1] "qt3"
Category.score 解1 解2 を使ってプロット。
> plot(a, pch=19, label.cex=0.5)
Sample.score の重ね書きをするには,
> plot(a, pch=19, label.cex=0.5,xlim=c(-2,2),ylim=c(-2,2),ann=F)
> par(new=T)
> plot(a, which="sample.score", label.cex=0.5, xlim=c(-2,2),ylim=c(-2,2))
(この例では,#1と#4が重なっている)
2013/11/19 サンプルスコアの重ね書きを加筆
http://personality-project.org/r/#anova
プログラム例と解説
Mixed (between and Within) designs
Example
2 between-subjects variables: Gender (2 levels) and Dosage (3 levels)
2 within-subjects variables: Task (2 levels) and Valence (3 levels)
簡単に正規分布を図示したいという院生の質問があったので。
3行目がグラフ,2行目は縦横軸の設定と縦軸のラベル,1行目は横軸の範囲。
x <- seq(-4, 4, 0.05)
plot(x,dnorm(x, mean=0, sd=1), type="n")
curve(dnorm(x, mean=0, sd=1), type="l")
一行で済ませるならば,
curve(dnorm(x, mean=0, sd=1),from=-4, to=4)
で可。
cf.
被験者内要因が 2つの場合 (1変量)(星野祐司氏 ( 立命館大学))
http://www.psy.ritsumei.ac.jp/~hoshino/spss/anova02x.html
多重比較はBonferroni を用いた。例題は森・吉田(1990, pp. 116-121)からで、要因A が 2水準,要因B が 4水準。
多重比較については,以下も参照。
入戸野 宏. (2004). 心理生理学データの分散分析. 生理心理学と精神生理学 22, 275-290.
http://home.hiroshima-u.ac.jp/nittono/ANOVA_JSPP2004.pdf
> x<-read.csv("anova02x.csv")
> x
FACTOR_A FACTOR_B SUBJECT DATA
1 1 1 1 3
2 1 2 1 4
3 1 3 1 6
4 1 4 1 5
5 1 1 2 3
6 1 2 2 3
7 1 3 2 6
8 1 4 2 7
9 1 1 3 1
10 1 2 3 4
11 1 3 3 6
12 1 4 3 8
13 1 1 4 3
14 1 2 4 5
15 1 3 4 4
16 1 4 4 7
17 1 1 5 5
18 1 2 5 7
19 1 3 5 8
20 1 4 5 9
21 2 1 1 3
22 2 2 1 2
23 2 3 1 3
24 2 4 1 2
25 2 1 2 5
26 2 2 2 6
27 2 3 2 2
28 2 4 2 3
29 2 1 3 2
30 2 2 3 3
31 2 3 3 3
32 2 4 3 3
33 2 1 4 4
34 2 2 4 6
35 2 3 4 6
36 2 4 4 4
37 2 1 5 6
38 2 2 5 4
39 2 3 5 5
40 2 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
>#記述統計
> tapply(DATA,list(A,B),length)
1 2 3 4
1 5 5 5 5
2 5 5 5 5
> tapply(DATA,list(A,B),mean)
1 2 3 4
1 3 4.6 6.0 7.2
2 4 4.2 3.8 3.6
> tapply(DATA,list(A,B),sd)
1 2 3 4
1 1.414214 1.516575 1.414214 1.483240
2 1.581139 1.788854 1.643168 1.516575
>#要因Aの水準1のデータをとりだす
> a1<-x[x$FACTOR_A=="1",]
> a1
FACTOR_A FACTOR_B SUBJECT DATA
1 1 1 1 3
2 1 2 1 4
3 1 3 1 6
4 1 4 1 5
5 1 1 2 3
6 1 2 2 3
7 1 3 2 6
8 1 4 2 7
9 1 1 3 1
10 1 2 3 4
11 1 3 3 6
12 1 4 3 8
13 1 1 4 3
14 1 2 4 5
15 1 3 4 4
16 1 4 4 7
17 1 1 5 5
18 1 2 5 7
19 1 3 5 8
20 1 4 5 9
> detach(x)
># 単純主効果
> attach(a1)
> A <- factor(FACTOR_A)
> B <- factor(FACTOR_B)
> S <- factor(SUBJECT)
> summary(aov(DATA~B+Error(S/B)))
Error: S
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 21.2 5.3
Error: S:B
Df Sum Sq Mean Sq F value Pr(>F)
B 3 49.200 16.400 15.375 0.000206 ***
Residuals 12 12.800 1.067
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>#多重比較
> pairwise.t.test(DATA,B,p.adj="bonferroni")
Pairwise comparisons using t tests with pooled SD
data: DATA and B
1 2 3
2 0.6113 - -
3 0.0299 0.8904 -
4 0.0019 0.0739 1.0000
P value adjustment method: bonferroni
> detach(a1)
以下,A2におけるBの効果,Bの各水準におけるAの効果を同様に繰り返す。
豊田秀樹 編著
データ・スクリプトファイルのダウンロード可。
http://www.tokyo-tosho.co.jp/books/ISBN978-4-489-02065-0.html