「R」のブログ記事一覧-Memorandums

R GUI Projects

2013-11-28 | R

R GUI Projects

http://www.sciviews.org/_rgui/

ヒストグラムと密度の推定

2013-05-03 | R

階級の設定についての工夫。

ヒストグラムと密度の推定 - RjpWiki

データの結合，抽出

2012-11-30 | R

複数（被験者別）のデータファイル（.txt スペース区切り）を結合し，さらにその一部（観測値）を抽出する。ここでの例は，９変数（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)

３つ以上のデータフレームの結合も同様

data123 = rbind(data1,data2,data3)

さらに，保存

write.table ( data123, "data123.txt" )

なお，データフレームに代入する　＝　は　＜ー　と同じ。

SEM by R

2012-06-16 | R

以下のサイトのデータを使い、Rのsemパッケージで検証的因子分析をおこなった結果を、GraphVizを使って図示してみた。

RでSEM（共分散構造分析／構造方程式モデリング）3

> 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

グループ別の基本統計量

2012-06-13 | R

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の２から４列目）ごとに求めるならば，グループ化変数をx[2:4]のように指定する。

cf. tapply()

R 　数量化3類

2011-09-02 | R

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

[,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　サンプルスコアの重ね書きを加筆

R 古典的MDS　

2011-08-18 | R

Rにあるサンプルデータ　eurodist をつかって古典的MDSを実行。

> x <-eurodist

> euro <-cmdscale(x)

> plot(euro,type="n")

> text(euro,rownames(euro))

R: ANOVA 混合計画

2011-08-06 | R

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)

正規分布の図示　by R

2011-05-25 | R

簡単に正規分布を図示したいという院生の質問があったので。

３行目がグラフ，２行目は縦横軸の設定と縦軸のラベル，１行目は横軸の範囲。

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)

で可。

CLASSICAL MDS by R

2011-02-17 | R

RによるMDS

http://www.statmethods.net/advstats/mds.html

Rによる心理統計入門１３：被験者内ANOVA

2010-09-14 | R

何度か取り上げた被験者内２要因をRで実行する例のまとめ。要因の各水準の組み合わせの記述統計，単純主効果，多重比較をあわせておこなった。単純主効果は，ある水準（A1）のデータをとりだして，もう一方の要因（B）の効果を検定する方法を採用した。詳しくは以下を参照。
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の効果を同様に繰り返す。

goo blog お知らせ

	ブログを読むだけ。毎月の訪問日数に応じてポイント進呈
	gooブロガーの今日のひとこと
	訪問者数に応じてdポイント最大1,000pt当たる！
	goo blogは20周年を迎えました！