スタートページJavaScript他言語R言語

lm(), predict() 回帰分析

下記の青線の部分をGoogle Colaboratryの「コード」部分にコピーアンドペースト(ペーストは Cntl+V)して実行すれば、右図の画像が表示されます。


lm() の概要

lm関数は、線形モデル(LinearModel)の分野での代表的関数です。基本的には説明変数、被説明変数が量的データで、主に回帰分析や重回帰分析に利用されます。
y = a + bx + ε (ε は誤差項)の形式でεを小さくすることを目的にしますが、εは正規分布に従うことが仮定されています。

εが任意の分布であるとき、説明変数に量的データを含む場合、被説明変数が2値の場合などでは、glm 関数(Generalized linear model、一般化線形モデル)を用います。

lm() の書式と引数

関数の一般形
    lm(formula,        // 必須。モデルの形式
       data,           // 必須。分析のためのデータ
       method = "qr",
        ...)

formula
    y~x             y = a + bx + ε (ε は誤差項)   
    y~x-1           y = bx + ε   (定数項を除外)
    y~x1+x2         y = a + b1x1 + b2x2 + ε
    y~1+x+I(x^2)    y = b0 + b1x1 + b2x2 + ε
    y~x1+x2+x1*x2   y = a + b1x1 + b2x2 + b3x1x2 + ε

data ベクトルで与える
    x1 <- c(  0, 10, 25, 30, 35, 50, 60)
    y  <- c( 65, 60, 50, 60, 65, 70, 80)

lm() からの戻り値

一般形式
    結果 <- lm(formula, data) 
       summary(結果)

主な戻り値   詳細は後述の実例説明参照
    summary(結果)          回帰分析の結果一覧
    print(結果)            簡略版を表示.
    coefficients(kekka)    回帰係数 (a, b1, b2 )
    deviance(結果)        重みつけられた残差平方和.
    plot(結果)             残差,当てはめ値などの 4 種類のプロットを生成.
    predict(結果, newdata) newdata を回帰式に当てはめた結果

単回帰

(x,y) の組から
  散布図
  相関係数 r
  回帰直線 y = a + bx の算出と表示
  信頼区間曲線の表示
を行います。
 

(x,y)の相関係数= 0.6167
回帰式: y= 56.078 + 0.2735 *x
決定係数 r2= 0.3803, 相関係数 r= 0.6167
p値= 0.140152

x <- c(  0, 10, 25, 30, 35, 50, 60)
y <- c( 65, 60, 50, 60, 65, 70, 80)

xmin <-  0;  xmax <- 60
ymin <- 30;  ymax <- 90     # 区間図がスケートアウトしないように広く設定

# ===== 散布図の作成
plot(x,y,                                  # plot のオプション指定
     xlim=c(xmin, xmax), ylim=c(ymin,ymax),  # 表示範囲
     pch=16, col="red", cex=3,               # 点の図形
     xaxt="n", yaxt="n", xlab="", ylab="")   # 後述で設定するので非表示にする

axis(side=1, at=seq(xmin, xmax, 10))         # X軸目盛り
axis(side=2, at=seq(ymin, ymax, 10))         # Y軸目盛り
abline(h=seq(ymin, ymax, 5),lty="dotted")    # 水平補助線
abline(v=seq(xmin, xmax, 5),lty="dotted")    # 垂直補助線
title(main="regression",                     # タイトル
      xlab="x-axis", ylab="y-axis")          # 軸ラベル

# ===== X-Y相関係数
r <- cor(x,y)
paste("(x,y)の相関係数=", r)

# ===== 回帰分析
kekka <- lm(y~x)                             # y = a + bx の関係
summary(kekka)
# 回帰直線の係数
係数 <- coefficients(kekka)
paste("回帰式:",係数[1], "+", 係数[2], "*x")

# ========== 決定係数、p値
r2 <- summary(kekka)$r.squared
r <- sqrt(r2)
para <- summary(kekka)$fstatistic
p <- pf(para[1], para[2], para[3], lower.tail=F)
paste("決定係数 r2=",r2, ", 相関係数 r=",r, ", p値=", p)

# ===== 信頼区間
newx <- data.frame(x = seq(xmin, xmax, 1))
              # 補間のためのX軸データの設定(データフレーム形式にする必要がある)
newyc <- predict(kekka, newdata=newx,  # newx による信頼区間上・下限値の計算
               interval='confidence',  # 信頼区間
               level=0.95)             # 95% 信頼
newyp <- predict(kekka, newdata=newx,  # newx による信頼区間上・下限値の計算
               interval='prediction',  # 予測区間
               level=0.95)             # 95% 信頼
