R言語で区切り文字による文字列の分割
今回は、区切り文字を指定して、文字列を分割してリスト化する方法を紹介します。
strsplit関数
strsplit関数を使うと任意の文字を区切り文字として、文字列を分割できます。例えば以下のように使います。
input <- "abc def ghi" res <- strsplit(input, " ")
区切り文字は正規表現
区切り文字はスペースやタブ(\t)などを指定しますが、実は、正規表現の扱いになっています。なので、+(1文字以上の繰り返し)のような表現を加えると、スペースやタブスペースが連続して入っている箇所を１回の区切りとして扱えます。
#abcの後にはスペースが一つ、defの後にはスペースが２つ入っています。 input2 <- "abc def ghi" res2 <- strsplit(input2, " ") #これだとabc, def,その後に空要素、そしてhgiが続きます。 res3 <- strsplit(input2, " +") #これだとabc, def, hgiだけになります。スペース２文字を区切り文字扱いするので、空要素ができません。
unlist関数でリスト内容にアクセス
strsplit関数の出力結果は、リストオブジェクトです。このオブジェクトだと、要素数が１で、分割した文字列にアクセスできません。そこでunlist関数を利用して、出力結果を配列にします。
input <- "abc def ghi" res <- strsplit(input, " ") length(res)#1がかえります res2 <- unlist(res) length(res2)#3がかえります。インデックスで[1]から[3]を指定すると各要素にアクセスできます。
研究者とは何をする人か？それに必要な資質は？
今回は、「研究者って何か？」を考えます。一応研究者であるので、職業としての研究者や自分の経歴を紹介する機会があります。しかし、先日のそのような機会では、そもそも研究者が何かをうまく語れなかったんです。なので、反省文も兼ねて、研究者やその資質について思うことをまとめておきます。
研究って何？
あらためて「研究とは？」と考えてみると、意外に定義が難しいのですが、なるべく短い言葉で言うなら「新たな知の創造」だと思います。今まで誰も知らなかった問題の答えを見つける。または、ある現象を今まで誰も考えなかった方法で説明するといった感じです。言葉としては短いですが、これを成し遂げるのはなかなか大変です。なぜなら、それが信頼できる知識であることを周りに認めてもらわなければならないからです。例えば、新しい知識を裏付けるきちんとした証拠が必要です。それが、実験や調査によって得られるものであれば、誰がやっても結果を再現できるだけの必要十分な手続きの記述がいります。また、結果の解析方法やその解釈に論理の飛躍があってはなりません。そういった一連の作業をまとめ上げたものが論文となります。
次に、知の創造をもう細かく考えると２種類別れると思います。１つは「新しい価値の創造」で、新たな知見が世の役に立つような場合です。もう一つは「新たな知識の創造」で、役に立つかどうかはおいておいて、今までわかっていなかったことを知識として定着させるような場合です。ただし、この２つは排他的なものではなく、両者を兼ね備えることも可能だと思っています。これについては、「パスツールの象限」という記事を以前に書きました。
また、「新たな」知の創造となるためには、「今までに何がわかっているか」をきちんと把握することも重要です。よく言われる「巨人の肩に乗る」という言葉が表しているように、歴史を作ってきた大研究者やその人たちが作ってきた知識体系に敬意を払い、さらに自分で新しい何かをそこに加えていくのが研究者なのかなと思います*1。
研究者ってどうやってなるもの？
研究さえしていれば、研究者なので、どうやってという問いは本来成り立たないかもしれないですね。しかし、職業としての研究者を目指すのであれば、基本的には大学学部の卒業に加えて、修士または博士課程の大学院を修了する必要があるでしょう。特に、大学でのポストを望むのであれば、最近は博士号取得がほぼ必須になっています。ただし分野によっては、一度専門職として社会に出てから、研究者として大学・研究所に入るようなキャリアパスも考えられます。
研究者として重要な資質は何？
世の研究者のイメージとしては「勉強が得意・大好きな人」なのかなと思います。確かにそういう人が多いです。でも、勉強が得意・好きなだけでは、「研究者」としてやっていけるかどうかはわからないと思います。私が、研究者になくてはならない資質の一つとして考えるのは、「批判的精神をもっていること」です。ここでいう批判とは、理屈もなく他を否定することでもないですし、揚げ足取りをすることでもないです。上記のような、新たな知となる可能性があるものを「本当にその説明でいいのか」疑ってかかる精神のことを指しています。これは、他に向かうばかりでなく、ときには自分の考え、データ、解析方法も批判的に見る必要があります。例えば、どんなに業績のある人が言ったことでも、教科書に載っていることでも「本当にそれでいいか」自分の頭で一度考えてみるのが、批判的精神です。常に疑ってかかるという意味では、かなり面倒臭い人です。でも、それがあることで、真摯な研究が成り立ちます。
ただ、誤解して欲しくないのは、批判的な精神は、もとからある性格のようなものではなく、トレーニングによって習得するものだと思います。なので、これから研究者を目指すあなたが現時点で批判的精神を持ち合わせていなくても大丈夫。徐々に的確な批判ができるように訓練していきましょう。
*1:google scholar という学術文献検索のサービスのページに「巨人の肩に乗る」という言葉があるのはこの意味でしょう
Dockerを利用してRStudioのRマークダウンファイルを使う
今回は、かなり技術的な内容です。Dockerを利用しながら、RStudioのRマークダウンファイルによる解析環境を構築します。ポイントは、日本語フォントが入っていてもきちんとPDF出力がなされるところです。
Dockerって？
Docker社が提供しているサービスです。 PC内に仮想のコンテナをつくって、コンテナ内でいろいろな情報処理を行います。いわゆる仮想PCと違うのは、必要なサービスだけを高速かつ少ない容量で仮想化しているところです。Docker自体の説明は、詳しいページがいくらでも検索できるのでそちらをごらんください。
使用環境
今回はMac OSXに Docker Toolboxをインストールした環境で実施しています。 Docker Toolboxは更新が頻繁にあるようで、しかも更新によって動かなかった機能が動くようになっているので、できるだけ最新を使ったほうが良さそうです。
- Mac OSX 10.10.5(Yosemite)
- Docker Toolbox 1.10.1a
Dockerfileで環境構築
Dockerfileやベースとなるイメージもすでにきちんと動いているものを利用します。この点もDockerのメリットの１つですね。
RStudioが動くイメージの一つにrocker/rstudioがあります。このイメージでは、ベースのRとRstudioがdebianというLinux OSの中で動きます。
https://hub.docker.com/r/rocker/rstudio/
さらにrocker/rstudioから派生したrocker/hadleyverseというイメージファイルがあります。このイメージのDockerfileを今回利用します。
このDockerfileをdocker-machineの中でビルドすれば、当然ちゃんと動きます。しかし、実際にやってみたところ、日本語フォントがうまく読めないようです。そこでDockerfileを修正しながら、日本語フォントでもきちんと動作することを目指します。
Hadleyverseというイメージとの差分
出来上がったDockerfile
Dockerfileは"Dockerfile"という名前で保存する必要がありそうです。
FROM rocker/rstudio
MAINTAINER "jnobuyuki"
RUN apt-get update
RUN apt-get install -y --no-install-recommends \
ibus-mozc \
manpages-ja
RUN apt-get install -y --no-install-recommends imagemagick \
lmodern \
texlive \
texlive-lang-cjk \
texlive-luatex \
texlive-xetex \
xdvik-ja \
dvipsk-ja \
gv \
texlive-fonts-recommended \
texlive-fonts-extra \
&& apt-get clean \
&& cd /usr/share/texlive/texmf-dist \
&& wget http://download.forest.impress.co.jp/pub/library/i/ipafont/10483/IPAfont00303.zip \
&& unzip IPAfont00303.zip \
&& echo "Map zi4.map" >> /usr/share/texlive/texmf-dist/web2c/updmap.cfg \
&& mktexlsr \
&& updmap-sys
## Install some external dependencies.
RUN apt-get update \
&& apt-get install -y --no-install-recommends -t unstable \
default-jdk \
default-jre \
gdal-bin \
icedtea-netx \
libatlas-base-dev \
libcairo2-dev \
libgsl0-dev \
libgdal-dev \
libgeos-dev \
libgeos-c1v5 \
librdf0-dev \
libssl-dev \
libmysqlclient-dev \
libpq-dev \
libsqlite3-dev \
libv8-dev \
libxcb1-dev \
libxdmcp-dev \
libxml2-dev \
libxslt1-dev \
libxt-dev \
netcdf-bin \
qpdf \
r-cran-rgl \
ssh \
&& R CMD javareconf \
&& apt-get clean \
&& rm -rf /var/lib/apt/lists/ \
&& rm -rf /tmp/downloaded_packages/ /tmp/*.rds
## Install the Hadleyverse packages (and some close friends).
RUN install2.r --error \
broom \
DiagrammeR \
devtools \
dplyr \
ggplot2 \
ggthemes \
haven \
httr \
knitr \
lubridate \
packrat \
pryr \
purrr \
reshape2 \
rmarkdown \
rmdformats \
rticles \
rvest \
readr \
readxl \
testthat \
tibble \
tidyr \
tufte \
shiny \
stringr \
xml2
## Manually install (useful packages from) the SUGGESTS list of the above packages.
## (because --deps TRUE can fail when packages are added/removed from CRAN)
RUN install2.r --error \
-r "https://cran.rstudio.com" \
-r "http://www.bioconductor.org/packages/release/bioc" \
base64enc \
BiocInstaller \
codetools \
covr \
data.table \
downloader \
gridExtra \
gtable \
hexbin \
Hmisc \
htmlwidgets \
jpeg \
Lahman \
lattice \
lintr \
MASS \
PKI \
png \
microbenchmark \
mgcv \
mapproj \
maps \
maptools \
mgcv \
nlme \
nycflights13 \
quantreg \
Rcpp \
rJava \
roxygen2 \
RMySQL \
RPostgreSQL \
RSQLite \
testit \
V8 \
XML \
&& r -e 'source("https://raw.githubusercontent.com/MangoTheCat/remotes/master/install-github.R")$value("mangothecat/remotes")' \
&& r -e 'remotes::install_github("wesm/feather/R")' \
&& rm -rf /tmp/downloaded_packages/ /tmp/*.rds
RUN echo "ja_JP.UTF-8 UTF-8" >> /etc/locale.gen
RUN locale-gen
ENV LANG ja_JP.UTF-8
ENV LANGUAGE ja_JP.UTF-8
ENV LC_ALL ja_JP.UTF-8
Dockerでの手順
1. Docker-machineを起動する
2. cd で Dockerfileがあるディレクトリに移動する
3. 以下のようにDockerfileからイメージをビルドする
docker build -t <レポジトリ名/イメージ名> .
4. 以下のようにイメージを起動する。-p オプションはホスト側とコンテナ側のポート番号のマッピングです。両方8787でオッケーです。
docker run -d -p 8787:8787 <レポジトリ名/イメージ名>
5. webブラウザを起動し、Docker-machineのIP:8787 でRstudioのログインが出てきます。デフォルトはID,PASSともにrstudioです。
以上でRSutdioが起動すると思います。
RStudioでの手順
1. New FileでRmarkdown形式のファイルを作る
2. Global optionのSweave項目でviewerをデフォルトのRstudio viewerからSystem Viewerに変更する *1。
3. preamble(ヘッダのような部分)を以下のようにします。タブスペースの個数がポイントです。
---
title: "Document Title"
author: "autho"
date: "****/**/**"
output:
pdf_document:
latex_engine: xelatex
monofont: IPAPGothic
mainfont: IPAPGothic
---このDockerイメージはIPAフォントをダウンロード・インストールしているので、これを使います。latex_engineはlualatexかxelatexが使えます。
5. knit PDFのボタンを押すか、ショートカット（Command + Shit + k）を実行します。
考察：Dockerを使うメリットとデメリット
手順は以上ですが、 Dockerを利用して解析することのメリットとデメリットを考えておきましょう。
メリット
- 突然母艦 PCが亡くなっても、１時間あれば解析環境が完全復活する。
- 共同研究者に同一解析環境を丸ごと渡せる。
- 授業で学生に開発環境をデータ込みで配布できる。
- 古いバージョンのアプリケーションによる解析環境を安全に保存できる。
*1:RStudio viewerは日本語フォントが文字化けする気がします
学会・研究会における保育サービスについて思うこと
今回は、子育て世代研究者あるあるのような話題で、それについて思うことを書きます。
学会における保育サービス
学会や研究会において、子連れでの参加を促すために保育サービスが提供されることがあります。子育て世代の研究者にとって、このサービスはその学会への参加を大いに助けてくれる貴重なものだと言えます。しかし、実際にサービスを利用できるかは、料金次第であるというのが本音です。
サービスの使用料金はいろいろ
まず、このようなサービスは学会・研究会参加者全員へのサービスではないので、受益者負担が生じるのはしかたないと思います。それにもかかわらず、日本教育心理学会第５７回総会や第７９回日本心理学会大会のように無料でサービスが提供されることもあります*1。
confit.atlas.jp
実際には、いくらか利用料がかかることがほとんどです。高くて利用できなかった例としては、１時間当たり１５００円と言われたことがあります*2。
学会でサービスが提供されない場合
大きな学会では、保育サービスは当たり前になりつつあると思いますが、小さな研究会ではまだ少ないです。それでもどうしても参加したい研究会があります。では、どうするか？私の場合は、所属している企業・大学の福利厚生サービスの利用を考えます。実際に、提携している認証保育所へ一時間３００円で預かってもらったこともあります。差額を企業の福利厚生サービスで負担してくれているわけです。１時間当たり３００円なら半日預けても２０００円もかからないので、負担可能な範囲かなと思います*3。
現状の問題とそれへの提案
学会が提供するサービスや会社の福利厚生サービスは、貴重なのですが、それが受けられないケースは結構あると思います。例えば、大学院生は福利厚生サービスを受けられないと思うので、学会が保育サービスを提供していないといきなり行き詰まります。あとは、福利厚生サービスを利用できたとしても、何をどうやって利用するかの手続きが煩雑な場合がほとんどです。学会発表、宿泊先の確保に加えて、保育サービスの手配とその助成手続きとなるので、手続きのための時間を確保できないことがあります。このような問題が生まれる背景として、大学や企業による保育サービスへの支援が、個人を対象としていることが挙げられます。
例えばこれを、保育サービスを提供しようとしている学会・研究会を対象にしたとします。すると、サービスを利用したいと考えている人たちをより確実にサポートできる気がします。では、どうすればこれを確立できるか？最近の学会では、使用費目を明らかにしたスポンサー企業があるようです。
http://www.anlp.jp/nlp2016/#sponsor
上記の懇親会・茶菓スポンサーと同様の扱いで「保育サービススポンサー」を募るのはどうでしょうか？スポンサー企業によっては、保育サービススポンサーになることの宣伝効果を期待できる場合がある気がします。もちろん文部科学省がこのような保育サービス支援助成金の制度を作ることもありだと思います。
10万PVになりました。ありがとうございます。
今回は、このブログを見てくださった皆様へのお礼です。
100,000PV
100,000PVは、このブログを始めた時、そして始めてしばらくの間には考えもつかなかった数字です。しかも2014年11月に10,000pv達成だったので、1年1〜2ヶ月の間に指数的にpvが伸びています。割と思いついたこと、書きたいこと、書いておいた方が良いと思うことを自由に書いています。いわゆるブロガーを志しているわけでもないですし、アフィリエイトなども一切やっていません。そう考えると、100,000pvはできすぎな数字と言えます。
どんな人が見てくださったのか？
これは10,000pv達成した時にも書きましたが、平日のアクセスと休日のアクセスの落差がかなりあります。休日は平日の半分以下のアクセスです。内容としては、基本的な内容を検索で見つけた人が多いように思います。本ブログのタイトルは「webbeginner」なのでその趣旨にかなりあった結果となっています。嬉しい限りです。
そろそろ「ビギナー」ではないのでは？
確かに、「ビギナー」とは言えないなと思う部分（例えばR言語の使い方）とまだまだ初心者同然の部分（例えばD3,p5などのJavaScriptライブラリ）が混在しています。しかし、R言語の使い方については、自分が初学者ということではなく、自分に関わりのある初学者の皆様がウェブから検索で偶然拾ってくれないかなという淡い期待を込めつつ記事を投稿しています。なので「初心者むけ」というポリシーに広い意味で合致していると思っています。
1,000,000pvは現実的にありえないので
あと一桁増やすためには、どう考えても増やすための努力が必要です。最初に書いたようにプロブロガーを目指してこのブログを運営しているわけではなく、このペースが保たれれば十分だと考えています。今後は、少し込み入った内容についても、本業に支障のでない限りで書いてみたいと思っています。そのような内容は、「初心者向け」ではないので、投稿先をqiitaに変えるかもしれません。というわけで、以下が私のqiitaへのリンクです。こちらもぜひよろしくお願いします。
R言語でサイズの大きいファイルの読み込み
今回は、ファイルの入力に関するヒントを書きます。
ファイルからの入力はread.table関数
R言語では、baseパッケージの中にread.table関数、read.csv関数などがあります。これを利用すれば、簡単にテキストファイルとして保存されているデータを読み込めます。しかし、read.table関数で作成されるdata.frameオブジェクトは容量制限が結構きびしく、２０変数で１００万行あるようなデ−タだと容量ぎりぎりになります。問題は、ファイルを中までしか読んでいなくても警告しか出ないことです。そのままコードが走ってしまうので、きちんと確認しないと実はファイルがうまく読み込めていなかったということになりかねません。
data.tableパッケージの利用
そこで、data.tableパッケージを利用します。data.tableオブジェクトは、data.frameよりも大きなデータを高速に読み込むことができます。
data.tableオブジェクトに入力するにはfread関数を用います。
#data.tableパッケージの起動 library(data.table) #ファイルの読み込み dat <- fread("filename.txt", sep = "\t")
なぜ分散分析には自由度が２つあるのか？
今回は、分散分析という解析方法の紹介をします。
平均値の差の検定です。
分散分析は、３つ以上の条件やグループがあるときの平均値の条件差・グループ間の差を比較するときに使います。確率の考え方を取り入れることで、平均の差が偶然と呼べる程度の差なのか、それとも偶然からはありえないほど大きな差になっているのかを明らかにできます。分散分析という名前（英語だとANalysis Of VArianceで大文字の部分をとってANOVAとも言われます ）から想像すると、分散の大きさの比較なのかなと思うかもしれませんが、実は平均の比較をしています。
平均の差の大きさをどうやって評価するのか。
条件やグループの平均の差は、同一条件、同一グループのばらつきに対する相対的な大きさとして評価します。統計的仮説検定としての分散分析は、条件、グループの平均値の差が０であることを帰無仮説*1にします。つまり、観測された条件間のばらつきは偶然得られたに過ぎず、その大きさは、条件内・グループ内の変動と等しいとなります。これを比で示したのがF値という統計量です。
そしてこのF値という統計量には２つの自由度が伴います。
なぜ2つの自由度があるのか
２つの自由度はそれぞれ「分子の自由度」と「分母の自由度」と呼ばれることがあります。実は、これが上記の式に対応しており、分子の自由度は条件・グループの平均の自由度、分母の自由度はグループ内のデータのばらつきの自由度を示しています。以下では３つの条件を比較する１要因分散分析の場合で考えてみます。
条件・グループの平均の自由度
では、分子の自由度の意味を少し考えてみます。分子の自由度は
です。例えば３つのグループを比較するときの自由度は２になります。自由度というのは、値の取り方の自由度です。グループを考慮に入れない全体の平均値がわかっているとき、３つのうち２つのグループの平均値がわかると残り一つのグループの平均値は自然に決まります。
図式的に考えれば以下のようになります。
各グループの全てがグループの平均値と一致している状態です。データがたくさんあっても、データの値の取り方は限られており、全体平均がわかっているなら、あと２つの値は自由にとることができます。
条件内・グループ内のばらつきの自由度
次は、条件内・グループ内のデータのばらつきです。ここでも自由度の値の取り方は決まっています。グループの平均値がもとまっていれば、グループ内に含まれるデータの最後の１つは、他の値がわかると自然に決まります。
これを全てのグループについて足し合わせるので、以下のようになります。
図式的に考えれば、以下のように、各データのグループ平均からのズレを考えていることになります。
Rを使ってF値を求める
ではR言語でF値を計算してみます。基本的には以下のlm関数またはaov関数を利用すれば簡単に計算可能です。
lm関数/ aov関数の利用
特に１つの被験者間要因分散分析ならlm関数かaov関数かで計算結果が変わらないので、どちらでも好きな方を使ってください。
#サンプルデータ dat <- data.frame(dv = c(3.188957, 3.188957, 3.188957, 3.188957, 3.188957, 6.375809, 6.375809, 6.375809, 6.375809, 6.375809, -9.016713, -9.016713, -9.016713, -9.016713, -9.016713), iv = c("red", "red", "red", "red", "red", "orange", "orange", "orange", "orange", "orange", "blue", "blue", "blue", "blue", "blue")) #lm関数の場合 res_lm <- lm(dv ~ iv, data =dat) anova(res_lm) #aov関数の場合 res_aov <- aov(dv ~ iv, data =dat) summary(res_aov)
条件間・グループ間の平均値の差と条件内・グループ内のデータのばらつきを比較する
次に、条件間・グループ間の平均値の差をグループ内のデータのばらつきを比較するすることを具体的な数値計算とともに考えてみましょう。
仮想のデータセットをつくる
ここでは、本来全く差がない３つのグループを考えます。R言語のrnorm関数では、任意の平均値と分散のデータをランダムに作成できるので、これを利用しましょう。
dat <- data.frame(dv = c(rnorm(5),rnorm(5),rnorm(5)),iv = c("red", "red", "red", "red", "red", "orange", "orange", "orange", "orange", "orange", "blue", "blue", "blue", "blue", "blue"))
このデータにはred,orange,blueの３つのグループがあり、それぞれ平均０、分散１の分布から５つのデータがサンプルとして取られています*2。
条件間・グループ間の平均値のばらつきの分布
条件間・グループ間の平均値が、全体の平均に比べてどの程度ばらつくのかは、平方和という値で考えます。各データが、グループの平均値に一致しているとみなして、それが全体の平均からのズレを計算し、二乗します。それを足し合わせたものが平方和です。二乗は平方とも呼ばれるので、二乗したものの合計という意味で平方和という言い方をしています。
R言語を利用するなら例えば以下のように計算できます。
#データ全体の平均値を計算する mean_dat <- mean(dat$dv) #各グループの平均値を計算する mean_red <- mean(dat[dat[,"iv"]== "red","dv"]) mean_orange <- mean(dat[dat[,"iv"]== "orange","dv"]) mean_blue <- mean(dat[dat[,"iv"]== "blue","dv"]) #平方和を求める MSS_group <- ((mean_red - mean_dat)^2*length(dat[dat[,"iv"]== "red","dv"]) + (mean_orange - mean_dat)^2*length(dat[dat[,"iv"]== "orange","dv"]) + (mean_blue - mean_dat)^2*length(dat[dat[,"iv"]== "blue","dv"]))
これを何度も繰り返すと、以下のように平方和の分布を求めることができます。
#上述の平方和を求める数式を関数にする getSS_between <- function(){ dat <- data.frame(dv = c(rnorm(5),rnorm(5),rnorm(5)),iv = c("red", "red", "red", "red", "red", "orange", "orange", "orange", "orange", "orange", "blue", "blue", "blue", "blue", "blue")) #データ全体の平均値を計算する mean_dat <- mean(dat$dv) #各グループの平均値を計算する mean_red <- mean(dat[dat[,"iv"]== "red","dv"]) mean_orange <- mean(dat[dat[,"iv"]== "orange","dv"]) mean_blue <- mean(dat[dat[,"iv"]== "blue","dv"]) #平方和を求める SS_between <- ((mean_red - mean_dat)^2*length(dat[dat[,"iv"]== "red","dv"]) + (mean_orange - mean_dat)^2*length(dat[dat[,"iv"]== "orange","dv"]) + (mean_blue - mean_dat)^2*length(dat[dat[,"iv"]== "blue","dv"])) return(MSS_group) } #100000回繰り返す res <- replicate(100000,getSS_between())
実際にやってみたところ、このような分布ができました。
理論的には、自由度２のカイ二乗分布となります。
条件内・グループ内のデータのばらつき
では、次に、条件内・グループ内のデータのばらつきを考えてみましょう。先ほどのグループ間の計算と同じように、平均値からのズレを二乗したものを合計して平方和を求めます。
#グループ内の平方和を求める関数 getSS_within <- function(){ dat <- data.frame(dv = c(rnorm(5),rnorm(5),rnorm(5)),iv = c("red", "red", "red", "red", "red", "orange", "orange", "orange", "orange", "orange", "blue", "blue", "blue", "blue", "blue")) #各グループの平均値を計算する mean_red <- mean(dat[dat[,"iv"]== "red","dv"]) mean_orange <- mean(dat[dat[,"iv"]== "orange","dv"]) mean_blue <- mean(dat[dat[,"iv"]== "blue","dv"]) #平方和を求める SS_within <- sum((dat[dat[,"iv"] == "red","dv"] - mean_red)^2) + sum((dat[dat[,"iv"] == "orange","dv"] - mean_orange)^2) + sum((dat[dat[,"iv"] == "blue","dv"] - mean_blue)^2) return(SS_within) } #100000回繰り返す res <- replicate(100000,getSS_within())
実際に１０００００回繰り返して分布を描き出すと以下のようになります。
これは理論として自由度１２（データ数１５から３を引いた）のカイ二乗分布に一致しました。
F分布
F分布は、２つの平方和のそれぞれを自由度で割った値（平均平方和）の比で求めます。分布の形状としては、分子であるグループ間の平均値のばらつきを示す平方和の分布形状に影響を受けやすく、類似します。
まとめ
以上、長くなりましたが、分散分析の２つの自由度は、条件間・グループ間の平均値の差に関する自由度と条件内・グループ内のデータのばらつきについての自由度をあらわてしていることをみてきました。途中で示したように、lm関数やaov関数を利用すればあっという間にF値とp値が求まるのですが、たまにはこうやって計算を分解して見ながら考えると、統計量をより深く理解できますね。
