R入門 ― インストールからデータ読み込み・要約統計・t検定・グラフ作成までの最初の一歩をやさしく解説 ―

本記事にはプロモーション(広告・アフィリエイトリンク)が含まれます。
記事の目次
Toggleこの記事でわかること
- RとRStudioの違いと入れる順番:RStudioはRの代わりではなく、Rを快適に動かすための作業環境。R本体を先に入れること
- 最初に覚える操作はこれだけ:
<-で入れて、c()でまとめて、mean()で計算する。この3つが分かればデータは動かせます - データの形を確認する作法:
head()・str()・summary()の3点セットと、summary()が出した数字をそのまま報告してはいけない理由 - 最初の統計解析:
t.test(len ~ supp, data = tg)の1行で出る p値0.06063 と95%信頼区間[-0.17, 7.57]を、どこを見てどう読むか - グラフとエラー:ヒストグラム・箱ひげ図・散布図を3行で描く方法と、初心者が必ず出会う4つのエラーメッセージの読み方
はじめに
統計を学び始めた方が最初にぶつかる壁は、統計の理論そのものよりも「で、それをどうやって計算するのか」という一点であることが少なくありません。教科書でt検定の式を追えるようになっても、手元のデータでそれを実行できなければ、統計は自分の道具になりません。そして製薬企業やCROで解析に関わる仕事に就けば、Excelの関数だけで済む場面はほとんどなくなります。
R は、その計算を引き受けてくれる無償のソフトウェアです。統計解析のために作られた言語なので、平均を出すのも、t検定をするのも、グラフを描くのも、たいてい1行で終わります。製薬業界では長らくSASが標準の位置を占めてきましたが、近年は解析の探索的な部分やシミュレーション、可視化でRを併用する現場が増え、規制当局への提出パッケージにRを用いた事例も出てきました。学生のうちに触れておく価値も、実務に入ってから学び直す価値も、どちらも十分にあります。
本記事は、Rを一度も起動したことがない方を対象に、インストールから最初の統計解析までを1本で通します。扱うのは、RとRStudioの違いとインストールの流れ、ベクトルとデータフレームという2つの基本形、データの読み込みと要約統計量、2群を比べるt検定、そしてヒストグラム・箱ひげ図・散布図です。最後に、初心者が必ず出会う4つのエラーメッセージの読み方を付けました。
題材には、Rに最初から入っている ToothGrowth という練習用データを使います。モルモット60匹に2種類の形態でビタミンCを与え、投与量を3水準に変えて歯の伸びを測った実験データです。臨床試験そのものではありませんが、「2つの群を比べる」「用量が上がると反応も上がるか見る」という問いの形は、臨床試験で毎日扱っているものと同じです。ダウンロードの手間もなく、あなたの手元でもまったく同じ数値が再現できます。
掲載しているコードはすべてそのまま実行でき、出力は R 4.5.1 で実際に得られた値です。上から順にコピーして貼り付けていけば、記事を読み終えたときには最初の解析が1つ終わっている、という構成にしました。統計手法そのものの理論には深入りせず、必要なところで詳しい既存記事へのリンクを置いています。まずは動かすこと、それがこの記事の目的です。
Rとは何か ― RStudioとの違いとインストールの流れ
統計解析を始めようと調べると、必ず「R」と「RStudio」という2つの名前が出てきます。名前は似ていますが役割はまったく違い、ここを取り違えたままだとインストールの最初でつまずきます。2つの関係と、何をどの順番で入れればよいのかを整理しておきましょう。
Rは統計のために作られた無償のソフトウェア
Rは、統計解析とグラフ作成のために設計されたプログラミング言語であり、それを実行するソフトウェアでもあります。配布元は CRAN(The Comprehensive R Archive Network、https://cran.r-project.org/)で、個人でも企業でも無償。ライセンス費用も使用許諾の申請も要りません。
Rの性格は3つに整理できます。1つ目は、統計に特化していることです。平均や標準偏差、t検定、回帰分析といった計算が最初から関数として用意されています。関数とは、名前を呼び出すと決まった処理をしてくれる小さな道具のこと。平均を出したければ mean を呼ぶというように、教科書の用語がそのまま命令の名前になっています。
2つ目は、パッケージで後からいくらでも広げられることです。パッケージとは、誰かが作った関数やデータをひとまとめにした追加部品のこと。CRANには膨大な数が公開されていて、生存時間解析、混合効果モデル、見栄えのよい作図まで、必要になったときに足していけます。新しい手法は、著者自身がパッケージを公開することも珍しくありません。
3つ目は、操作の手順が文字として残ることです。マウス操作は記録に残りませんが、Rでは「何をしたか」がそのままコードとして残り、半年後に同じファイルを実行すれば同じ数字が出ます。この再現性は、研究や申請業務では決定的な価値を持ちます。
位置づけを、現場でよく使われるSASやExcelと並べて整理しておきます。
| 観点 | R | SAS | Excel |
|---|---|---|---|
| 費用 | 無償 | 有償のライセンス契約 | Officeに付属(有償) |
| 主な使われ方 | 学術研究やデータ解析全般。製薬の現場でも広く使われる | 承認申請の解析業務で長く使われてきた | 集計、表づくり、簡単なグラフ |
| 手順の再現 | スクリプトに残り、何度でも再現できる | 同じくプログラムとして残る | マウス操作は残らず、後から追跡しにくい |
| 新しい手法への対応 | パッケージを入れればすぐ試せる | 公式に実装されるのを待つことが多い | 標準機能では限界がある |
| 学び始めやすさ | すぐ入手できるが打鍵に慣れが必要 | 個人では入手しにくい | 最も直感的 |
どのソフトを選ぶかは目的と職場の事情で決まります。選び方そのものは生物統計家に必要なスキル(統計解析ソフトウェア)で整理しています。
RStudioはRを動かすための作業環境
RStudio は、Posit Software, PBC が提供する統合開発環境です。統合開発環境とは、コードを書く、実行する、結果や図を見るといった作業を1つの画面にまとめた道具のこと。無償のオープンソース版が公開されていて、Windows、macOS、Linux で動きます。
押さえてほしいのは、RStudioはRの代わりではないということです。車にたとえれば、計算をしているエンジンがR、ハンドルやメーターにあたるのがRStudio。RStudioだけを入れても、Rが入っていなければ何も計算できません。
したがって、インストールは次の順番になります。
まずR本体です。CRANのトップページには「Download R for Windows」「Download R for macOS」「Download R for Linux」という案内が並んでいます。自分のパソコンに合うものを選び、Windowsならさらに「base」に進んでインストーラをダウンロードします。設定は既定のままで問題ありません。
次にRStudioです。Posit の公式サイト(https://posit.co/products/open-source/rstudio/)のダウンロード案内から、自分のOS用のファイルを選びます。macOSでは搭載しているCPUに応じて選ぶ場合があります。ページの文言やボタンの位置は更新されることがあるため、細かい表記は公式サイトの案内に従ってください。大事なのは「Rが先、RStudioが後」という順番だけです。
インストールが終わったら、以後起動するのはRStudioです。R本体のアイコンは基本的に触りません。
起動直後の画面と、コンソール・スクリプトの違い
RStudioを初めて起動すると、画面がいくつかの区画に分かれています。呼び名と役割を先に知っておくと迷子になりません。
- コンソール:打ち込んだ命令をその場で実行し、結果を返す場所
- スクリプト(ソース):命令を書きためておく手順書。初回起動時は表示されていないことがありますが、メニューから新しいRスクリプトを作れば現れます
- 環境(Environment):いま手元にあるデータや計算結果の一覧
- プロット・ファイル・ヘルプ:作った図の表示、フォルダの中身、関数の説明書
初心者が最初に迷うのが、コンソールに直接打つのとスクリプトに書いて実行するのは何が違うのか、という点です。結果そのものは同じで、違うのは記録が残るかどうかだけです。
コンソールは電卓で、打ち込めばすぐ答えが返りますが、閉じてしまえば手順は消えます。スクリプトは手順書で、ファイルとして保存でき、明日でも来月でも同じ順番で実行し直せます。行にカーソルを置いて Ctrl と Enter(macOSでは Command と Enter)を押せば、その行だけがコンソールに送られて実行されます。ちょっと試すだけならコンソール、後で見返す可能性があるならスクリプト。実務のデータ解析では、例外なくスクリプトに書きます。手順が残らない解析は、検証のしようがないからです。
もうひとつ、初日に必ず出会うのがパッケージの扱いです。使うには2段階の手順が必要で、この2つを混同するとエラーの原因になります。
1つ目が install.packages("ggplot2") です。インターネットからパッケージをダウンロードして保存します。丸かっこの中に書いた "ggplot2" を引数と呼びます。引数とは、関数に渡す材料のこと。この作業はパソコン1台につき原則1回だけで済みます。
2つ目が library(ggplot2) です。保存されたパッケージを、いま使えるように読み込みます。こちらはRを起動するたびに毎回必要です。本にたとえるなら、install.packages() は本を買って本棚に入れる作業、library() は本棚から机に出す作業です。なお、R本体だけで install.packages() を実行すると、ダウンロード元のミラーサイトを選ぶ画面が出ることがあります。近い地域のものを選べば問題ありません。
インストールは1回、読み込みは毎回。この区別を覚えておくと、エラーの切り分けが一気に楽になります。「そんなパッケージはない」というエラーはインストールがまだ済んでいない合図、「その関数が見つからない」は
library() を忘れている合図です。スクリプトの先頭に必要な library() をまとめて書くのが定番です。学習の進め方全体を見渡したい方は、生物統計家を目指すための勉強法もあわせて読んでみてください。
Rの基本操作 ― ベクトルとデータフレームを触ってみる
ここからは実際に手を動かします。RStudioを起動し、コンソールに1行ずつ打ち込んでください。コピー&ペーストでも構いませんが、最初は自分で打ったほうが指と目に残ります。
まずRが正しく入っているかを確認します。
R.version.string
[1] "R version 4.5.1 (2025-06-13 ucrt)"
このように文字列が返ってくれば、Rは正常に動いています。本記事は R 4.5.1 で動作を確認しています。表示される番号は導入した時期によって違いますが、以降のコードは基本的な機能しか使っていないので、番号が違っても同じ結果が得られます。
電卓として使い、数値のかたまりに名前を付ける
まずはRを電卓として使い、そのうえでRらしい「複数の数値をひとまとめにして名前を付ける」に進みます。
1 + 2
sqrt(16)
x <- c(120, 135, 128, 142, 119)
x
mean(x)
sd(x)
length(x)
[1] 3
[1] 4
[1] 120 135 128 142 119
[1] 128.8
[1] 9.833616
[1] 5
上から順に見ていきます。1 + 2 は [1] 3 を返しました。sqrt(16) は平方根を求める関数で、答えは4です。出力の先頭の [1] は「この行は1番目の要素から表示しています」という位置の目印であって、値ではありません。要素が多くて折り返すときは、2行目の先頭に [6] のような番号が付きます。
続く x <- c(120, 135, 128, 142, 119) が、この記事でいちばん大事な1行です。c() は combine(つなげる)の頭文字で、複数の値を1本のベクトルにまとめる関数。ベクトルとは、同じ種類の値が順番に並んだ1列のデータのことです。ここでは5人分の収縮期血圧(mmHg)に見立てた5つの数値を並べました。
矢印の記号 <- は代入を表します。右側で作ったものに左側の名前を付けて保存する、という意味です。こうして保存されたものをオブジェクト、つまりデータや計算結果に付けた名札付きの箱と呼びます。以後は x と打つだけで、5つの数値をまとめて呼び出せます。RStudioの環境ペインにも x が現れているはずです。
そして mean(x) が平均の128.8、sd(x) が標準偏差の9.833616、length(x) が要素の個数5を返しました。数値をそのつど書き並べる必要はなく、x という名札を関数に渡すだけで済みます。
中身を要約し、必要な部分だけ取り出す
ベクトルを作ったら、全体を眺めたり一部だけを抜き出したりします。ここがRの操作感の中心です。
summary(x)
x[2]
x[x > 130]
x / 10
Min. 1st Qu. Median Mean 3rd Qu. Max.
119.0 120.0 128.0 128.8 135.0 142.0
[1] 135
[1] 135 142
[1] 12.0 13.5 12.8 14.2 11.9
summary(x) は、最小値、第1四分位数、中央値、平均、第3四分位数、最大値をまとめて返します。数値を1つずつ眺めなくても、分布のおおよその形をひと目で確認できます。
x[2] は角かっこで「2番目の要素」を指定した書き方で、135が返りました。Rの番号は1から始まります。
x[x > 130] はもう一段おもしろい書き方です。角かっこの中に条件を書くと、その条件を満たす要素だけを取り出せます。ここでは130より大きい135と142が返りました。データの絞り込みは、解析で最も頻繁に行う操作です。
そして x / 10 です。5つの要素すべてに一度に割り算が適用され、答えが5つ返っています。多くの言語なら「1個ずつ順番に処理する」繰り返しの命令を書く場面ですが、Rでは要りません。1つの数値ではなくデータのかたまりを一度に扱えること、これがRという道具の本質です。
表の形のデータを扱う ― データフレーム
実際のデータは1列では足りません。被験者番号があり、割り付けられた群があり、測定値がある、という表の形をしています。Rでこの形を受け持つのがデータフレームで、行が1人分(1件分)、列が項目に対応するExcelの表と同じ構造です。
df <- data.frame(
subjid = c("S01", "S02", "S03", "S04", "S05", "S06"),
arm = c("Active", "Active", "Active", "Placebo", "Placebo", "Placebo"),
sbp = c(128, 135, 121, 142, 139, 147)
)
df
str(df)
dim(df)
df$sbp
df[df$arm == "Active", ]
mean(df$sbp[df$arm == "Active"])
subjid arm sbp
1 S01 Active 128
2 S02 Active 135
3 S03 Active 121
4 S04 Placebo 142
5 S05 Placebo 139
6 S06 Placebo 147
'data.frame': 6 obs. of 3 variables:
$ subjid: chr "S01" "S02" "S03" "S04" ...
$ arm : chr "Active" "Active" "Active" "Placebo" ...
$ sbp : num 128 135 121 142 139 147
[1] 6 3
[1] 128 135 121 142 139 147
subjid arm sbp
1 S01 Active 128
2 S02 Active 135
3 S03 Active 121
[1] 128
data.frame() に3本のベクトルを渡して、6人分の小さなデータを作りました。被験者番号、割り付け群、収縮期血圧の3列です。列の長さは必ずそろえてください。1つでも違うとエラーになります。df と打つと表が表示され、左端に付く1から6の番号は行番号であって、データの中身ではありません。
str() は structure(構造)の略で、中身の設計図を教えてくれる関数です。6 obs. of 3 variables は6行3列という意味。続く行の chr は文字列(character)、num は数値(numeric)を表します。この型の違いは初心者がつまずく代表格で、数値のつもりの列が文字列になっていると、平均を計算した瞬間にエラーになります。データを読み込んだら必ず str() で型を確かめる習慣を付けてください。dim() は行数と列数だけを返します。
df$sbp のドル記号は「このデータフレームの、この列」を指します。1列だけ取り出すとベクトルに戻るので、mean() や summary() がそのまま使えます。
最後の2行が山場です。df[df$arm == "Active", ] の角かっこの中にはカンマがあり、前が行、後ろが列を指定する場所です。ここでは前に「armがActiveと等しい行」という条件を書き、後ろは空欄=「すべての列」にしています。結果として、Active群の3行が全列そろって返りました。等号を2つ重ねた == は「等しいかどうかを判定する」記号で、1つの = は代入になるため区別が必要です。
mean(df$sbp[df$arm == "Active"]) は、Active群だけの平均で128です。こちらは df$sbp というベクトルに条件を書いているので、カンマは要りません。データフレームに角かっこを使うときだけカンマが必要、と整理しておくと混乱しません。
角かっこのカンマは省略しないでください。
df[1, ] は1行目の全列、df[, 1] は1列目の全行、df[1, 1] は1行1列目の値です。ところがカンマを忘れて df[1] と書くと「1列目」と解釈され、エラーにならないまま意図と違う結果が返ります。黙って通るぶん発見が遅れます。行を選ぶときはカンマを必ず打つ、と覚えてください。ここまでに出てきた書き方を一覧にしておきます。覚えようとせず、必要になったら見返す早見表として使ってください。
| 書き方 | 何をするか | 例 |
|---|---|---|
<- | 右側の結果に名前を付けて保存する | x <- c(120, 135) |
c() | 複数の値を1本のベクトルにまとめる | c(120, 135, 128) |
mean() / sd() | 平均と標準偏差を計算する | mean(x) |
length() | 要素がいくつあるかを数える | length(x) |
summary() | 要約統計量をまとめて出す | summary(x) |
| 角かっこ | 位置や条件で要素を取り出す | x[2]、x[x > 130] |
data.frame() | 複数の列を表の形にまとめる | data.frame(sbp = x) |
str() / dim() | 列の型と、行数・列数を確認する | str(df) |
$ | データフレームから1列を取り出す | df$sbp |
これで、値をまとめる、要約する、条件で取り出す、表として扱うという4つの基本動作が身につきました。ここまでの操作は、扱うデータが6行でも6万行でも同じです。手を動かしながら統計の考え方も固めたい方は、生物統計学を学ぶためのおすすめ書籍で入門書を選んでみてください。
データを読み込み、要約統計量で全体像をつかむ
ここからは、もう少し「データらしい大きさ」のものを扱います。とはいえ、いきなり手元のファイルを読み込もうとすると、ファイルの場所や文字コードでつまずいて肝心の解析にたどり着けません。まずはRに最初から入っている練習用のデータを使います。誰の環境でも同じ結果が出るので、自分の書き方が正しいのかを確かめながら進められます。
Rに最初から入っている練習用データを開く
この記事では ToothGrowth というデータを使います。モルモット60匹に、ビタミンCを2種類の投与形態(オレンジジュース=OJ、アスコルビン酸=VC)で、3水準の用量(0.5・1・2 mg/day)に分けて与え、歯の長さを測った小さな実験データです。1つの組み合わせにつき10匹ずつ、合計60行。「2つの与え方を比べる」「量を増やすと反応も大きくなるか見る」という問いの形が臨床試験とまったく同じなので、練習台としてはうってつけです。
新しいデータを受け取ったら、いきなり検定を走らせてはいけません。まず「何行あるか」「どんな列があるか」「各列はどんな型か」を確認します。この確認には決まった4つの関数を使います。
head(ToothGrowth)
str(ToothGrowth)
dim(ToothGrowth)
summary(ToothGrowth)
len supp dose
1 4.2 VC 0.5
2 11.5 VC 0.5
3 7.3 VC 0.5
4 5.8 VC 0.5
5 6.4 VC 0.5
6 10.0 VC 0.5
'data.frame': 60 obs. of 3 variables:
$ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
$ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
$ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
[1] 60 3
len supp dose
Min. : 4.20 OJ:30 Min. :0.500
1st Qu.:13.07 VC:30 1st Qu.:0.500
Median :19.25 Median :1.000
Mean :18.81 Mean :1.167
3rd Qu.:25.27 3rd Qu.:2.000
Max. :33.90 Max. :2.000
4つの出力を順に見ていきます。head() は先頭6行だけを表示する関数で、全部を画面に流さずに「だいたいこんな形か」をつかめます。列は len(歯の長さ)、supp(投与形態)、dose(用量)の3つです。
str() は前の節でも使った、各列の型を教えてくれる関数です。num は数値、Factor w/ 2 levels "OJ","VC" は「OJとVCの2種類しかない分類の列」という意味。数値の列と分類の列を区別しておくことは、あとで検定やグラフを書くときに効いてきます。dim() の [1] 60 3 は60行3列を表します。
summary() は各列の要約です。数値の列には最小値から最大値までの6つの値が並びます。len は最小4.20、最大33.90、平均18.81。8倍近い開きがあり、かなり散らばったデータだとわかります。分類の列である supp については、OJが30匹、VCが30匹と個数が表示されます。
ここで、初心者が最も引っかかりやすい落とし穴に触れておきます。
出力の右端に
dose の平均 1.167 と出ていますが、この数字に意味はありません。用量は 0.5・1・2 という3つの水準しかなく、それを20匹ずつ平均しただけの値だからです。「1.167 mg を投与した群」はどこにも存在しません。summary() が出した数字を、全部そのまま報告してはいけないということです。どの列が本当は分類で平均に意味がないのかは、人間が判断するしかありません。CSVに書き出して、読み込んでみる
練習用データはRの中にありますが、実際の仕事で扱うのは外部のファイルです。ここでは一度CSVとして書き出し、それを読み直します。往復させれば、書き出しと読み込みの両方を一度に体験できます。
write.csv(ToothGrowth, "toothgrowth.csv", row.names = FALSE)
tg <- read.csv("toothgrowth.csv")
head(tg, 3)
nrow(tg)
len supp dose
1 4.2 VC 0.5
2 11.5 VC 0.5
3 7.3 VC 0.5
[1] 60
write.csv() の row.names = FALSE は「行の名前をファイルに書かない」という指定です。これを付け忘れると、1・2・3…という連番だけの列が余分に増えてしまい、Excelで開いたときに列がひとつずれたように見えます。付ける癖をつけてください。
読み込みは read.csv("ファイル名.csv") の1行だけで、自分のデータを読むときもファイル名を差し替えるだけです。読み込んだ結果は tg という短い名前に入れました。以降はこの tg を使っていきます。head(tg, 3) のように2つめに数字を渡すと表示する行数を変えられ、nrow() は行数だけを返します。60行が正しく読めたことが確認できました。
読み込みで失敗する原因のほとんどは、統計の話ではなく「ファイルの置き場所」です。Rには作業ディレクトリ(いま作業している基準のフォルダ)という考え方があり、ファイル名だけを書いた場合、Rはそのフォルダの中を探します。
getwd()
getwd() を実行すると、いまの作業ディレクトリがフルパスで表示されます。ここに表示されたフォルダに、読みたいCSVを置いてください。RStudioなら「Session → Set Working Directory」から変更できます。プロジェクト機能を使えば、プロジェクトのフォルダが自動で作業ディレクトリになり、この手間自体がなくなります。
平均・度数・グループ別平均で全体像をつかむ
データが手元に入ったら、次は要約です。まずは全体の平均と標準偏差から。標準偏差は、データが平均のまわりにどれくらい散らばっているかを表す指標です。
mean(tg$len)
sd(tg$len)
round(mean(tg$len), 2)
table(tg$supp)
table(tg$supp, tg$dose)
aggregate(len ~ supp, data = tg, FUN = mean)
aggregate(len ~ supp + dose, data = tg, FUN = mean)
[1] 18.81333
[1] 7.649315
[1] 18.81
OJ VC
30 30
0.5 1 2
OJ 10 10 10
VC 10 10 10
supp len
1 OJ 20.66333
2 VC 16.96333
supp dose len
1 OJ 0.5 13.23
2 VC 0.5 7.98
3 OJ 1.0 22.70
4 VC 1.0 16.77
5 OJ 2.0 26.06
6 VC 2.0 26.14
上から順に読みます。全体の平均は 18.81333、標準偏差は 7.649315 でした。Rは桁数を丸めずに出しますが、報告に使うなら round(mean(tg$len), 2) のように桁を指定します。2つめの引数 2 が小数点以下の桁数で、結果は 18.81 です。なお、標準偏差とよく混同される標準誤差は別物なので、違いが曖昧な方は【やさしく解説】標準偏差(SD)と標準誤差(SE)の違いを先に読んでおくと、このあとの信頼区間の話がすっきり入ります。
table() は、分類の列が何件あるかを数える関数です。table(tg$supp) はOJが30件、VCが30件。カンマでつないで2つ渡すと、縦と横の組み合わせで数えるクロス集計になります。すべての枠が10で、きれいに10匹ずつ割り当てられた実験だとわかります。この「箱の中身が想定どおりか」の確認は、解析前に必ずやってください。
そして aggregate() が、グループ別に集計する関数です。len ~ supp は「len を supp で分けて見る」という指定で、FUN = mean は「各グループで平均を計算せよ」という意味です。結果は OJ が 20.66333、VC が 16.96333。差はおよそ3.7で、オレンジジュースのほうが長い、と読めます。
ところが、supp + dose のように2つ並べて用量ごとに分けると、話が変わります。
| 用量 | OJの平均 | VCの平均 | 差(OJ-VC) |
|---|---|---|---|
| 0.5 mg/day | 13.23 | 7.98 | 5.25 |
| 1.0 mg/day | 22.70 | 16.77 | 5.93 |
| 2.0 mg/day | 26.06 | 26.14 | -0.08 |
| 全体(用量をまとめた場合) | 20.66 | 16.96 | 3.70 |
0.5 mg と 1 mg では、たしかにOJのほうが5以上高い。ところが 2 mg では 26.06 と 26.14 で、ほとんど同じです。全体をひとまとめにして出した「3.7の差」は、低い用量での差が平均に紛れ込んだ結果であって、どの用量でも同じように3.7の差があるわけではありません。
この観察が、この記事でいちばん覚えて帰ってほしいことです。平均をひとつ出して終わりにせず、意味のある単位で分けて見る。手順としては aggregate() に列をもう1つ足すだけですが、結論はまったく変わります。
臨床試験の解析でも同じことが起こります。全体では効果があるように見えても、年齢層や重症度で分けると特定の層でしか差がない、という場面は珍しくありません。だからこそ、どの単位で分けて見るかは解析計画の段階で決めておきます。
全体平均だけを見ると「OJのほうが効く」ですが、用量別に見ると「差が出るのは低用量のときだけで、2 mgでは追いつく」と読めます。同じデータから出た、まったく違う結論です。集計の単位を変えると結論が変わりうることを知っているかどうかが、解析者の分かれ目になります。
最後に、実データでは避けて通れない欠測値の話をしておきます。欠測値とは「測れなかった・記録がない」ことを表す印で、Rでは NA と書きます。Rに入っているもう一つの練習用データ airquality(ニューヨークの大気環境データ)で確かめます。
sum(is.na(airquality$Ozone))
mean(airquality$Ozone)
mean(airquality$Ozone, na.rm = TRUE)
[1] 37
[1] NA
[1] 42.12931
is.na() は各要素が欠測かどうかを判定し、sum() でその個数を数えられます。37件ありました。そのまま mean() を取ると答えは NA です。エラーではなく、「欠測が混じっているので平均が何かは決められない」とRが正直に答えているのです。黙って37件を捨てて平均を返さないところが、Rの親切さです。
計算したいときは na.rm = TRUE を付けます。na.rm は「NAを取り除く(remove)」の意味で、これで 42.12931 が得られます。ただし、これは37件を無条件に捨てた値です。捨てた37日が特定の季節や、測定装置が止まった日に偏っていないかを考えずに使ってはいけません。おまじないのように付ける癖だけは避けてください。
最初の統計解析 ― t検定で2群を比べる
ここまでで、OJのほうが平均3.7ほど高いことがわかりました。しかし、この3.7は本物の差でしょうか。それとも、たまたま選んだ60匹のばらつきで生じただけでしょうか。この問いに答える最も基本的な道具が、2群の平均を比べるt検定です。Rでは1行で実行できます。
t.test() を1行書いてみる
まずは走らせてみましょう。
t.test(len ~ supp, data = tg)
Welch Two Sample t-test
data: len by supp
t = 1.9153, df = 55.309, p-value = 0.06063
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
-0.1710156 7.5710156
sample estimates:
mean in group OJ mean in group VC
20.66333 16.96333
コードの中身を分解します。len ~ supp の ~(チルダ)は、Rで「左を右で説明する」関係を書く記号で、ここでは「歯の長さを投与形態で分けて比べる」と読みます。data = tg は「その列は tg の中にあります」という指定で、これがあるおかげで tg$len のように毎回データ名を書かずに済みます。
この ~ は、t検定だけの特別な文法ではありません。分散分析でも回帰でもRのモデルはほぼ同じ形で書けるので、ここで覚えておくとこの先ずっと使い回せます。
なお、2群をベクトルで別々に渡す t.test(a, b) という書き方もありますが、実務のデータは1行1個体の表で届くので、表と列名をそのまま渡せる ~ に慣れておくほうが手戻りがありません。
出力は英語で7行ほど並びますが、構造は単純です。上から順に5つのブロックに分かれています。
| 出力の行 | 何が書いてあるか |
|---|---|
| Welch Two Sample t-test | 実行された検定の名前。Rが選んだ方法がここに出る |
| data: len by supp | 何を、何で分けて比べたか。意図どおりか必ず確認する |
| t = 1.9153, df = 55.309, p-value = 0.06063 | 検定統計量・自由度・p値 |
| alternative hypothesis: … not equal to 0 | 対立仮説。2群の平均の差が0ではない、という両側の検定 |
| 95 percent confidence interval | 2群の平均の差の95%信頼区間 |
| sample estimates | 各群の平均 20.66333 / 16.96333。aggregate() の結果と一致 |
3行目の t は、2群の平均の差を、そのばらつきの大きさで割った値です。式で書くと次の形になります。
\[ t=\frac{\bar{x}_1-\bar{x}_2}{SE} \]
分子は2群の平均の差、分母はその差がどれくらいふらつくかを表す標準誤差です。差が大きいほど、あるいはふらつきが小さいほど、t は大きくなります。ここでは 1.9153 でした。
そしてp値は 0.06063。慣習的な5%の基準は下回っていません。ここで結論を急がないでください。
p = 0.06063 は「差がないことを示した」わけではありません。95%信頼区間は -0.1710156 から 7.5710156 で、下限がわずかに0を下回っているだけです。この区間は「差が0から7.6くらいまでのどこかにありそうだ」と言っており、大きな差の可能性を十分に含んでいます。有意でなかった=差がない、ではない。p値と信頼区間はセットで読んでください。
なぜ有意にならなかったのか。理由のひとつは、前の節で見た層別の構造です。用量をまとめて2群にしてしまったため、用量による大きなばらつきが差を見えにくくしています。p値や信頼区間そのものの考え方をもう一段深く知りたい方は、信頼区間とp値の関係を図解で理解するが図付きで整理しています。
結果は画面に出て終わりではない
Rの解析結果は、画面に印刷されて消えるものではありません。名前を付けて保存しておくと、あとから中身を部品ごとに取り出せます。
res <- t.test(len ~ supp, data = tg)
res$p.value
res$conf.int
round(res$estimate, 2)
[1] 0.06063451
[1] -0.1710156 7.5710156
attr(,"conf.level")
[1] 0.95
mean in group OJ mean in group VC
20.66 16.96
res <- t.test(...) として結果を res に入れると、res$p.value でp値だけ、res$conf.int で信頼区間だけを取り出せます。$ は「この中のこの部品」を指す記号です。round(res$estimate, 2) のように、取り出したうえで丸めることもできます。
これは地味ですが決定的に重要な性質です。数値を目で読んで報告書に打ち直す作業は、必ずどこかで写し間違えます。取り出せるということは、報告用の表をコードで組み立てられ、データが差し替わっても同じコードを流すだけで全部が更新できるということです。「表を手で転記しない」自動化の第一歩が、この $ にあります。
どんな部品が入っているかは names(res) で一覧できます。t.test() に限らず、Rの多くの解析関数が結果をひとまとまりのオブジェクトとして返すので、この確認のしかたはそのまま使い回せます。
Rの既定はWelch法である
出力の1行目が Welch Two Sample t-test だったことに気づいたでしょうか。教科書で最初に習うのは2群の分散(ばらつき)が等しいと仮定するStudentのt検定ですが、Rの t.test() は既定でその仮定を置かないWelch法を使います。教科書どおりにするには引数を1つ足します。
t.test(len ~ supp, data = tg, var.equal = TRUE)
Two Sample t-test
data: len by supp
t = 1.9153, df = 58, p-value = 0.06039
alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
95 percent confidence interval:
-0.1670064 7.5670064
sample estimates:
mean in group OJ mean in group VC
20.66333 16.96333
見出しが Two Sample t-test に変わり、自由度が 55.309 から 58 になりました。Welch法では自由度が小数になるのが特徴です。p値は 0.06063 から 0.06039 へほとんど動いていません。今回は各群30匹ずつでばらつきも近いため、どちらでもほぼ同じ答えになります。
とはいえ、群の人数が大きく違い、ばらつきも違う場面では結果がはっきり変わります。「手計算やExcelと数値が合わない」と悩む初心者の原因は、たいていこの既定の違いです。まず出力の1行目でどちらが走ったかを確認してください。どちらの方法をいつ選ぶべきかという判断や、自由度の決まり方はt検定の数理的導出とRによる実装例で扱っています。ここでは「Rの既定はWelchである」という事実だけ持ち帰れば十分です。
グラフでデータを見る ― ヒストグラム・箱ひげ図・散布図
数値の要約は便利ですが、平均が同じでも散らばり方や山のかたちはまったく違うことがあります。そこで役に立つのが図です。Rには追加インストールなしで使える作図の関数がそろっており、まず覚えるべきものは3つ。分布を見る hist()、2つのグループを比べる boxplot()、2つの数値の関係を見る plot() です。ここまで使ってきた tg で順に描きます。
分布のかたちを見る ― ヒストグラム
ヒストグラムは、数値を一定の幅の区間に区切り、各区間に何個のデータが入ったかを棒の高さで表した図です。まずは歯の長さ len の分布を見てみましょう。
hist(tg$len,
breaks = 10,
col = "#AED6F1",
border = "white",
main = "Distribution of tooth length",
xlab = "len (tooth length)",
ylab = "Frequency")
1行目の tg$len が描く対象で、2行目以降は見た目の調整です。breaks = 10 は「だいたい10本くらいの棒に」という目安で、Rはこれを参考にきりのよい区切り幅に調整します。col は棒の色、border は枠線の色、main はタイトル、xlab と ylab は軸のラベル。カッコの中に引数をカンマで区切って並べるのが、Rの作図の基本のかたちです。
出てきた図を見てください。len は4.2から33.9まで広く散らばり、10〜15付近と20〜27付近に山がある二こぶ気味の形で、左右対称の釣鐘型ではありません。平均18.81、標準偏差7.65という2つの数字だけでは、平均のまわりに1つの山があるデータを思い浮かべてしまいます。
この二こぶには理由があります。用量0.5、1、2 mg/dayの3群が混ざったデータだからです。用量が低い群は歯が短く、高い群は長い。性質の違う集団を1枚にまとめると、こうして山が複数できます。分布のかたちがおかしいと思ったら、まず「別の集団が混ざっていないか」を疑ってください。
2群を並べて比べる ― 箱ひげ図
箱ひげ図は、データの中央値・四分位数・おおよその範囲を1つの箱とヒゲで表した図で、グループごとの分布を横に並べて比べるのに向いています。投与形態 supp で len を比べてみましょう。
boxplot(len ~ supp,
data = tg,
col = c("#AED6F1", "#F9E79F"),
main = "Tooth length by supplement type",
xlab = "supp (supplement type)",
ylab = "len (tooth length)")
注目してほしいのは1行目の len ~ supp です。t.test() にも出てきた「左側の値を右側のグループで分けて見る」という書き方で、Rでは検定でも作図でも同じ文法が使えます。col に c() で2色を渡しているのは、箱が2つできるからです。
図を見ると、OJの箱のほうが上にあり、中央値を示す太い横線も高い位置にあります。平均20.66と16.96、OJが3.7高いという結果と向きは一致します。しかし同時に、2つの箱は大きく重なっています。箱の縦の長さ、つまり群内のばらつきのほうが、2群の中央値の差より大きいのです。
この見え方は、t検定の結果とぴったり対応しています。p値0.06063 は5%を下回らず、95%信頼区間 -0.17〜7.57 は下限がわずかに0を割り込んだ状態でした。図の重なり具合と数値が同じことを言っています。数字と図の印象が食い違ったときは、どちらかの読み違いを疑ってください。 ヒゲの長さや外れ値の点の決め方といった細かい仕様は箱ひげ図をSASとRで実装する方法で詳しく扱っています。
3つ目の変数を重ねる ― 散布図
散布図は、2つの数値の関係を点で表した図です。ここでは点の色で3つ目の情報も重ねます。横軸に用量 dose、縦軸に len、supp を色で塗り分けましょう。
plot(tg$dose, tg$len,
pch = 19,
col = ifelse(tg$supp == "OJ", "#2E86C1", "#E67E22"),
main = "Dose vs tooth length",
xlab = "dose (mg/day)",
ylab = "len (tooth length)")
legend("topleft", legend = c("OJ", "VC"),
col = c("#2E86C1", "#E67E22"), pch = 19, bty = "n")
pch = 19 は点の形の指定で、19番は塗りつぶした丸です。色を決める ifelse() は「条件、正しいときの値、正しくないときの値」を順に書く関数で、60個の点を「supp がOJなら青、そうでなければ橙」と判定しています。legend() は凡例を追加する関数で、"topleft" は左上、bty = "n" は枠なしの指定。plot() のあとに legend() を重ねる順番も覚えておいてください。
この図が、記事全体でいちばん情報量の多い1枚です。用量が上がるほど len が上がる傾向がはっきり見え、用量0.5と1では青(OJ)の点が橙(VC)より上に位置します。ところが2 mgでは青と橙が混ざり合ってしまいます。用量別の平均が26.06(OJ)と26.14(VC)でほぼ同じだった集計結果と完全に一致します。
全体平均で見た3.7の差は、主に低い用量で作られたものでした。それを用量で分けずにt検定にかけたので差が薄まり、p値が0.06まで上がった、という筋書きです。図は解析の飾りではなく、結果を解釈するための材料です。
最後に、描いた図をファイルとして残す方法です。作図のコードを png() と dev.off() ではさむと画像ファイルとして保存されます。本記事の3枚も、この設定で生成した実物です。
png("fig1_hist.png", width = 1200, height = 800, res = 150)
hist(tg$len,
breaks = 10,
col = "#AED6F1",
border = "white",
main = "Distribution of tooth length",
xlab = "len (tooth length)",
ylab = "Frequency")
dev.off()
png() でファイルを開き、作図し、dev.off() で閉じる。この3行がセットです。width と height は画像の大きさ(画素数)、res は解像度で、資料に貼るなら150前後にすると文字がつぶれません。dev.off() を忘れると、ファイルが開きっぱなしになって中身が書き込まれず、画面にも図が出なくなります。「グラフが急に表示されなくなった」ときは、まず dev.off() を1回実行してみてください。つまずきやすいエラーと実務でのポイント
Rを始めた人が最初の1週間で挫折する理由の大半は、統計が難しいからではありません。赤い英語のエラーが出て、何を言われているのか分からず手が止まるからです。しかしRのエラーは書式が決まっていて、種類も多くありません。読み方さえ知っていれば、エラーは「どこで、なぜ止まったか」を教えてくれる案内板になります。ここでは必ず出会う4つを、わざと発生させて読んでみます。
エラーは「Error in ◯◯ :」から読む
Rのエラーは、多くの場合「Error in 失敗した場所 : 失敗した理由」という形です。まず in のうしろを見て、次にコロンのうしろを見る。この順番だけ守れば、英語が苦手でも8割は自力で解決できます。
1つ目は、いちばん多いタイプミスです。
mena(x)
Error in mena(x) : could not find function "mena"
in のうしろに mena(x) とあり、コロンのうしろは「mena という関数が見つかりません」。mean を打ち間違えただけです。関数名が見つからないと言われたら、まずスペルを疑ってください。 Rは似た名前を勝手に補ってはくれません。
2つ目は、数値のつもりが文字列だった場合です。
"5" + 1
Error in "5" + 1 : non-numeric argument to binary operator
「二項演算子に数値でない引数が渡された」という意味です。二項演算子とは + のように左右2つの値を受け取る記号のこと。ダブルクォートで囲むとRは文字として扱うので、足し算ができません。Excelから読み込んだ列に空白や記号が紛れると、列全体が文字列になってこのエラーが出ます。
3つ目は、パッケージまわりです。
library(tidyverseee)
Error in library(tidyverseee) : there is no package called 'tidyverseee'
「そんな名前のパッケージはありません」。原因は名前の打ち間違いか、そもそもインストールしていないかです。install.packages("パッケージ名") は一度だけ、library() は毎回、という二段構えを思い出してください。
4つ目は、ファイルが見つからない場合です。
read.csv("no_such_file.csv")
Error in file(file, "rt") : cannot open the connection
Calls: read.csv -> read.table -> file
In addition: Warning message:
In file(file, "rt") :
cannot open file 'no_such_file.csv': No such file or directory
Error in file(file, "rt") は、read.csv() の内部でファイルを開く関数まで進んで止まった、という意味です。2行目の Calls: は「read.csv が read.table を呼び、それが file を呼んだ」という道順。最後の行に本当の原因、「そのようなファイルもディレクトリもありません」が書いてあります。ファイル名の打ち間違いか、ファイルが作業ディレクトリの外にあるかです。getwd() で現在地を確認するのが最短の対処になります。
4つ目の出力には
Warning message という行も混ざっています。Error と Warning はまったく別物です。Error は処理がそこで止まった合図で、結果は作られていません。Warning は「処理は進んだけれど気になることがあった」という報告で、結果自体は作られています。怖いのはむしろ後者。止まらないぶん気づかずに進み、あとで数値が合わない原因になります。内容を読んで納得してから次へ進んでください。ここまでの4つを、原因と対処の形でまとめておきます。
| エラーの中身 | よくある原因 | 最初にやること |
|---|---|---|
| could not find function | 関数名のタイプミス。または必要なパッケージ未読み込み | スペルを見直す。パッケージの関数なら library() を確認 |
| non-numeric argument to binary operator | 数値のつもりの値が文字列になっている | class() や str() で型を確認。読み込み元の空白や記号を疑う |
| there is no package called … | パッケージ名の間違い、またはインストール漏れ | 名前を確認し、install.packages() を一度だけ実行 |
| cannot open the connection | ファイル名の間違い、または作業ディレクトリの外にある | getwd() で現在地を確認し、置き場所か指定を直す |
型が変わっていないかを確かめる
エラーが出ないのに結果がおかしい、という場合の犯人はたいてい「型」です。型とは、その列が数値なのか文字なのかというデータの種類のこと。ここで気になる点を確認しておきましょう。
class(tg$supp)
class(tg$dose)
class(factor(tg$supp))
levels(factor(tg$supp))
[1] "character"
[1] "numeric"
[1] "factor"
[1] "OJ" "VC"
tg$supp は character、つまり文字列だと出ました。しかし前半で str(ToothGrowth) を見たとき、supp は Factor w/ 2 levels だったはずです。理由は単純で、一度CSVに書き出して読み直したから。CSVは中身の文字だけを保存する形式なので、型の情報は残りません。 読み込み時にRが改めて推測し、文字列だと判断したわけです。
factor() を通すと factor 型になり、levels() でその水準、つまり「取りうる値の一覧と並び順」が確認できます。ここでは "OJ" "VC" の順。今回は文字列のままでも箱ひげ図もt検定も動きましたが、どのグループを基準に比較するかなど水準の順序が意味を持つ場面では factor() で明示すべきです。次でその場面が出てきます。
次の一歩 ― 用量も一緒にモデルに入れる
最後に、ここまでの結果が次にどうつながるかを示しておきます。t検定では投与形態だけを見ていましたが、散布図で分かったとおり len は用量にも大きく左右されていました。両方を同時に扱う道具が回帰モデルです。
tg$dose_f <- factor(tg$dose)
fit <- lm(len ~ supp + dose_f, data = tg)
summary(fit)
Call:
lm(formula = len ~ supp + dose_f, data = tg)
Residuals:
Min 1Q Median 3Q Max
-7.085 -2.751 -0.800 2.446 9.650
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.4550 0.9883 12.603 < 2e-16 ***
suppVC -3.7000 0.9883 -3.744 0.000429 ***
dose_f1 9.1300 1.2104 7.543 4.38e-10 ***
dose_f2 15.4950 1.2104 12.802 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.828 on 56 degrees of freedom
Multiple R-squared: 0.7623, Adjusted R-squared: 0.7496
F-statistic: 59.88 on 3 and 56 DF, p-value: < 2.2e-16
suppVC の行を見てください。推定値は-3.70、p値は0.000429です。t検定では0.061でしたから、同じデータなのに結論が変わって見えます。差の大きさ3.7自体は変わっていません。変わったのは、その3.7がどれくらい確かに見えるかです。用量を一緒にモデルへ入れたことで用量による大きな上下動が取り除かれ、残りのばらつきが小さくなりました。背景のノイズが小さくなれば、同じ3.7でもはっきり浮かび上がります。t検定で有意でなかったからこの要因は効いていない、と決めつけてはいけない、という教訓がここにあります。なお factor() を挟んだのは、dose を数値のまま入れると「1 mg増えるごとに一定量ずつ増える」という直線の仮定が入るためです。
回帰モデルの中身、係数の決まり方や残差の点検は、この記事の範囲を超えます。考え方から知りたい方は「回帰分析って何をしているの?」線を引くだけで理解する入門を、3群以上の平均の比較に進みたい方は一元配置分散分析(One-way ANOVA)をご覧ください。ここでは「Rでは t.test() も lm() も len ~ supp と同じ書き方で呼べる」ことだけ持ち帰れば十分です。
参考書籍
Rは本を1冊やり切ると景色が変わります。おすすめは、まずR自体の操作に慣れる本、次にRを動かしながら統計の考え方を身につける本、最後に統計そのものを腰を据えて学ぶ本、という3段構えです。

read.csv() と aggregate() を手作業で並べた部分が、この本では一貫した書き方に置き換わります。作図の章はグラフの節の続きで、plot() より柔軟なggplot2まで進めます。プロジェクト機能の解説も、エラーの節と相性が良いところです。

関連記事・次のステップ
Rでデータを読み込み、要約し、t検定を実行し、グラフを描くところまで来たら、次は「その数値が何を意味しているか」を深める段階です。気になったところから読み進めてください。
t検定の周辺をもう一歩深く理解したい方は、同じ被験者を前後で比べる場合を扱った対応のあるt検定とウィルコクソン符号順位検定が続きになります。p値0.06063 の読み方に引っかかりを覚えた方には、p値とは何か?をおすすめします。
標本平均のばらつきという考え方そのものを納得したい方は、Rのシミュレーションで確かめる中心極限定理とは ― なぜ標本平均は正規分布に近づくのか ―が、本記事と同じく手を動かしながら読める内容です。結果が0か1かの二値データを扱いたくなったら2値変数とロジスティック回帰へ進めます。
臨床試験の実務で使われる解析まで見通しておきたい方は、経時的に繰り返し測定したデータを扱うmmrmをRで実装するや、診断精度を評価するROC曲線とAUCをRで実装するが、本記事の延長線上にあります。どちらも「データを読み込む→関数に渡す→出力を読む」という同じ流れで動いています。
まとめ
本記事では、Rを初めて使う方に向けて、インストールから最初の統計解析までを一通り実装しました。R本体とRStudioの役割の違いを整理し、c() でベクトルを作り、data.frame() で表を組み立て、head()・str()・summary() でデータの形を確認しました。その上で aggregate() による集計から、全体平均では3.7の差があるのに用量2 mgでは 26.06 と 26.14 でほぼ差がない、という発見にたどり着き、t.test() で p値0.06063 と95%信頼区間[-0.17, 7.57]を求め、3種類のグラフで同じデータを目で確かめました。
覚えていただきたいのは、個々の関数名ではありません。データを読み込み、形を確認し、要約し、検定し、図にするという一連の流れです。この流れはmmrmでも生存時間解析でもロジスティック回帰でも同じで、変わるのは真ん中で呼ぶ関数の名前だけ。この記事の範囲を自分の手で再現できれば、どの解析にも同じやり方で入っていけます。
もうひとつ強調したいのは、summary() の dose の平均1.167 に意味がなかったこと、全体平均だけでは用量2 mgでの逆転を見落とすこと、t検定のp値が0.061 でも用量を一緒にモデルに入れれば 0.000429 になったことです。Rは指示した計算を正確に返しますが、その数字を出す意味があるかどうかは判断してくれません。出力を鵜呑みにせず、分けて見る、図にして確かめる。この習慣が、ソフトウェアの操作よりも長く役に立ちます。
エラーが出ても、あなたがRを壊したのではなく、Rが「どこで困ったか」を教えてくれているだけです。Error in の後ろを読み、タイプミスを直し、もう一度実行する。この往復を繰り返すうちに、Rは確実に手に馴染んできます。まずは本記事のコードを上から順に打ち込むところから始めてみてください。











