R Language
線形モデル（回帰）

構文

• lm（数式、データ、サブセット、重み、na.action、メソッド= "qr"、モデル= TRUE、x = FALSE、y = FALSE、qr = TRUE、singular.ok = TRUE、コントラスト= NULL、オフセット、.. 。）

パラメーター

mtcarsデータセットの線形回帰

ビルトインのmtcars データフレームには、重量、燃費（ガロンあたりのマイル数）、速度などの32台の車に関する情報が含まれています（詳細については、 help(mtcars)使用しhelp(mtcars) ）。

燃費（ mpg ）と重量（ wt ）の関係に興味がある場合、次のようにこれらの変数をプロットすることがあります。

plot(mpg ~ wt, data = mtcars, col=2)

プロットは（線形の）関係を示しています。線形回帰を実行して線形モデルの係数を求める場合は、 lm関数を使用します。

fit <- lm(mpg ~ wt, data = mtcars)

ここでの~は「説明する」を意味するので、 mpg ~ wtという式はwtで説明されるmpgを予測していることを意味します。出力を表示する最も有益な方法は次のとおりです。

summary(fit)

それは出力を与える：

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5432 -2.3647 -0.1252  1.4096  6.8727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,    Adjusted R-squared:  0.7446 
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

これは以下の情報を提供します：

• mpgの最適予測を示唆する各係数（ wtとy切片）の推定勾配は37.2851 + (-5.3445) * wt
• 切片と体重はおそらく偶然によるものではないことを示唆する各係数のp値
• R ^ 2や調整されたR ^ 2などのフィットの全体的な推定値は、 mpgの変動のどれがモデルによって説明されるかを示します

予測されたmpgを表示する最初のプロットにラインを追加することができます：

abline(fit,col=3,lwd=2)

そのプロットに方程式を追加することもできます。まず、 coef係数をcoefます。その後、 paste0を使用して係数を適切な変数と+/-でpaste0 、方程式を構築します。最後に、 mtextを使用してプロットに追加しmtext 。

bs <- round(coef(fit), 3) 
lmlab <- paste0("mpg = ", bs[1],
             ifelse(sign(bs[2])==1, " + ", " - "), abs(bs[2]), " wt ")
mtext(lmlab, 3, line=-2) 

結果は次のとおりです。

回帰をプロットする（ベース）

mtcars例を続けると、ここでは、出版に適した線形回帰のプロットを作成する簡単な方法があります。

最初に線形モデルをフィットさせ、

fit <- lm(mpg ~ wt, data = mtcars)

次に、関心のある2つの変数をプロットし、定義ドメイン内に回帰直線を追加します。

plot(mtcars$wt,mtcars$mpg,pch=18, xlab = 'wt',ylab = 'mpg')
lines(c(min(mtcars$wt),max(mtcars$wt)),
as.numeric(predict(fit, data.frame(wt=c(min(mtcars$wt),max(mtcars$wt))))))

ほぼそこに！最後のステップは、プロット、回帰方程式、rsquare、および相関係数を追加することです。これは、 vector関数を使用して行われます。

rp = vector('expression',3)
rp[1] = substitute(expression(italic(y) == MYOTHERVALUE3 + MYOTHERVALUE4 %*% x), 
          list(MYOTHERVALUE3 = format(fit$coefficients[1], digits = 2),
                        MYOTHERVALUE4 = format(fit$coefficients[2], digits = 2)))[2]
rp[2] = substitute(expression(italic(R)^2 == MYVALUE), 
             list(MYVALUE = format(summary(fit)$adj.r.squared,dig=3)))[2]
rp[3] = substitute(expression(Pearson-R == MYOTHERVALUE2), 
             list(MYOTHERVALUE2 = format(cor(mtcars$wt,mtcars$mpg), digits = 2)))[2]

legend("topright", legend = rp, bty = 'n')

ベクトル関数を適応させることによって、RMSEなどの他のパラメータを追加できます。 10要素の伝説が必要だと想像してみてください。ベクトル定義は次のようになります。

rp = vector('expression',10)

r[1] .... r[10]を定義する必要がありますr[10]

出力は次のとおりです。

加重

時には、モデルが他のものよりもいくつかのデータ点や例に大きな重みを付けることを望みます。これは、モデルの学習中に入力データの重みを指定することで可能になります。例に対して不均一な重みを使用するシナリオには、一般的に2種類あります。

• 分析重み：異なる観測の異なる精度レベルを反映します。たとえば、各観測値が地理的エリアからの平均結果であるデータを分析する場合、分析加重は推定分散の逆数に比例します。データの平均値を扱う場合に、観測数を指定して比重量を指定すると便利です。 ソース

