Rのモデル式「~」をマスター!回帰分析からグラフ作成まで実例で学ぶ
最も基本的な使い方は、以下のように「YはXによって説明される」または「YはXの関数である」といった関係を表す場合です。
Y ~ X
この式は、統計モデルにおいて「目的変数(あるいは応答変数)Yは、説明変数(あるいは予測変数)Xによってモデル化される」という意味を持ちます。
具体的な例を挙げながら、いくつかの使い方を説明します。
-
単回帰モデル
もし、身長(height
)が体重(weight
)によって説明されるというモデルを考えたい場合、以下のように書きます。model <- lm(height ~ weight, data = my_data)
ここで、
lm()
関数は線形モデルをフィットさせる関数で、height
が目的変数、weight
が説明変数です。 -
重回帰モデル
複数の説明変数を使いたい場合は、+
記号でそれらを連結します。model <- lm(y ~ x1 + x2 + x3, data = my_data)
これは、「
y
はx1
、x2
、x3
によって説明される」というモデルを意味します。 -
交互作用項(Interaction Terms)
2つの変数の効果が互いに影響し合う場合、交互作用項をモデルに含めることができます。x1:x2
:x1
とx2
の交互作用のみを含めます。x1*x2
: これはx1 + x2 + x1:x2
の短縮形です。つまり、x1
、x2
の主効果と、それらの交互作用を含めます。
例:
model <- lm(y ~ age * gender, data = my_data)
これは、「
y
はage
とgender
の主効果、およびage
とgender
の交互作用によって説明される」というモデルです。 -
すべての変数を含める場合
データフレーム内の他のすべての変数を説明変数として含めたい場合、.
(ドット)を使うことができます。model <- lm(y ~ ., data = my_data)
これは、「
y
以外のmy_data
内のすべての変数がy
を説明する」というモデルです。 -
特定変数をモデルから除外する場合
my_data
内のほとんどの変数を使いたいが、特定の変数だけは除外したいという場合は、-
を使います。model <- lm(y ~ . - exclude_var, data = my_data)
これは、「
y
以外のすべての変数を使うが、exclude_var
は除く」というモデルです。 -
カテゴリカル変数の扱い
Rは、モデル式の中でチルダを使うと、因子型(factor)の変数を自動的にダミー変数に変換して扱ってくれます。これは、統計モデリングにおいて非常に便利です。
Rにおけるチルダ(~
)関連のよくあるエラーとトラブルシューティング
チルダ(~
)は主にlm()
, glm()
, lmer()
などのモデル関数で使われる**モデル式(formula)**の一部として機能します。そのため、チルダ関連のエラーは、通常、モデル式の記述ミスや、式内で参照されているデータの問題に起因します。
変数名の間違い (object not found または could not find variable)
エラーの例
# 間違った変数名
data_df <- data.frame(height_cm = c(160, 170, 180), weight_kg = c(55, 65, 75))
lm(height ~ weight_kg, data = data_df) # 'height' が存在しない
# エラーメッセージの例:
# Error in eval(predvars, data, env) : object 'height' not found
原因
モデル式の中で指定した変数名が、実際に使用するデータフレーム内に存在しない場合に発生します。これは、スペルミス、大文字・小文字の区別(Rは大文字・小文字を区別します)、またはデータフレームへの読み込みミスなどが原因です。
トラブルシューティング
- データフレームの確認
目的のデータが正しくデータフレームに読み込まれているかを確認します。 - スペルミスの修正
モデル式内の変数名がデータフレーム内の変数名と完全に一致するように修正します。 - 変数名の確認
names(data_df)
やhead(data_df)
を使って、データフレーム内の実際の変数名を確認してください。
データが指定されていない、または間違ったデータが指定されている
エラーの例
# data引数が欠けている、または間違っている
# height_cm, weight_kg がグローバル環境にない場合
lm(height_cm ~ weight_kg)
# エラーメッセージの例:
# Error in eval(predvars, data, env) : object 'height_cm' not found
原因
モデル式を評価する際に、Rが指定された変数を見つけられない場合に発生します。これは主に、data
引数にデータフレームが正しく指定されていないか、あるいはそもそもdata
引数を使用せずに、変数がグローバル環境に存在しない場合に起こります。
トラブルシューティング
- 変数の存在確認
ls()
やexists()
を使って、変数が現在のRセッションに存在するかを確認します。もしデータフレームを使わずに個別のベクトルとして変数がある場合は、その変数が利用可能であることを確認します。 - data引数の指定
モデル関数(例:lm()
,glm()
)のdata
引数に、目的のデータフレームを明示的に指定してください。lm(height_cm ~ weight_kg, data = data_df)
モデル式の構文エラー
エラーの例
# 不要なカンマ、または演算子の誤用
lm(Y ~ X1, + X2, data = my_data) # カンマが不要
# エラーメッセージの例:
# Error: unexpected ',' in "lm(Y ~ X1, + X2"
トラブルシューティング
- ドット(.)の誤用
Y ~ .
は「Y以外のすべての変数」を意味します。これを意図しない場合や、別の意味でドットを使おうとしている場合は注意が必要です。 - 括弧のバランス
特に複雑な式の場合、括弧の開始と終了が一致しているか確認します。 - 演算子の確認
- 説明変数を複数指定する場合は
+
を使用します(例:Y ~ X1 + X2
)。 - 交互作用項は
:
または*
を使用します(例:Y ~ X1:X2
またはY ~ X1*X2
)。 - 除外したい変数は
-
を使用します(例:Y ~ . - X3
)。
- 説明変数を複数指定する場合は
変数の型が不適切 (non-numeric argument to a binary operator など)
エラーの例
# 数値型であるべき変数が文字列型になっている
data_bad_type <- data.frame(score = c(10, 15, "twenty"), study_hours = c(2, 3, 4))
lm(score ~ study_hours, data = data_bad_type)
# エラーメッセージの例:
# Error in model.frame.default(formula = score ~ study_hours, data = data_bad_type, :
# variable 'score' is not numeric
原因
モデル関数は、通常、数値データまたは因子(factor)データを期待します。文字列(character)型や、予期しない型の変数がモデル式に含まれている場合にエラーが発生します。
トラブルシューティング
- 型の変換
必要な場合、as.numeric()
,as.factor()
,as.character()
などを使って変数の型を変換します。- 数値に変換できない文字が含まれていると、
as.numeric()
はNA
(欠損値)を生成することがあります。これには注意が必要です。
- 数値に変換できない文字が含まれていると、
- 変数の型を確認
str(data_df)
やclass(data_df$variable_name)
を使って、各変数の型を確認します。
欠損値(NA)の扱い
エラーの例
直接的なエラーは少ないですが、モデルの結果が意図しないものになることがあります。
data_na <- data.frame(Y = c(1, 2, NA, 4), X = c(10, 20, 30, 40))
lm(Y ~ X, data = data_na) # NAを含む行は自動的に除外される
原因
多くのRのモデル関数は、デフォルトで欠損値(NA
)を含む行をモデル構築から自動的に除外します(na.action = na.omit
がデフォルトの場合が多い)。これにより、意図しないデータ数の減少や、モデルの解釈に影響が出ることがあります。
トラブルシューティング
- na.action引数の利用
モデル関数にna.action
引数を指定して、欠損値の処理方法を明示的に制御できます。lm(Y ~ X, data = data_na, na.action = na.exclude)
- 欠損値の処理方法
na.omit(data_df)
: 欠損値を含む行をデータフレームから完全に削除します。na.exclude(data_df)
:na.omit
と似ていますが、予測などの後続処理で欠損値の場所を保持します。na.fail
: 欠損値があればエラーを出します。- 欠損値の補完(Imputation):
mice
やAmelia
などのパッケージを使って、欠損値を統計的に補完することも検討します。
- 欠損値の確認
summary(data_df)
やcolSums(is.na(data_df))
を使って、データフレーム内の欠損値の数を把握します。
パッケージがロードされていない
エラーの例
# lmer() 関数はlme4パッケージにある
# library(lme4) を実行していない場合
lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
# エラーメッセージの例:
# Error in lmer(Reaction ~ Days + (1 | Subject), data = sleepstudy) :
# could not find function "lmer"
原因
チルダ自体がエラーの原因ではありませんが、特定のモデル関数(例: lmer
や glm.nb
など)を使用する場合、それらが定義されているパッケージを事前にロードしていないと、「関数が見つからない」というエラーが発生します。
- パッケージのインストール
もしパッケージがインストールされていない場合は、install.packages("package_name")
でインストールします。 - パッケージのロード
使用する関数が含まれるパッケージをlibrary()
関数でロードします。library(lme4) # lmer()を使う場合
- 検索エンジンの活用
エラーメッセージをそのままコピーしてGoogleなどで検索すると、同じ問題に直面した他の人の解決策が見つかることが多いです。Stack Overflowなどは特に役立ちます。 - コードの簡略化
複雑なモデル式でエラーが出る場合、まずシンプルなモデル式(例:Y ~ X1
)で試してみて、徐々に要素を追加していくと、どこで問題が発生したか特定しやすくなります。 - エラーメッセージを読む
Rのエラーメッセージはしばしば具体的な情報を含んでいます。何が間違っているのかを理解する手がかりになります。
ここでは、チルダ(~
)を使ったRプログラミングの具体例をいくつかご紹介します。
線形回帰モデル (lm())
最も一般的な用途の一つが、線形回帰モデルの定義です。lm()
関数は、目的変数(反応変数)が一つ以上の説明変数(予測変数)によってどのように説明されるかをモデル化します。
基本形: 目的変数 ~ 説明変数
例1: 単純な線形回帰
データセットmtcars
(自動車に関する組み込みデータセット)を使って、車の馬力(hp
)が燃費(mpg
)にどう影響するかをモデル化します。
# データをロード
data(mtcars)
# 単純な線形回帰モデルの作成
# mpg(燃費)を目的変数、hp(馬力)を説明変数とする
model_simple <- lm(mpg ~ hp, data = mtcars)
# モデルの要約を表示
summary(model_simple)
# 結果の解釈:
# hpが増加するとmpgが減少するという負の関係が見られます。
# Adjusted R-squaredは、モデルがデータの変動をどの程度説明しているかを示します。
この例では、mpg ~ hp
という式が「mpg
はhp
によって説明される」という意味になります。
例2: 重回帰モデル
複数の説明変数を使ってモデルを作成する場合、+
記号で連結します。
# mpgを目的変数、hpとwt(車重)を説明変数とする重回帰モデル
model_multiple <- lm(mpg ~ hp + wt, data = mtcars)
# モデルの要約を表示
summary(model_multiple)
# 結果の解釈:
# hpとwtの両方がmpgに独立して影響を与えることがわかります。
# 複数の要因を同時に考慮することで、より精密なモデルを構築できます。
mpg ~ hp + wt
は、「mpg
はhp
とwt
によって説明される」ことを示します。
交互作用項の追加
変数が互いに影響し合う(交互作用がある)場合、モデル式に交互作用項を含めることができます。
変数A * 変数B
:変数A
の主効果、変数B
の主効果、そして変数A
と変数B
の交互作用を含める(A + B + A:B
の短縮形)。変数A : 変数B
:変数A
と変数B
の交互作用のみを含める。
例3: 交互作用項を含むモデル
mtcars
データセットで、hp
とam
(オートマチックかマニュアルか)の交互作用がmpg
に与える影響を見てみましょう。
# amを因子型に変換(カテゴリカル変数として扱うため)
mtcars$am <- as.factor(mtcars$am) # 0 = auto, 1 = manual
# hpとamの交互作用を含むモデル
# mpg ~ hp + am + hp:am と同等
model_interaction <- lm(mpg ~ hp * am, data = mtcars)
# モデルの要約を表示
summary(model_interaction)
# 結果の解釈:
# hpとamの交互作用項が有意であれば、馬力と燃費の関係が、
# オートマチック車とマニュアル車で異なることを示唆します。
mpg ~ hp * am
という式は、「mpg
はhp
とam
のそれぞれの効果に加え、hp
とam
が組み合わさった時の特別な効果(交互作用)によって説明される」という意味です。
すべての説明変数を含める (.)
データフレーム内の目的変数以外のすべての変数を説明変数として含めたい場合、ドット(.
)を使います。
例4: すべての変数を使用
# mpg以外のmtcars内のすべての変数を説明変数とする
model_all_vars <- lm(mpg ~ ., data = mtcars)
# モデルの要約を表示
summary(model_all_vars)
# 注意: 変数が多いとモデルが複雑になりすぎることがあります。
# また、多重共線性の問題にも注意が必要です。
mpg ~ .
は非常に便利ですが、必要以上に複雑なモデルを作成したり、解釈を困難にしたりする可能性があるため、注意して使用する必要があります。
特定の変数を除外する (-)
すべての変数を含めたいが、特定の変数だけは除外したい場合、ハイフン(-
)を使います。
例5: 特定の変数を除外
mtcars
データセットでmpg
を目的変数とし、disp
(排気量)以外のすべての変数を使いたい場合。
# mpg以外の変数すべてを使用し、dispだけを除外する
model_exclude_var <- lm(mpg ~ . - disp, data = mtcars)
# モデルの要約を表示
summary(model_exclude_var)
mpg ~ . - disp
は、「mpg
以外のすべての変数を使うが、その中からdisp
は除く」という意味です。
チルダは、統計モデルだけでなく、データの可視化でも変数の関係を指定するために使われることがあります。特にggplot2
パッケージのfacet_wrap()
やfacet_grid()
関数でよく使われます。
例6: facet_wrap()
での使用
変数に基づいて複数のサブプロットを作成する際にチルダを使います。
# ggplot2パッケージをロード
library(ggplot2)
# hpとmpgの関係をam(オートマチック/マニュアル)で分割して表示
ggplot(mtcars, aes(x = hp, y = mpg)) +
geom_point() +
facet_wrap(~ am) # amの値ごとにグラフを分割
# 結果の解釈:
# amが0(オートマチック)と1(マニュアル)の2つのパネルが生成され、
# それぞれのパネルでhpとmpgの関係が見られます。
~ am
は、「am
の値に基づいてデータを分割し、それぞれのサブプロットを作成する」という意味です。
Rのチルダ(~
)は、データ分析において変数間の関係を簡潔かつ強力に表現するための記号です。モデル式の定義からグラフ作成まで、幅広いRの機能で利用されており、Rを使ったデータ分析の基本的なスキルとして習得しておくことが重要です。
モデル式を使わない回帰分析(行列演算)
統計モデルの根底には線形代数と行列演算があります。Rでは、lm()
などの関数が内部的にこれらを処理してくれますが、直接行列を使って回帰係数を計算することも可能です。これは特に、カスタムの統計手法を実装したい場合や、モデルの動作原理を深く理解したい場合に有用です。
代替方法の例: 最小二乗法の直接計算
線形回帰モデル Y=Xbeta+epsilon における係数 beta は、最小二乗法により hatbeta=(XTX)−1XTY で計算できます。
# データ準備
data(mtcars)
Y <- mtcars$mpg # 目的変数
X <- as.matrix(cbind(1, mtcars$hp)) # 説明変数 (切片のために1の列を追加)
# 係数の計算
beta_hat <- solve(t(X) %*% X) %*% t(X) %*% Y
# 結果の表示
print(beta_hat)
# lm()の結果と比較
model_lm <- lm(mpg ~ hp, data = mtcars)
print(coef(model_lm))
この方法では、チルダを使ったモデル式は一切登場せず、純粋な行列演算で回帰係数を導出しています。ただし、これは非常に基礎的なケースであり、交互作用項やカテゴリ変数などを扱う場合は、X
行列を適切に構築する必要があります。
データフレームの列名を直接指定
単純な操作や、モデルを構築するのではなく、データフレーム内の特定の列を直接操作したい場合は、チルダを使用する必要はありません。
代替方法の例1: 変数間の相関計算
data(mtcars)
# チルダを使わない方法
cor(mtcars$mpg, mtcars$hp)
# チルダを使った方法 (これは相関分析関数では一般的ではない)
# cor(~ mpg + hp, data = mtcars) のように使われる関数もありますが、
# 汎用的ではありません。
代替方法の例2: 新しい変数の作成
data(mtcars)
# チルダを使わない方法
mtcars$hp_per_weight <- mtcars$hp / mtcars$wt
head(mtcars)
# チルダは変数作成には通常使われません。
プログラムによる動的なモデル式の構築
モデル式は通常、Rの構文として直接記述しますが、文字列としてモデル式を構築し、それをRの関数に渡すことも可能です。これは、ユーザーの入力に基づいてモデルの変数を動的に変更したい場合や、多数のモデルをループで生成したい場合に非常に役立ちます。
代替方法の例: paste()
とas.formula()
の組み合わせ
data(mtcars)
# 説明変数リストを文字列として定義
explanatory_vars <- c("hp", "wt", "cyl")
# モデル式を文字列として構築
formula_string <- paste("mpg ~", paste(explanatory_vars, collapse = " + "))
print(formula_string)
# 文字列をモデル式オブジェクトに変換
dynamic_formula <- as.formula(formula_string)
print(dynamic_formula)
# モデルの実行
model_dynamic <- lm(dynamic_formula, data = mtcars)
summary(model_dynamic)
この方法では、直接mpg ~ hp + wt + cyl
と書く代わりに、as.formula()
関数と文字列操作を組み合わせてモデル式を生成しています。これにより、例えば特定の条件に基づいて説明変数を追加・削除するような柔軟なプログラミングが可能になります。
汎用的なデータ操作パッケージ(Tidyverseなど)
dplyr
やtidyr
といったTidyverseエコシステムのパッケージは、データ操作の強力なツールを提供しますが、これらは通常、チルダを直接使うことはありません。代わりに、パイプ演算子(%>%
)と関数を使って、データの変換や集計を行います。
代替方法の例: dplyr
でのグループ化と要約
library(dplyr)
data(mtcars)
# am(オートマチック/マニュアル)でグループ化し、mpgの平均を計算
# チルダは使用しない
summary_mpg_by_am <- mtcars %>%
group_by(am) %>%
summarise(
mean_mpg = mean(mpg),
sd_mpg = sd(mpg)
)
print(summary_mpg_by_am)
ここでは、group_by()
やsummarise()
といった関数内で列名を直接指定しており、チルダは使われていません。Tidyverseは、データ操作をより直感的かつ読みやすい方法で記述するためのものです。
- カスタムアルゴリズムの実装や、Rの内部動作の理解を深めたい場合は、行列演算などによる直接的な実装が役立ちます。
- 基本的なデータ操作や探索的データ分析では、チルダを使わない直接的な列アクセスやTidyverseの関数が適しています。
- 動的なモデル構築が必要な場合は、
as.formula()
と文字列操作を組み合わせる方法が非常に強力です。 - ほとんどの場合、チルダを使ったモデル式が最適です。 統計モデルの記述には、簡潔で読みやすく、Rの標準的な方法だからです。Rの多くの統計関数は、このモデル式を前提として設計されています。