R Language
線形モデル(回帰)
サーチ…
構文
- lm(数式、データ、サブセット、重み、na.action、メソッド= "qr"、モデル= TRUE、x = FALSE、y = FALSE、qr = TRUE、singular.ok = TRUE、コントラスト= NULL、オフセット、.. 。)
パラメーター
パラメータ | 意味 |
---|---|
式 | ウィルキンソン=ロジャース表記の数式。 response ~ ... where ... は、環境内の変数またはdata 引数で指定されたdata フレームに対応する項が含まれています |
データ | 応答変数と予測変数を含むデータフレーム |
サブセット | 使用される観測のサブセットを指定するベクトルは、 data 内の変数の観点から論理的なステートメントとして表現することができdata |
重み | 分析重量(上記の重量部を参照) |
na.action | 欠損( NA )値の処理方法:参照: ?na.action |
方法 | フィッティングの実行方法。選択肢は"qr" または"model.frame" (後者は、モデルをフィットせずにモデルフレームを返します。これはmodel=TRUE を指定するのと同じです) |
モデル | フィットしたオブジェクトにモデルフレームを格納するかどうか |
バツ | 適合したオブジェクトにモデル行列を格納するかどうか |
y | 適合したオブジェクトにモデル応答を格納するかどうか |
qr | フィットしたオブジェクトにQR分解を格納するかどうか |
singular.ok | 特異フィットを可能にするかどうか、共線予測子を持つモデル(この場合、係数のサブセットは自動的にNA に設定されます |
対照 | モデル内の特定の要因に使用されるコントラストのリスト。 ?model.matrix.default contrasts.arg 引数を参照してください。コントラストは、 options() ( contrasts 引数を参照)またはあるcontrast 属性を割り当てることによって設定することもできoptions() ?contrasts 参照) |
オフセット | モデル内の先験的な既知のコンポーネントを指定するために使用されます。数式の一部として指定することもできます。 See ?model.offset |
... | 下位フィッティング関数( lm.fit() またはlm.wfit() )に渡される追加の引数 |
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つの前提をチェックします。
- 予測変数(この場合は
mpg
)の期待値は、予測変数の線形結合(この場合はwt
)によって与えられます。この推定値は偏っていると予想しています。したがって、残差は、プレディクタのすべての値に対して平均値を中心にしなければなりません。この場合、残差は両端で正であり、中間で負である傾向があり、変数間に非線形の関係があることがわかります。 - 実際の予測された変数は、その推定値の周りに正規分布していること。したがって、残差は正規分布する必要があります。正規分布データの場合、通常の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
というラベルの列があり、同じ名前の新しいデータを必ず呼び出すことになりました。
あぶない
いくつかの一般的な落とし穴を見てみましょう:
新しいオブジェクトにdata.frameを使用しないでください:
predict(my_mdl, newdata) Error in eval(predvars, data, env) : numeric 'envir' arg not of length one
新しいデータフレームで同じ名前を使用しない:
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