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

被験者内要因の分散分析(repeated-measure analysis)をRで行う場合、実験要因（固定（効果）モデル）と被験者要因（変量（無作為効果）モデル）の指定がややわかりづらい。以下にそのくわしい解説がある。

Notes on the use of R for psychology experiments and questionnaires
Jonathan Baron: Department of Psychology, University of Pennsylvania　
Yuelin Li: Department of Psychiatry & Behavioral Sciences
http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#SECTION00078000000000000000


特にError() の書き方、変量モデル（Random vs. Fixed Effects ）の解説が参考になる。SPSS、SASの出力との比較にもふれている。
試みに、岩原(1965) p.281 の分散分析（混合モデル）の例を処理してみた。Bが被験者（２人）で実験要因Aのすべての水準（３水準）を経験した、とする。セルの大きさは3。
> data.ex=read.csv("Iwahara23-6.csv",header=T)
> data.ex
A B SCORE
1 a1 b1 7
2 a1 b1 10
3 a1 b1 12
4 a2 b1 15
5 a2 b1 14
6 a2 b1 16
7 a3 b1 2
8 a3 b1 3
9 a3 b1 4
10 a1 b2 5
11 a1 b2 6
12 a1 b2 4
13 a2 b2 12
14 a2 b2 10
15 a2 b2 9
16 a3 b2 10
17 a3 b2 9
18 a3 b2 8

> aov.ex=aov(SCORE~A+Error(B/A),data.ex)
> summary(aov.ex)

Error: B
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 1 5.5556 5.5556

Error: B:A
Df Sum Sq Mean Sq F value Pr(>F)
A 2 149.333 74.667 1.3125 0.4324
Residuals 2 113.778 56.889

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 12 25.3333 2.1111
>

交互作用のMSを分母にして主効果の検定をおこなうので、F=1.31となる。

比較のために、被験者間２要因（被験者3×2×3）とした処理（p.274)を以下に。

> aov.ex=aov(SCORE~A*B,data.ex)
> summary(aov.ex)
Df Sum Sq Mean Sq F value Pr(>F)
A 2 149.333 74.667 35.3684 9.309e-06 ***
B 1 5.556 5.556 2.6316 0.1307
A:B 2 113.778 56.889 26.9474 3.647e-05 ***
Residuals 12 25.333 2.111
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


上の74.667/56.889がrepreated-measure の主効果の検定となる。

ANOVA君

フリーの統計ソフトウェア「R」で動作する分散分析関数

一般化線形モデルの交互作用項：R

R の glm() を使って解析した場合：

生態学のデータ解析 - FAQ 一般化線形モデル
http://hosho.ees.hokudai.ac.jp/~kubo/ce/FaqGlm.html#toc3

Rによる心理統計再入門7：分散分析と多重比較

一元配置（３水準）の分散分析と多重比較（Tukey法）の例

Excelなどでcsv形式データファイルを作成。

> data.ex=read.csv("data1.csv",header=T)
> data.ex
A score
1 a1 12
2 a1 11
3 a1 13
4 a1 14
5 a1 11
6 a2 13
7 a2 15
8 a2 16
9 a2 17
10 a2 15
11 a3 21
12 a3 19
13 a3 18
14 a3 19
15 a3 22

> aov.ex=aov(score~A,data.ex)
> summary(aov.ex)
Df Sum Sq Mean Sq F value Pr(>F)
A 2 146.533 73.267 33.303 1.266e-05 ***
Residuals 12 26.400 2.200
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> TukeyHSD(aov.ex,"A" , ordered = TRUE)
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered

Fit: aov(formula = score ~ A, data = data.ex)

$A
diff lwr upr p adj
a2-a1 3.0 0.4973221 5.502678 0.0194356
a3-a1 7.6 5.0973221 10.102678 0.0000092
a3-a2 4.6 2.0973221 7.102678 0.0009792

>

Rによる心理統計再入門6：２つの平均値の差の検定と信頼限界

よく使われる２つの平均値の差の検定と、差の95％の信頼区間をもとめる例。


> x1<-c(3,4,5,6,7)
> x2<-c(0,1,2,3,4)
> t.test(x1,x2)

Welch Two Sample t-test

