今回は、線形混合モデル（Linear Mixed Model, LMM）のお話です。学生の皆さんが線形混合モデルを使った分析のイメージをつかむのに役立てばうれしいです。
線形混合モデルは、
「たくさんの人からデータを集めたとき、個人差も考慮しながら、いろいろな要因（例：活動量や気温）の影響を分析したい」
ときに使える分析手法です。線形混合モデルでは、全員に共通する傾向を表す固定効果と、個人やグループで異なる特性を表すランダム効果を区別して、評価します。
固定効果とランダム効果の違いとは？
線形混合モデルは、次の2つの種類の効果を同時に扱うことができる統計モデルです：
- 固定効果（Fixed Effects）：全体に共通する傾向
- ランダム効果（Random Effects）：個人やグループごとのばらつき（共通傾向からのずれ）
固定効果（Fixed Effects）
定義：
全ての観測対象（例：全対象者）に共通して当てはまる影響や傾向を表します。
例：
ランダム効果（Random Effects）
定義：
観測対象ごとに異なる「個人差・個体差」をモデル化します。
例：
固定効果とランダム効果のまとめ
| 項目 | 固定効果 | ランダム効果 |
|---|---|---|
| 意味 | 全体に共通する傾向 | 個人やグループに依存したばらつき |
| 対象 | すべての観測対象に共通 | 各個人・各グループごとに異なる |
| 例 | 活動量や気温の効果 | 個人ごとの安静時心拍数の違い |
| 使う目的 | 「共通の傾向」を見つける | 「個人差」を適切に扱う |
とにかく試してみる：身体加速度を用いた心拍数の予測
仮想データの生成
実際に分析してみるためには具体的なデータが必要なので、ここでは乱数を使ってデータを作ります。ここでは、実験対象者は、ランダムに決めたおおよそのペースで、ジョギングを行い、そのときの身体加速度と心拍数を計測したとします。身体加速度は、高周波成分を抽出して、その標準偏差を指標として使います。10人の対象者に対し、一人当たり、20回計測をしたとしています。
| 変数名 | 内容 |
|---|---|
ID |
実験対象の個人の識別番号 |
Measurement |
個人ごとの計測回数 |
SD_Acc_HF |
加速度センサの高周波成分の標準偏差（単位：適当） |
HR |
心拍数（beats per minute） |
# 実験対象者数と計測回数
n_subjects <- 10
n_measurements <- 20
n_total <- n_subjects * n_measurements
# 個人ID（1〜10を繰り返し）
ID <- rep(1:n_subjects, each = n_measurements)
# 計測回数（1〜20を繰り返し）
Measurement <- rep(1:n_measurements, times = n_subjects)
# 加速度の高周波成分の標準偏差（0.01〜0.06 m/s²）
SD_Acc_HF <- round(runif(n_total, min = 0.01, max = 0.06), 3)
# 個人ごとのベース心拍数（平均70±5 bpm）
subject_intercepts <- rnorm(n_subjects, mean = 70, sd = 5)
intercept_per_ID <- subject_intercepts[ID]
# 個人ごとの傾き（加速度への感度、平均230±20）
subject_slopes <- rnorm(n_subjects, mean = 230, sd = 20)
slope_per_ID <- subject_slopes[ID]
# 測定誤差（平均0、標準偏差2）
error <- rnorm(n_total, mean = 0, sd = 2)
# 心拍数の生成（個人ごとの切片と傾きを使用）
HR <- round(intercept_per_ID + slope_per_ID * SD_Acc_HF + error, 1)
# データフレームの作成
df <- data.frame(
ID = factor(ID),
Measurement = Measurement,
SD_Acc_HF = SD_Acc_HF,
HR = HR
)
# 最初の10行を表示
head(df, 10)/code>Rで線形混合モデル分析（その１：ランダム切片）
先ほど生成したデータを用いてRで分析してみます。
この分析では、実験で得られた加速度センサの情報のうち、高周波成分の変動の大きさ（標準偏差）である SD_Acc_HF を用いて、心拍数 HRがどのように変化するかを調べます。高周波成分は、身体の細かい動きや強い振動を捉える指標であり、その変動が大きいほど活動が激しいと考えられます。
心拍数は個人の体質や体調によって基準値が異なるため、すべての人を同じ基準で評価してしまうと、分析結果に偏りが出る可能性があります。そこで、このモデルでは、各対象者の個人差を考慮するために、個人ごとの心拍数の「基準値のずれ」をランダム効果として導入します。このように、全体に共通する影響（固定効果）と、個人ごとの違い（ランダム効果）の両方を取り入れることができるのが、線形混合モデルの大きな特徴です。
た
今回は、 lme4 パッケージを使います。lme4 をインストールしていない場合は、
install.packages("lme4")
を実行して、パッケージをインストールして下さい。
Rで線形混合モデル分析をするには、まず lme4 パッケージを読み込んでおき、lmer() 関数を使います。モデル式は
HR ~ SD_Acc_HF + (1 | ID)
と記述します。これは、心拍数（HR）が、加速度の高周波成分（SD_Acc_HF）によって説明されることを意味し、同時に個人（ID）ごとに異なるベースライン（切片）を許すという構造になっています。
このモデルをRで実行するには、以下のようにスクリプトを書きます。
# パッケージの読み込み
library(lme4)
model <- lmer(HR ~ SD_Acc_HF + (1 | ID), data = df)
summary(model)出力結果には、固定効果としての SD_Acc_HF の回帰係数や、実験対象者ごとのばらつきを表すランダム効果の推定値、また残差の標準偏差などが表示されます。たとえば、SD_Acc_HF の係数が230であれば、それは「高周波加速度成分の標準偏差が 0.01 だけ増加すると、心拍数は平均で 2.3 bpm 増加する」と解釈できます。
一方、ランダム効果として表示される「ID の Intercept の標準偏差」が大きければ、個人ごとに心拍数の基準値に大きな差があることを示しています。このように、固定効果はすべての対象者に共通する関係性を表し、ランダム効果はそれぞれの個人がもつ個別の特徴をモデルに組み込んでいます。
このような線形混合モデルは、集団全体の傾向と個人ごとの差異の両方を同時に扱うことができるため、生理データのように「人による違い」が大きいデータの解析に適しています。
結果の読み方
上に示した命令 summary(model) を実行すると、固定効果とランダム効果に関する結果が表示されます。ここでは以下のような出力が得られたと仮定して、その意味を解説します。
Linear mixed model fit by REML ['lmerMod'] Formula: HR ~ SD_Acc_HF + (1 | ID) Data: df REML criterion at convergence: 883.2 Scaled residuals: Min 1Q Median 3Q Max -2.39196 -0.73648 0.06453 0.65878 2.34852 Random effects: Groups Name Variance Std.Dev. ID (Intercept) 19.697 4.438 Residual 4.072 2.018 Number of obs: 200, groups: ID, 10 Fixed effects: Estimate Std. Error t value (Intercept) 71.825 1.453 49.44 SD_Acc_HF 242.107 9.966 24.29 Correlation of Fixed Effects: (Intr) SD_Acc_HF -0.239
結果の出力は、大きく分けて以下の3つの部分に分かれています。
1. 残差と収束情報
REML criterion at convergence: 883.2 Scaled residuals: Min 1Q Median 3Q Max -2.39196 -0.73648 0.06453 0.65878 2.34852
この部分は、モデルが正しく計算できたことを示す情報です。REML criterionはモデルの適合度（あてはまりの良さ）を示す値ですが、単体ではあまり意味を持ちません。他のモデルと比較する際に使用されます。
Scaled residualsは、予測値と実測値の差（残差）を標準化したもので、中央値がほぼ0に近く、±3の範囲に収まっているため、残差の分布に大きな異常はないと読み取れます。
2. ランダム効果（個人差のばらつき）
Random effects: Groups Name Variance Std.Dev. ID (Intercept) 19.697 4.438 Residual 4.072 2.018 Number of obs: 200, groups: ID, 10
ここはランダム効果の部分で、個人（ID）ごとの違いに関する情報です。
ID (Intercept)の標準偏差 = 4.438 という値は、「心拍数の平均的な高さ」に約±4.4 bpmの個人差があることを意味します。
また、Residualの標準偏差 = 2.018 は、同じ個人でも、計測ごとに約±2 bpmのばらつき（誤差や日ごとの変動）があることを示しています。
「Number of obs: 200, groups: ID, 10」より、10人の対象者に対してそれぞれ20回（合計200件）の測定が行われていることが確認できます。
3. 固定効果（全体に共通する関係）
Fixed effects: Estimate Std. Error t value (Intercept) 71.825 1.453 49.44 SD_Acc_HF 242.107 9.966 24.29
この部分は、**全体に共通する傾向（固定効果）**の推定結果です。
切片（Intercept） = 71.825 は、SD_Acc_HF = 0 のときの心拍数の理論値を意味します。実際には SD_Acc_HF が0になることはありませんが、モデルの数式上の出発点として解釈されます。
SD_Acc_HF の係数 = 242.107 は、加速度センサの高周波成分の標準偏差が 1.00 増加すると、心拍数が平均で約 242 bpm 増加すると推定されています。ただし、SD_Acc_HF の値は実際には 0.01 ～ 0.06 程度なので、現実的には「0.01増加すると約2.42 bpm 心拍数が増える」という意味になります。
t値 = 24.29 は、p値の評価に使える統計量であり、この効果が統計的に非常に有意であることを示唆しています。
結論：このモデルからわかること
この分析結果から、次のようなことが読み取れます：
- 心拍数のベースラインには約±4.4 bpmの個人差がある。
- 高周波加速度成分（SD_Acc_HF）の変動が心拍数に強く影響しており、わずかな増加でも心拍が有意に上昇する。
- 同一たいたいでも、日ごとに約±2 bpmのばらつきがあることから、一定の自然変動が存在する。
これらをふまえると、身体の細かい動き（高周波の加速度成分）を使って、心拍数の変化を予測するモデルが有効であることが確認できます。そして、このモデルは単純な回帰ではなく、個人差を吸収しつつ、共通の傾向を捉えるという点で非常に実用的です。
Rで線形混合モデル分析（その２：ランダム切片＆傾き）
上の例では、全ての実験対象者が「高周波加速度成分（SD_Acc_HF）」に対して同じ傾き（影響）を持つと仮定していました。しかし、実際にはある人は加速度の変化に敏感に心拍が反応し、別の人はそうでもない場合があります。
そこで今回は、個人ごとに SD_Acc_HF に対する心拍数の応答（傾き）も異なると仮定し、切片と傾きの両方にランダム効果を導入したモデルを構築します。
以下がRで命令です。
# パッケージの読み込み
library(lme4)
# 傾きにもランダム効果を導入
model_random_slope <- lmer(HR ~ SD_Acc_HF + (SD_Acc_HF | ID), data = df)
# 結果の表示
summary(model_random_slope)モデル式
HR ~ SD_Acc_HF + (SD_Acc_HF | ID)
は次のような構造を持っています：
HR ~ SD_Acc_HF：全体に共通する固定効果（平均的な傾きと切片）(SD_Acc_HF | ID)：個人ごとに「切片と傾きの両方」にランダムなずれを許す（つまり個人差を考慮）
これにより、以下が可能になります：
- Aさんは「安静時心拍数が高く、加速度に敏感」
- Bさんは「安静時心拍数は低く、加速度に鈍感」 といった**個人ごとのパターンをモデルが評価します。
以下に、あなたが提供された出力結果に基づき、読み方の解説を新たに書き直しました。内容はランダム効果と固定効果の順に説明し、結論で要点をまとめています。
出力結果の読み方
1. ランダム効果（個人差）
モデル出力の中で「Random effects」として表示される部分は、作業者ごとに異なる心拍数の傾向や、加速度への反応の違いを示しています。
Random effects: Groups Name Variance Std.Dev. Corr ID (Intercept) 14.711 3.835 SD_Acc_HF 2365.355 48.635 0.17 Residual 3.601 1.898 Number of obs: 200, groups: ID, 10
この出力から、まず (Intercept) の標準偏差が約 3.8 であることが分かります。これは作業者によって安静時の平均心拍数が ±3.8 bpm 程度ばらついていることを意味します。また、SD_Acc_HF に対する傾きの標準偏差が非常に大きく、約 48.6 あることから、加速度の変化に対する心拍数の反応には個人差が大きいことが読み取れます。
切片と傾きの相関は 0.17 と比較的小さく、わずかに正の関係があります。これは、安静時の心拍数が高い人ほど、加速度の影響を受けやすい傾向があることを示していますが、強い相関とは言えないため、慎重に解釈する必要があります。
また、Residual の標準偏差は 1.898 であり、同じ作業者内での計測誤差や日々の心拍数の自然なばらつきがこの程度あることを意味しています。
2. 固定効果（全体の傾向）
続いて、「Fixed effects」の出力は全体に共通する平均的な傾向を示しています。
Fixed effects: Estimate Std. Error t value (Intercept) 71.931 1.264 56.91 SD_Acc_HF 237.752 18.050 13.17
切片（Intercept）の推定値は約 71.9 であり、これは加速度が0だった場合の心拍数の基準値に相当します。実際の加速度が0になることはほぼありませんが、モデルの数式上の出発点としての意味を持ちます。
SD_Acc_HF の係数は 237.8 となっており、これは加速度の標準偏差が1.0（m/s²）増えると心拍数が平均で237.8 bpm 上がることを示します。実際の変化は0.01〜0.06程度なので、より現実的には「0.01増加すると心拍数が約2.38 bpm 増える」と読み替えるとよいでしょう。
t value = 13.17 という大きな値から、加速度の影響が統計的に有意であり、偶然の変動によるものではないことが示唆されます。
結論：このモデルからわかること
この結果から、加速度が増加することで心拍数も上がるという明確な傾向が確認されます。さらに、作業者ごとに心拍数のベースラインにも、加速度に対する反応の仕方にも個人差があることがモデルにより明らかになっています。相関が強くないため、安静時心拍数と加速度感度には一貫した関係は見られませんが、わずかな傾向は認められます。
このような線形混合モデルを使うことで、全体の傾向だけでなく、個人の違いを丁寧に捉えながら、生理データの構造をより正確に理解することができます。ウェアラブルデバイスのように繰り返しデータが得られる場面では、非常に有用な分析手法だといえるでしょう。
おまけ：切片を固定効果とし、傾きだけにランダム効果を仮定
この場合は、以下のようにします。
# パッケージの読み込み
library(lme4)
# 切片は固定、傾きのみランダム効果
model_random_slope_only <- lmer(HR ~ SD_Acc_HF + (0 + SD_Acc_HF | ID), data = df)
# モデル結果を表示
summary(model_random_slope_only)【発展】曜日の効果を検証
（注意：サンプルデータを生成するスクリプトは最後に掲載しました。）
ここでは、心拍数（HR）が加速度センサの高周波成分（SD_Acc_HF）だけでなく、曜日の違いによっても変化するのではないかと想定して調べます。これは、ソーシャルジェットラグの評価を意識した問題設定です。
曜日については、月曜から金曜までの5つの値がありますが、金曜日を基準（＝比べる対象）として設定します。こうすることで、他の曜日（月〜木）が金曜日と比べて心拍数にどれだけ影響を与えているのかが分かりやすくなります。
Rでこのモデルを組むには、まず曜日を「因子型」というカテゴリ変数にして、金曜日を基準に設定します。加速度と曜日を固定効果として、さらに個人ごとに切片（安静時の心拍数）、加速度への反応、曜日への反応の違いをランダム効果として含めます。以下がそのRスクリプトです。
# パッケージの読み込み
library(lme4)
# 曜日を因子に変換し、金曜日を基準にします
df$Weekday <- factor(df$Weekday, levels = c("Mon", "Tue", "Wed", "Thu", "Fri"))
df$Weekday <- relevel(df$Weekday, ref = "Fri")
# 線形混合モデルの構築
model <- lmer(HR ~ SD_Acc_HF + Weekday + (1 + SD_Acc_HF + Weekday | ID), data = df)
# 結果の表示
summary(model)結果の読み方
上のスクリプトを実行した結果、以下の結果がえられたとして、結果の読み方を説明します。
Linear mixed model fit by REML ['lmerMod'] Formula: HR ~ SD_Acc_HF + Weekday + (1 + SD_Acc_HF + Weekday | ID) Data: df REML criterion at convergence: 871.1 Scaled residuals: Min 1Q Median 3Q Max -2.66148 -0.58260 -0.06069 0.58481 2.72529 Random effects: Groups Name Variance Std.Dev. Corr ID (Intercept) 10.1855 3.1915 SD_Acc_HF 330.3732 18.1762 0.11 WeekdayMon 0.8771 0.9365 -0.15 -1.00 WeekdayTue 2.0909 1.4460 0.17 -0.33 0.32 WeekdayWed 1.2394 1.1133 0.01 -0.39 0.39 0.99 WeekdayThu 4.0000 2.0000 -0.31 -0.45 0.47 0.88 0.95 Residual 3.6654 1.9145 Number of obs: 200, groups: ID, 10 Fixed effects: Estimate Std. Error t value (Intercept) 67.8929 1.1135 60.974 SD_Acc_HF 225.7206 12.1646 18.556 WeekdayMon -1.6610 0.5236 -3.172 WeekdayTue -1.1840 0.6293 -1.881 WeekdayWed -0.2140 0.5559 -0.385 WeekdayThu -0.4799 0.7659 -0.627 Correlation of Fixed Effects: (Intr) SD_A_H WkdyMn WekdyT WkdyWd SD_Acc_HF -0.238 WeekdayMon -0.222 -0.310 WeekdayTue -0.012 -0.142 0.417 WeekdayWed -0.129 -0.156 0.458 0.721 WeekdayThu -0.323 -0.214 0.452 0.726 0.715 optimizer (nloptwrap) convergence code: 0 (OK) boundary (singular) fit: see help('isSingular')
このモデルは、心拍数（HR）が加速度センサの高周波成分（SD_Acc_HF）と曜日（Weekday）によってどう変化するかを調べるために構築されています。また、実験対象者（ID）ごとに安静時心拍数、加速度への反応、曜日による影響が異なることも考慮に入れ、それらをランダム効果として扱っています。
はじめに、ランダム効果の出力を見てみましょう。Rの出力では、以下のように表示されます。
Random effects: Groups Name Variance Std.Dev. Corr ID (Intercept) 10.1855 3.1915 SD_Acc_HF 330.3732 18.1762 0.11 WeekdayMon 0.8771 0.9365 -0.15 -1.00 WeekdayTue 2.0909 1.4460 0.17 -0.33 0.32 WeekdayWed 1.2394 1.1133 0.01 -0.39 0.39 0.99 WeekdayThu 4.0000 2.0000 -0.31 -0.45 0.47 0.88 0.95 Residual 3.6654 1.9145
この出力からわかるように、安静時心拍数（切片）には標準偏差が約3.19あり、個人差があることが確認できます。加速度への反応の違いも大きく、標準偏差は18.18と推定されています。これは、人によって加速度に対して心拍が大きく上がる場合もあれば、あまり変化しない場合もあることを表しています。
曜日の影響にも個人差があり、たとえば月曜日（WeekdayMon）の標準偏差は約0.94、木曜日（WeekdayThu）は2.00と、曜日によってばらつきの大きさが異なっています。相関係数（Corr）の部分を見ると、-1.00や0.99など極端な値が含まれており、モデルの推定がやや不安定になっていることがうかがえます。ただし、今回のデータはシミュレーションにより人工的に生成されたものであるため、このような相関が現れるのは特に問題ではありません。
次に、固定効果の出力を確認します。Rでは以下のように表示されます。
Fixed effects: Estimate Std. Error t value (Intercept) 67.8929 1.1135 60.974 SD_Acc_HF 225.7206 12.1646 18.556 WeekdayMon -1.6610 0.5236 -3.172 WeekdayTue -1.1840 0.6293 -1.881 WeekdayWed -0.2140 0.5559 -0.385 WeekdayThu -0.4799 0.7659 -0.627
この部分は、全体としての平均的な影響を表しています。まず、切片（Intercept）は約67.89であり、これは金曜日に加速度がゼロだったときの心拍数の理論値です。加速度の影響（SD_Acc_HF）は係数が225.72となっており、加速度が1単位増えると心拍数が 225.7 bpm増えるという意味になります。加速度の値は0.01〜0.06程度なので、0.01単位の変化で約 2.26 bpm の変化と考えると現実的です。
曜日の効果は金曜日を基準として、他の曜日との違いを表しています。たとえば月曜日（WeekdayMon）は -1.66 と推定されており、金曜日より心拍数が平均で 1.66 bpm低いことがわかります。この値は t 値が -3.17 と大きく、有意な差と判断できます。火曜日は -1.18 でやや低めですが、統計的に有意かどうかはやや微妙です。水曜と木曜については、差が小さく、はっきりとした傾向は見られません。
出力の最後には、次のようなメッセージが表示されています。
optimizer (nloptwrap) convergence code: 0 (OK) boundary (singular) fit: see help('isSingular')
これはモデルの収束自体は成功しているものの、一部のランダム効果の推定が不安定、または冗長であることを示しています。とくに「singular fit」は、モデルが複雑すぎて一部のパラメータがうまく区別できていない可能性を意味します。必要に応じて、ランダム効果の数を減らすか、モデルを簡素化して比較検討することが望ましいです。
このように、今回のモデルは加速度と曜日の影響を分析しながら、個人ごとの違いもきちんと捉えており、日常的な生理データの理解にとても役立つアプローチです。複雑なモデルではありますが、個人差や曜日ごとの変化をよりリアルに表現できることが、この分析の大きな強みと言えるでしょう。
おまけ：数値データの生成
データの生成は省略しようと思いましたが、念のため掲載しておきます。
# -------------------------------
# 線形混合モデル検証用データ生成
# SD_Acc_HF × Weekday × 個人のランダム効果を含む
# -------------------------------
# 人数と計測数の設定
n_subjects <- 10 # 人数
n_days <- 20 # 計測回数/人
n_total <- n_subjects * n_days
# 個人IDと計測番号
ID <- rep(1:n_subjects, each = n_days)
Measurement <- rep(1:n_days, times = n_subjects)
# 曜日（均等に割り当て、金曜を基準とする）
weekdays <- c("Mon", "Tue", "Wed", "Thu", "Fri")
Weekday <- rep(weekdays, length.out = n_total)
# 加速度の高周波成分（0.01～0.06の一様分布）
SD_Acc_HF <- round(runif(n_total, min = 0.01, max = 0.06), 3)
# 個人ごとのランダム切片（安静時HRの違い）
subject_intercepts <- rnorm(n_subjects, mean = 70, sd = 5)
intercept_ID <- subject_intercepts[ID]
# 個人ごとのランダム傾き（SD_Acc_HFへの感度の違い）
subject_slopes <- rnorm(n_subjects, mean = 230, sd = 20)
slope_ID <- subject_slopes[ID]
# 曜日ごとの固定効果（金曜日 = 0）
weekday_effects_fixed <- c(Mon = -2, Tue = -1.5, Wed = 0.5, Thu = -1, Fri = 0)
# 個人ごとの曜日効果（ランダムにばらつかせる）
weekday_random_effects <- matrix(rnorm(n_subjects * 4, mean = 0, sd = 1),
nrow = n_subjects, ncol = 4)
colnames(weekday_random_effects) <- c("Mon", "Tue", "Wed", "Thu")
# 曜日ごとの個人効果（IDによって異なる）
weekday_effect_total <- numeric(n_total)
for (i in 1:n_total) {
day <- Weekday[i]
subj <- ID[i]
fixed <- if (day == "Fri") 0 else weekday_effects_fixed[day]
random <- if (day == "Fri") 0 else weekday_random_effects[subj, day]
weekday_effect_total[i] <- fixed + random
}
# 残差（観測誤差など）
error <- rnorm(n_total, mean = 0, sd = 2)
# 心拍数の計算
HR <- round(intercept_ID + slope_ID * SD_Acc_HF + weekday_effect_total + error, 1)
# データフレーム作成
df <- data.frame(
ID = factor(ID),
Measurement = Measurement,
Weekday = factor(Weekday, levels = c("Mon", "Tue", "Wed", "Thu", "Fri")),
SD_Acc_HF = SD_Acc_HF,
HR = HR
)
# データ表示（上位）
head(df, 10)
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
[blog:g:11696248318754550903:banner][blog:g:11696248318754550877:banner] [blog:g:11696248318755496791:banner]