• サンプリングウェイト（Inverse Probability Weights - IPW）：データが収集された集団とは異なる集団に対して標準化された統計を計算するための統計的手法。異なる標本集団および標的推定の集団（標的集団）を有する研究デザインは、適用において一般的である。欠損値を持つデータを扱う場合に便利です。 ソース

lm()関数は解析的な重み付けを行います。サンプリングウェイトの場合、 surveyパッケージを使用して調査設計オブジェクトを作成し、 svyglm()を実行します。デフォルトでは、 surveyパッケージはサンプリング加重を使用します。 （注： lm()およびsvyglm()家族とgaussian()どちらも、加重最小二乗を最小化することにより、係数を解くため、すべてが、同じ点推定値が生成されます彼らはどのように計算されるか標準誤差が異なります。）

テストデータ

data <- structure(list(lexptot = c(9.1595012302023, 9.86330744180814, 
8.92372556833205, 8.58202430280175, 10.1133857229336), progvillm = c(1L, 
1L, 1L, 1L, 0L), sexhead = c(1L, 1L, 0L, 1L, 1L), agehead = c(79L, 
43L, 52L, 48L, 35L), weight = c(1.04273509979248, 1.01139605045319, 
1.01139605045319, 1.01139605045319, 0.76305216550827)), .Names = c("lexptot", 
"progvillm", "sexhead", "agehead", "weight"), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -5L))

アナリティックウェイト

lm.analytic <- lm(lexptot ~ progvillm + sexhead + agehead, 
                            data = data, weight = weight)
summary(lm.analytic)

出力

Call:
lm(formula = lexptot ~ progvillm + sexhead + agehead, data = data, 
    weights = weight)

Weighted Residuals:
         1          2          3          4          5 
 9.249e-02  5.823e-01  0.000e+00 -6.762e-01 -1.527e-16

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.016054   1.744293   5.742    0.110
progvillm   -0.781204   1.344974  -0.581    0.665
sexhead      0.306742   1.040625   0.295    0.818
agehead     -0.005983   0.032024  -0.187    0.882

Residual standard error: 0.8971 on 1 degrees of freedom
Multiple R-squared:  0.467, Adjusted R-squared:  -1.132 
F-statistic: 0.2921 on 3 and 1 DF,  p-value: 0.8386

サンプリングウェイト（IPW）

library(survey)
data$X <- 1:nrow(data)             # Create unique id
        
# Build survey design object with unique id, ipw, and data.frame
des1 <- svydesign(id = ~X,  weights = ~weight, data = data)
        
# Run glm with survey design object
prog.lm <- svyglm(lexptot ~ progvillm + sexhead + agehead, design=des1)

出力

Call:
svyglm(formula = lexptot ~ progvillm + sexhead + agehead, design = des1)

Survey design:
svydesign(id = ~X, weights = ~weight, data = data)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 10.016054   0.183942  54.452   0.0117 *
progvillm   -0.781204   0.640372  -1.220   0.4371  
sexhead      0.306742   0.397089   0.772   0.5813  
agehead     -0.005983   0.014747  -0.406   0.7546  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2078647)

Number of Fisher Scoring iterations: 2

多項式回帰による非線形性のチェック

場合によっては、線形回帰で作業するときに、データの非線形性をチェックする必要があります。これを行う1つの方法は、多項式モデルに適合させ、線形モデルよりもデータに適合するかどうかをチェックすることです。変数関係が本質的に多項式であると考えられるため、2次以上のモデルに適合することを示す理論的なような他の理由がある。

mtcarsデータセットに二次モデルを適合させましょう。線形モデルについては、mtcarsデータセットの線形回帰を参照してください。

まず変数mpg （Miles / gallon）、 disp （Displacement（cu.in.））、 wt （Weight（1000 lbs））の散布図を作成します。 mpgとdispの関係は非線形に見えます。

plot(mtcars[,c("mpg","disp","wt")])

線形近似はdispが重要でないことを示す。

fit0 = lm(mpg ~ wt+disp, mtcars)
summary(fit0)

# Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 34.96055    2.16454  16.151 4.91e-16 ***
#wt          -3.35082    1.16413  -2.878  0.00743 ** 
#disp        -0.01773    0.00919  -1.929  0.06362 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.917 on 29 degrees of freedom
#Multiple R-squared:  0.7809,    Adjusted R-squared:  0.7658

次に、2次モデルの結果を得るために、 I(disp^2)を追加しました。新しいモデルはR^2を見るとよりよく見え、すべての変数が重要です。

fit1 = lm(mpg ~ wt+disp+I(disp^2), mtcars)
summary(fit1)