data: x1 and x2
t = 3, df = 8, p-value = 0.01707
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.6939959 5.3060041
sample estimates:
mean of x mean of y
5 2

同様にデータに対応がある場合：
> x2<-c(1,1,2,3,4)
> x1<-c(3,4,5,6,6)
> t.test(x1,x2,paired=TRUE)

Paired t-test

data: x1 and x2
t = 10.6145, df = 4, p-value = 0.0004460
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
1.919913 3.280087
sample estimates:
mean of the differences
2.6

Rによる心理統計再入門5：母平均の推定

母平均の区間推定および１つの平均値の検定

t.test()を使う。今回はc()でデータを直接読み込んでみる(cf.)。
９つのデータから母平均がmu=6といえるか検定する。また95％の信頼区間をもとめる。

> x <-c(1,2,3,4,5,6,7,8,9)
> x
[1] 1 2 3 4 5 6 7 8 9
> mean(x)
[1] 5
> var(x)
[1] 7.5


> t.test(x,mu=6)

One Sample t-test

data: x
t = -1.0954, df = 8, p-value = 0.3052
alternative hypothesis: true mean is not equal to 6
95 percent confidence interval:
2.894916 7.105084
sample estimates:
mean of x
5

>



cf.
>help(c)
c():The default method combines its arguments to form a vector.

Rによる心理統計（再）入門４：母比率の推定

正答率、政党支持率、視聴率などの母数を推定する例


　たとえばサイコロを６００回ふって「１の目」が９５回出たとき、このサイコロの１の目の出る比率を推定する。
Rではbinom.test()関数が利用できる。この関数は、母比率の関する検定と区間推定を行う。

> binom.test(95,600,1/6)

Exact binomial test

data: 95 and 600
number of successes = 95, number of trials = 600, p-value = 0.6221
alternative hypothesis: true probability of success is not equal to 0.1666667
95 percent confidence interval:
0.1300309 0.1900465
sample estimates:
probability of success
0.1583333

仮定される正反応（仮定される「1の目の出る」）確率を1/6としてみると、その比と観測比(95/600)の検定が行われる。この例では、95/600と1/6は有意に異なるとはいえない（p=0.6221）。
デフォルトで95%の信頼区間がもとめられる。
0.1300309　 ＜　p^　＜0.1900465
したがって、このサイコロは「マトモである」と結論する。

　比較のために、６００００回サイコロをふって９５００回１の目がでた場合で母比率を推定すると、以下のようになる。
> binom.test(9500,60000,1/6)

Exact binomial test

data: 9500 and 60000
number of successes = 9500, number of trials = 60000, p-value = 3.58e-08
alternative hypothesis: true probability of success is not equal to 0.1666667
95 percent confidence interval:
0.1554206 0.1612792
sample estimates:
probability of success
0.1583333

この場合、「このサイコロはイカサマである」となる。


　なお、母比率の推定のわかりやすい説明は、以下などを参照する。試行回数が多ければ、二項分布が正規分布で近似できることを利用する。ただし、問題とする確率が0.5から大きく外れるとき（たとえば0.01など）、二項分布はポアソン分布で近似される。

http://www.kwansei.ac.jp/hs/z90010/sugakuc/toukei/hiritsu/hiritsu.htm

Rによる心理統計入門３：分散分析

基本的な分散分析の例。
２要因完全無作為計画（被験者間）で、要因A、要因B、交互作用ABを検定する。関数aov()を使用する。

たとえば、データファイルが独立変数A、B、従属変数SCOREの３列で入力されているData.csvならば、、
> data.ex=read.delim("Data.csv",header=T)
> aov.ex=aov(SCORE~A*B,data.ex)
結果を参照するには、以下で。
> summary(aov.ex)

Rによる心理統計入門２：χ2乗検定

最も基礎的な検定として、分割表で独立性の検定をおこなう。関数chisq.test()を使用する。

まず、分割表を行列として読み込む。
> ct<- matrix(c(25,35,67,33),ncol=2,byrow=T)
> ct
[,1] [,2]
[1,] 25 35
[2,] 67 33

> chisq.test(ct)

Pearson's Chi-squared test with Yates' continuity correction

