R Markdownで文章中に統計量を表示させる。

いつの間にかR StudioでR Markdownのノートブック形式が扱えるようになっていたので，久々にR Studioをいじってますが，どうにもMacの日本語入力との相性が悪いですね。コンソールでの日本語入力には問題ないのに，エディタでは入力を確定するまでカーソルが移動しなかったり，文節区切りを示すラインが表示されなかったりするのが残念です。

さて，ノートブック形式を使用するかどうかは別として，Markdown記法とRを融合したR Markdownは，ちょっとしたレポートの作成には非常に便利です。ただ，ざっと調べて見たところ，R（のggplot2）で作成したグラフをそのまま表示させるような使い方はWebでたくさん紹介されているのに対し，分析結果を文章として報告する方法についてはあまり紹介されていない感じがしたので，今更かもしれませんが自分の備忘録も兼ねてここに書いておきます。

## 分析用Rスクリプト

以下のスクリプトは，ごくごく単純な１要因３水準の分散分析を実行するためのものです。分析結果そのものを表示させるわけではないので，ここでは`include = FALSE`として結果を非表示にしています。なお，以下のコードでは，```` ``` ````を表示させるためにわざと最初の行と最後の行を４文字分インデントしてありますが，実際にはそれらは不要です。

```
    ```{r, include = FALSE}

require(lsr) # 効果量を求めるパッケージ
A<-c(5,4,5,5)
B<-c(3,4,5,4)
C<-c(2,1,3,3)

my.data <- data.frame(測定値=c(A,B,C),
                       条件=c(rep('A',length(A)), rep('B',length(B)),rep('C',length(C))))

my.res <- aov(測定値~条件, data=my.data)

    ```
```

## 分析結果の報告

このデータの分析結果を論文に記載する場合，たとえば心理学系の論文では次のような書き方になります。

> A，B，Cの３条件間で平均値に有意な差が見られるかどうかを見るために，各条件の測定値を従属変数とし，測定条件を独立変数とする１要因分散分析を行った。その結果，測定条件の主効果が有意であった ($F_{(2,9)}=10.77,\ p=.004,\ \eta_p^2=0.705$)。

ここで，分析の結果得られた統計量を文中に直接入力したのでは，データに修正があった場合，それらの数値を手作業で修正しなくてはならなくなってしまいます。

そこで，`R`ブロック中でまず必要な値を別の変数に代入し，それらを文章中にインライン表示させるようにします。この時，`aov()`など，要約結果がリストとして返されるものについては`unlist()`でベクトルに変換しておきます。

```
    ```{r, include = FALSE}

out <- unlist(summary(my.res)) # 分散分析の結果
df1 <- out['Df1'] # 主効果の自由度
df2 <- out['Df2'] # 誤差（残差）の自由度
F.value <- out['F value1'] # F値
p.value <- out['Pr(>F)1'] # 有意確率
eta.sq <- etaSquared(my.res)[2] # 偏イータ２乗

    ```
```
 
そのうえで，必要な値を`` `r ` ``を使って文中に挿入します。ここで注意が必要なのは，$F$値や$p$値はそのままでは小数点以下６桁程度まで表示されてしまいますので，`sprintf()`などを使用して表示桁数を揃える必要があります。また，ここでは統計量を数式モードで表示するようにしてあります。

```
A，B，Cの３条件で平均値に有意な差が見られるかどうかを見るために，
各条件の測定値を従属変数とし，測定条件を独立変数とする１要因分散分析を行ったところ，
測定条件の主効果が有意であった ($F_{(`r df1`,`r df2`)}=`r sprintf('%.3f', F.value)`,
\ p=`r sprintf('%.3f',p.value)`,\ \eta_p^2=`r sprintf('%3.3f',eta.sq)`$)。
```

これをコンパイルすると次のようになります。

> A，B，Cの３条件間で平均値に有意な差が見られるかどうかを見るために，各条件の測定値を従属変数とし，測定条件を独立変数とする１要因分散分析を行った。その結果，測定条件の主効果が有意であった ($F_{(2,9)}=10.77,\ p=0.004,\ \eta_p^2=0.705$)。

このままでも問題ないかもしれませんが，$p$値については`0.004`ではなく`.004`と記載せよという決まりがあるかもしれません。その場合，`substring()`を使って小数点の前の`0`を削除します。

```
A，B，Cの３条件で平均値に有意な差が見られるかどうかを見るために，
各条件の測定値を従属変数とし，測定条件を独立変数とする１要因分散分析を行ったところ，
測定条件の主効果が有意であった ($F_{(`r df1`,`r df2`)}=`r sprintf('%.3f', F.value)`,
\ p=`r substring(sprintf('%.3f',p.value),2)`,\ \eta_p^2=`r sprintf('%3.3f',eta.sq)`$)。
```

これをコンパイルすると次のようになります。

> A，B，Cの３条件間で平均値に有意な差が見られるかどうかを見るために，各条件の測定値を従属変数とし，測定条件を独立変数とする１要因分散分析を行った。その結果，測定条件の主効果が有意であった ($F_{(2,9)}=10.77,\ p=.004,\ \eta_p^2=0.705$)。

### まとめ
データの分析から報告までを１つのファイルで行い，その結果をMS Word形式で書き出すR Markdownファイルは次のようになります。

```
---
title: "R Markdownのサンプル"
output: word_document
---

# 結果

    ```{r, include = FALSE}
require(lsr) # 効果量を求めるパッケージ
A<-c(5,4,5,5)
B<-c(3,4,5,4)
C<-c(2,1,3,3)

my.data <- data.frame(測定値=c(A,B,C),
                       条件=c(rep('A',length(A)), rep('B',length(B)),rep('C',length(C))))

my.res <- aov(測定値~条件, data=my.data)

out <- unlist(summary(my.res)) # 分散分析の結果
df1 <- out['Df1'] # 主効果の自由度
df2 <- out['Df2'] # 誤差（残差）の自由度
F.value <- out['F value1'] # F値
p.value <- out['Pr(>F)1'] # 有意確率
eta.sq <- etaSquared(my.res)[2] # 偏イータ２乗
    ``` 

A，B，Cの３条件で平均値に有意な差が見られるかどうかを見るために，
各条件の測定値を従属変数とし，測定条件を独立変数とする１要因分散分析を行ったところ，
測定条件の主効果が有意であった ($F_{(`r df1`,`r df2`)}=`r sprintf('%.3f', F.value)`,
\ p=`r substring(sprintf('%.3f',p.value),2)`,\ \eta_p^2=`r sprintf('%3.3f',eta.sq)`$)。
```