# 信頼区間の表示
lines(newx$x, newyc[,"fit"], lwd=4, col='red')    # 回帰式
lines(newx$x, newyc[,"upr"], lwd=2, col="blue")   # 信頼区間上限の曲線
lines(newx$x, newyc[,"lwr"], lwd=2, col="blue")   # 信頼区間下限の曲線
lines(newx$x, newyp[,"upr"], lwd=2, col="green")  # 予測上限の曲線
lines(newx$x, newyp[,"lwr"], lwd=2, col="green")  # 予測下限の曲線

cor(x,y)

点 (x,y) の相関関係です。回帰式とは無関係。

lm(y~x)

回帰分析により回帰直線の係数を求める関数です。
    # 回帰式の構成
    #   y~x         y= a0 + a1*x            xはベクトル
    #   y~x-1       y= a1*x                 切片0(比例)
    #   y~x+I(x^2)  y= a0 + a1*x + a2*x^2   高次式
    #   y~x1+x2     y= a0 + a1*x1 + a2*x2   多項式

kekka <- lm(y~x);  summary(kekka)  とすれば、次の情報が得られます
    # 
    # Call:
    # lm(formula = y ~ x)
    # 
    # Residuals: 残差(回帰式との差)
    #      Min      1Q    Median     3Q     Max 
    #   8.9218   1.1860 -12.9178  -4.2857  -0.6536   0.2426   7.5067 
    # 
    # Coefficients:
    #             係数      標準偏差   t値     確率  有意差水準
    #             Estimate Std. Error t value Pr(>|t|)    
    # (Intercept)  56.0782     5.5835  10.043 0.000167 ***
    # x             0.2736     0.1562   1.752 0.140152    
    # ---           └ y = 56.0782 + 0.2736x
    # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    # 
    # Residual standard error: 8.038 on 5 degrees of freedom
    # Multiple R-squared:  0.3804,	Adjusted R-squared:  0.2565
    #                       └ r2(決定係数)
    # F-statistic:  3.07 on 1 and 5 DF,  p-value: 0.1402 ← P値
    #                └ 自由度が1と5でのF値    └ <0.05 なら回帰式とよく合致、>0.05 なら合致しない

coefficients(kekka)  とすれば、上表から係数をの取出せます
    # (Intercept)      x 
    # 56.0781671   0.2735849 

summary(kekka)$r.squared 上のR-squared(決定係数)を取り出す
    # 回帰式と点との関係を示します。
    # 決定係数が0.38とは、62%の残差があるということで、あまりよい近似ではありません。

summary(kekka)$fstatistic 上のF-statisticを取り出す
    # p値:帰無仮説「母集団が回帰式に合致しない」 の発生確率が14%  
    # 0.05以下ならば、仮説は棄却され「母集団が回帰式に95%合致する」と結論できる。
    # ここでは 0.1402>0.05 なので「棄却できない」(どちらともいえない)ことになる。

predict(kekka, newdata=newx)

kekka <- lm(y~x) で得た回帰式に、x=newdata の値を代入した結果を得る関数です。

    level=0.95                # 有意水準(変更できます)
  interval="confidence"     # 信頼区間(平均予測値の周りの不確実性を反映)
             "prediction"     # 予測区間(単一の値の周りの不確実性を反映)

  信頼区間(図の青線)
   回帰式の信頼区間。母集団の母数(標本から測定できない)に対して「どの範囲にあると推定できるか」を示す。
   測定・解析を同一方法で100回繰り返したら、100回のうち95回は回帰直線はこの範囲を通ると考えられる。
  予測区間(図の緑線)
   母集団を仮定した上で、将来観察されるであろう標本値に対して「どの範囲にあると予測されるか」を示す。
      予測区間が残差も考慮にいれているので、区間幅は広くなる。

    全体をまとめると次のようになります。
  #                信頼区間     予測区間
  #  x   y    fit    lwr.x  upr.x   lwr.y  upr.y
  #          回帰式  下限   上限    下限   上限
  #  0  65   56.07   41.72  70.43   30.91  81.23
  # 10  60   58.81   47.61  70.01   35.31  82.31
  # 25  50   62.91   54.85  70.98   40.73  85.09
  # 30  60   64.28   56.19  71.83   41.91  86.10
  # 35  65   65.65   57.58  73.71   43.47  87.83
  # 50  70   69.75   58.55  80.95   46.25  93.26
  # 60  80   72.49   58.14  86.84   47.33  97.65

二次式 y = a + bx + cx^2 の当てはめ


x-y 相関係数 = 0.6167
回帰式: y= 64.426 -0.6870*x +0.0160*x^2
決定係数 r2= 0.8415, 相関係数 r= 0.9173
p値= 0.0251

kekka <- lm(y~x+I(x^2)) と変えただけです。


