はじめに

線形回帰分析では残差が正規分布に従うことを前提としています。しかし、データが常に正規分布に従うという保証はありません。例えば、目的変数が「正解・不正解」「はい・いいえ」のような2値データがこれに該当します。線形回帰分析を拡張させた一般化線形回帰分析では、このように正規分布に従っていないデータに対応することが可能です。目的変数を$Y$、説明変数を$X$、係数を$A$、誤差を$E$とすると、線形回帰分析は$Y=AX+E$と表すことができ、一般化線形回帰分析では非線形関数を$g(μ)=AX$に変換することで線形モデルとして扱うことができます。なお、$μ$は目的変数の平均、$g(μ)$はリンク関数と呼ばれています。

線形回帰分析はRでlm関数、一般化線形回帰分析ではglm関数を使うことで分析することができます。関数glmの基本的な書式は以下のように表されます。

glm(formula, family, data)

formulaは「目的変数 ~ 説明変数」のようにlm関数と同じように使うことができます。lm関数と異なるのはfamilyの部分です。familyにはリンク関数を指定します。デフォルトではgaussian (正規分布) が指定されています。

分布	リンク $g(μ)$	リンク関数
正規 (gaussian)	$μ$	link = “identity”
二項 (binomial)	log[$μ/(1-μ)$]	link = “logit”
ポアソン (poisson)	log($μ$)	link = “log”
ガンマ (Gamma)	1/$μ$	link = “inverse”
逆正規 (Inverse.gaussian)	1/$μ^2$	link = “1/mu^2”

しかし、glm関数で対応できない複雑なデータも存在します。以下では平滑化回帰と加法モデルについて見ていきましょう。

平滑化回帰

まずはデータをシミュレーションで作成してみましょう。

set.seed(123)
x1 <- seq(-10, 10, 0.1)
y1 <- 50 * sin(x1) + x1^2 + 10 * rnorm(length(x1), 0, 1)
plot(x1, y1)

作成されるデータは以下の通りです。

Untitled

それでは非線形回帰である多項式によるモデルを作成し、当てはまりを確認してみます。今回は2次、5次、8次の多項式を作成しました。

library(tidyverse)
lm_p2 <- lm(y1 ~ poly(x1, 2)) %>%　fitted()
lm_p5 <- lm(y1 ~ poly(x1, 5)) %>%　fitted()
lm_p8 <- lm(y1 ~ poly(x1, 8)) %>%　fitted()

続いてはこれらを可視化してみます。

# 元のデータをプロット
plot_data <- tibble(x1, y1)

# 各多項式フィットのデータを追加
plot_data <- plot_data %>%
mutate(
　fit_p2 = lm_p2,
　fit_p5 = lm_p5,
　fit_p8 = lm_p8
 )

# 可視化
ggplot(plot_data, aes(x = x1)) +
 geom_point(aes(y = y1), color = "black") +
 geom_line(aes(y = fit_p2, colour = "2次多項式")) +
 geom_line(aes(y = fit_p5, colour = "5次多項式")) +
 geom_line(aes(y = fit_p8, colour = "8次多項式")) +
 theme_minimal() +
 labs(y = "y1", title = "多項式によるフィットの比較") +
      scale_colour_manual(
      name = "モデルタイプ",
      values = c("2次多項式" = "blue", "5次多項式" = "red", "8次多項式" = "green")) +
 theme(text = element_text(family= "HiraKakuProN-W3"))

Untitled

プロットを見ると、2次式よりは8次式の方が当てはまりが良さそうですが、いずれもデータにフィットしているとは言い難いです。

そこで複雑に変化するデータの回帰方法である平滑化回帰を見ていきましょう。平滑化回帰では局所的な重み付き最小二乗回帰を用いて、データに滑らかなカーブをフィットさせる手法です。

# 平滑化回帰分析
loess_fit <- loess(y1 ~ x1, span = 0.3) # spanは平滑度を調整
summary(loess_fit)

結果は以下の通りです。

> summary(loess_fit)
Call:
loess(formula = y1 ~ x1, span = 0.3)

Number of Observations: 201
Equivalent Number of Parameters: 9.92
Residual Standard Error: 10.53
Trace of smoother matrix: 10.96 (exact)

Control settings:
span : 0.3
degree : 2
family : gaussian
surface : interpolate cell = 0.2
normalize: TRUE
parametric: FALSE
drop.square: FALSE

それぞれの解釈は以下の通りです。

等価パラメータ数 (Equivalent Number of Parameters): モデルの複雑さを示す指標です。この値が大きいほどモデルは複雑です。
残差標準誤差 (Residual Standard Error): モデルの予測が実際のデータからどの程度離れているかを示す指標です。この値が小さいほど、モデルの予測が実際のデータに近いことを意味します。
平滑行列のトレース (Trace of smoother matrix): これは平滑化の程度を示す指標です。
制御設定 (Control settings): モデルをフィットさせる際に使用されたさまざまな設定を示しています。例えば、span = 0.3 は平滑化の程度を、degree = 2 は使用された多項式の次数を示しています。

得られた結果を可視化してみましょう。