data: ct
X-squared = 8.8389, df = 1, p-value = 0.002949

連続性の修正をおこなわないならば、
> chisq.test(ct,correct=F)

Rによる心理統計入門１：データの読み込み

Rをインストールした後は、データファイルの作成と読み込みが最初のステップになる。
　カンマ（,）　Tab　空白　などで区切られたテキストファイルの形式によって以下の関数にいずれかでデータを読み込む。

read.csv
read.delim
read.table

たとえば、データファイルData.csvを読み込むには、
>data.ex=read.csv("Data.csv",header=T)
あるいは、
> data.ex=read.delim("Data.csv",header=T)
> data.ex

References
自由書式テキストファイルからデータを読む（群馬大学：青木研究室）
http://aoki2.si.gunma-u.ac.jp/R/io/free.html

biostat Wiki

医学統計処理(特に共分散構造分析)の解説
http://akimichi.homeunix.net/hiki/biostat/


R言語による医学統計
Mxによる共分散構造分析

SEM with R

Rで共分散構造分析を実行する。

AN INTRODUCTION TO STRUCTURAL EQUATION MODELING
John Fox
(Department of Sociology, McMaster University, Canada)

http://socserv.mcmaster.ca/jfox/Courses/Oxford-2006/index.html
NotesとR Script fileが用意されている。
Structural Equation Models (sem)
http://socserv.mcmaster.ca/jfox/Misc/sem/index.html

彼のサイトにあるように、sem package を CRAN (the Comprehensive R Archive Network)からダウンロードする。
CRAN (the Comprehensive R Archive Network)
↓
CRAN - Contributed Packages
http://cran.r-project.org/src/contrib/PACKAGES.html#sem

CRAN - Package sem
http://cran.r-project.org/src/contrib/Descriptions/sem.html
OSX, Windows版あり。

Rと項目反応理論

cf.
R: Statistical Software for Psychology Research
http://www.personality-project.org/r/#irt
http://www.personality-project.org/r/html/irt.person.rasch.html

正規分布の図示３

有意水準α（両側検定）の棄却域を図示する関数。
無駄に長い。必要がなければしない作業。

normd2 <- function (α)
{
z <- seq(-3,3,0.01)
x <-z
plot(z,dnorm(z,mean=0,sd=1.0),type="n")
curve(dnorm(x,mean=0,sd=1.0),type="l",add=T)
alpha <- α
title("Alpha=0.05")
zmin <- -3
zmax <- 3
critical.left <- qnorm(alpha/2, mean=0, sd=1.0)
xaxis <- seq(zmin, critical.left, length=100)
yaxis <- c(dnorm(xaxis, mean=0, sd=1.0), 0, 0)
yaxis <- c(dnorm(xaxis, mean=0, sd=1.0), 0, 0)
xaxis <- c(xaxis, critical.left, xmin)
polygon(xaxis, yaxis, density=25)
critical.right <- qnorm(alpha/2, mean=0,sd=1.0,lower.tail=F)
xaxis <- seq(critical.right, zmax, length=100)
yaxis <- c(dnorm(xaxis, mean=0, sd=1.0), 0, 0)
xaxis <- c(xaxis, zmax, critical.right)
polygon(xaxis, yaxis, density=25)
ypos <- dnorm(critical.left, mean=0, sd=1.0)
text(zmin, ypos, "rejection?nregion", adj=0)
text(zmax, ypos, "rejection?nregion", adj=1)
text((critical.left+critical.right)/2, 2*ypos+0.02, "acceptance region")
xaxis <- c(rep(critical.left,2), rep(critical.right,2))
yaxis <- c(2*ypos-0.02, 2*ypos, 2*ypos, 2*ypos-0.02)
lines(xaxis,yaxis)

}

> normd2(0.05)

双峰分布

平均の異なる正規分布を２つ重ねた双峰分布の作図。たとえば男と女の身長の分布を合わせたような分布。

bimodal <- function(d)
{
x <- seq(-4, 4, 0.05)
plot(x,dnorm(x, mean=0, sd=0.8), type="n")
curve(dnorm(x, mean=-d, sd=0.8)+dnorm(x, mean=d, sd=0.8), type="l",add=T)
}

> bimodal(2)