# Coefficients:
#              Estimate Std. Error t value Pr(>|t|)    
#(Intercept) 41.4019837  2.4266906  17.061  2.5e-16 ***
#wt          -3.4179165  0.9545642  -3.581 0.001278 ** 
#disp        -0.0823950  0.0182460  -4.516 0.000104 ***
#I(disp^2)    0.0001277  0.0000328   3.892 0.000561 ***
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 2.391 on 28 degrees of freedom
#Multiple R-squared:  0.8578,    Adjusted R-squared:  0.8426 

3つの変数があるので、当てはめられたモデルは、

mpg = 41.4020-3.4179*wt-0.0824*disp+0.0001277*disp^2

多項式回帰を指定するもう1つの方法は、パラメータraw=TRUE polyを使用する方法です。そうでない場合は直交多項式を考慮します（詳細はhelp(ploy)を参照）。同じ結果が得られます：

summary(lm(mpg ~ wt+poly(disp, 2, raw=TRUE),mtcars))

最後に、推定された表面のプロットを表示する必要がある場合はどうなりますか？さて、 R 3Dプロットを作るための多くのオプションがあります。ここでは、 p3dパッケージのFit3dを使用します。

library(p3d)
Init3d(family="serif", cex = 1)
Plot3d(mpg ~ disp+wt, mtcars)
Axes3d()
Fit3d(fit1)

品質評価

回帰モデルを作成した後、結果を確認し、モデルが適切であり、手元のデータとうまく適合しているかどうかを判断することが重要です。これは、残差プロットおよび他の診断プロットを調べることによって行うことができます。

# fit the model
fit <- lm(mpg ~ wt, data = mtcars)
# 
par(mfrow=c(2,1))
# plot model object
plot(fit, which =1:2)

これらのプロットは、モデルを構築する際に行われた2つの前提をチェックします。

1. 予測変数（この場合はmpg ）の期待値は、予測変数の線形結合（この場合はwt ）によって与えられます。この推定値は偏っていると予想しています。したがって、残差は、プレディクタのすべての値に対して平均値を中心にしなければなりません。この場合、残差は両端で正であり、中間で負である傾向があり、変数間に非線形の関係があることがわかります。
2. 実際の予測された変数は、その推定値の周りに正規分布していること。したがって、残差は正規分布する必要があります。正規分布データの場合、通常のQQプロット内の点は対角線上にあるか、または対角線に近くなければなりません。ここにはいくらかのスキューがあります。

'predict'関数を使用する

モデルが構築されると、新しいデータでテストする主な機能がpredictされます。この例では、 mtcarsビルトインデータセットを使用して、変位に対してガロンあたりのマイルを後退させます。

my_mdl <- lm(mpg ~ disp, data=mtcars)
my_mdl

Call:
lm(formula = mpg ~ disp, data = mtcars)

Coefficients:
(Intercept)         disp  
   29.59985     -0.04122

私がディスプレースメント付きの新しいデータソースを持っていれば、ガロンあたりの推定マイル数を知ることができました。

set.seed(1234)
newdata <- sample(mtcars$disp, 5)
newdata
[1] 258.0  71.1  75.7 145.0 400.0

newdf <- data.frame(disp=newdata)
predict(my_mdl, newdf)
       1        2        3        4        5 
18.96635 26.66946 26.47987 23.62366 13.11381

このプロセスの最も重要な部分は、元のデータと同じ列名を持つ新しいデータフレームを作成することです。この場合、元のデータにはdispというラベルの列があり、同じ名前の新しいデータを必ず呼び出すことになりました。

あぶない

いくつかの一般的な落とし穴を見てみましょう：

1. 新しいオブジェクトにdata.frameを使用しないでください：

   predict(my_mdl, newdata)
   Error in eval(predvars, data, env) : 
      numeric 'envir' arg not of length one
   

2. 新しいデータフレームで同じ名前を使用しない：

   newdf2 <- data.frame(newdata)
   predict(my_mdl, newdf2)
   Error in eval(expr, envir, enclos) : object 'disp' not found
   

正確さ

予測の精度を確認するには、新しいデータの実際のy値が必要です。この例では、 newdfは 'mpg'と 'disp'の列が必要です。

newdf <- data.frame(mpg=mtcars$mpg[1:10], disp=mtcars$disp[1:10])
#     mpg  disp
# 1  21.0 160.0
# 2  21.0 160.0
# 3  22.8 108.0
# 4  21.4 258.0
# 5  18.7 360.0
# 6  18.1 225.0
# 7  14.3 360.0
# 8  24.4 146.7
# 9  22.8 140.8
# 10 19.2 167.6

p <- predict(my_mdl, newdf)

#root mean square error
sqrt(mean((p - newdf$mpg)^2, na.rm=TRUE))
[1] 2.325148