# ========== 係数算出
kekka <- lm(y~x+I(x^2))  # 説明変数が二次式のときの指定
summary(kekka) # 抜粋
    # Coefficients:
    #              Estimate     回帰式
    # (Intercept) 64.426186   a    y = 64.426 - 0.687x + 0.016x^2
    # x           -0.687009   b
    # I(x^2)       0.016010   c
    #                        ┌ 決定係数
    # Multiple R-squared:  0.8415,	Adjusted R-squared:  0.7623 
    # F-statistic: 10.62 on 2 and 4 DF,  p-value: 0.02512
    #                                              └ <0.05 95%で有意差あり

predict(kekka)    # xに対応する回帰式の値
  #                信頼区間     予測区間
  #  x   y    fit    lwr.x  upr.x   lwr.y  upr.y
  #          回帰式  下限   上限    下限   上限
  #  0  65   64.42   53.33  75.51   47.62  81.22
  # 10  60   59.15   52.31  66.00   44.79  73.51
  # 30  60   58.22   51.36  65.08   43.85  72.58
  # 50  70   70.10   63.25  76.94   55.74  84.45
  # 60  80   80.84   69.75  91.93   64.04  97.64


重回帰分析

y = a + b*x1 + c*x2 のように、説明変数 x が複数の場合です。lm(y~x1+x2) になります。
X軸が定まらないのでグラフは描けません。その代わりに、x1,x2,y, 回帰式の値、区間の値 を表形式で出力します。 そのため、入力データをデータフレームの形式で与えるのが便利です。

# ========== 入力データ
df <- data.frame(          # データフレームの形式で与えると便利
         x1 = c(  0, 10, 25, 30, 35, 50, 60),
         x2 = c( 30, 20, 10, 15, 20, 20, 25),
         y  = c( 65, 60, 50, 60, 65, 70, 80))
print("■■入力データ■■")
df
# ========== 回帰分析
kekka <- lm(y~x1+x2)
print("■■lm() simmary■■")
summary(kekka)

print("■■回帰式■■")
係数 <- coefficients(kekka)
paste("回帰式:", 係数[1], "+", 係数[2], "*x1", 係数[3], "*x2")

# ========== 決定係数、p値
r2 <- summary(kekka)$r.squared
r <- sqrt(r2)
para <- summary(kekka)$fstatistic
p <- pf(para[1], para[2], para[3], lower.tail=F)
print("■■決定係数等■■")
paste("決定係数 r2=",r2, ", 相関係数 r=",r, ", p値=", p)

# ===== 信頼区間・予測区間の計算
newx <- data.frame(x = seq(1,7, 1))
newyc <- predict(kekka, newdata = newx,
                interval = 'confidence', level = 0.95)     # 新X点に対応する信頼区間95%限界点
newyp <- predict(kekka, newdata = newx,
                interval = 'prediction', level = 0.95)     # 新X点に対応する予測区間95%限界点

# ========== 表の作成
df$回帰式   <- newyc[, "fit"]            # 回帰式を計算した結果 "fit" を列に追加
df$信頼下限 <- newyc[, "lwr"]
df$信頼上限 <- newyc[, "upr"]
df$予測下限 <- newyp[, "lwr"]
df$予測上限 <- newyp[, "upr"]
print("■■結果の表■■")
df

結果の表示

決定係数 = 0.979、p値 = 0.0004 など、かなりよい近似になりました。

"■■入力データ■■"
   x1   x2   y     ←列名 これに回帰式が信頼区間などの列が追加される
    0   30   65
   10   20   60
   25   10   50
   30   15   60
   35   20   65
   50   20   70
   60   25   80

■■lm() simmary■■"
Call:
lm(formula = y ~ x1 + x2)

Residuals:
      1       2       3       4       5       6       7 
-1.0673  2.0361 -1.4409  1.3465 -0.8662 -0.6076  0.5993 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 32.27422    2.47284  13.051 0.000199 ***
x1           0.31609    0.03194   9.895 0.000585 ***
x2           1.12644    0.10400  10.831 0.000412 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.632 on 4 degrees of freedom
Multiple R-squared:  0.9796,	Adjusted R-squared:  0.9694 
F-statistic:  95.9 on 2 and 4 DF,  p-value: 0.0004174

■■回帰式■■"
'回帰式: 32.2742200328407 + 0.316091954022988 *x1 1.1264367816092 *x2'

■■決定係数等■■"
'決定係数 r2= 0.979570146433632 , 相関係数 r= 0.989732361011618 , p値= 0.000417378916743219'

■■結果の表■■"
    #   x1   x2   y    回帰式   信頼下限   信頼上限   予測下限   予測上限
    #    0   30   65   66.067    62.010     70.124     59.985     72.149
    #   10   20   60   57.963    55.498     60.429     52.805     63.122
    #   25   10   50   51.440    48.008     54.873     45.756     57.125
    #   30   15   60   58.653    56.413     60.893     53.599     63.707
    #   35   20   65   65.866    64.097     67.635     61.002     70.730
    #   50   20   70   70.607    68.141     73.073     65.449     75.765
    #   60   25   80   79.400    75.789     83.011     73.606     85.194