Tex, python, illusrator, VPSの学生ノート

latotex-blog

R

【R】線形回帰の推定区間(上限値・下限値)に直線を引く

投稿日:2020年10月9日 更新日:







ggplot2を使えば95%信頼区間を自動的に塗りつぶしてくれますが、
どうしても直線で表現したかったので、ggplot2を使う場合と使わない場合で試しました。

使うデータは、ある異次元スポーツ漫画のキャラクターの身長体重です。
データの詳細は以下の記事に載せていますので、ご覧ください。

スポンサーリンク

該当コード

詳しい説明は後の項でしますので、取り敢えずコードを先に載せます。

まずは、ggplot2を使う場合です。
勿体ないですが、推定区間を自動的に塗りつぶす機能をオフにします。

#ライブラリの宣言
library(tidyverse) #データ抽出・グラフ作図など
library(latex2exp) #TeX記法を有効化

#csvの読み込み
input <- read_csv('C:/Users/user/Estimation of Intervals/input.csv')

#回帰直線の決定係数、p値を求める
model_input <- lm(weight ~ height, data = input)
summary(model_input)

#データフレームを昇順にし、X軸、y軸に対応する値に整理、
x <- input %>% arrange(height) %>% select(height)
y <- input %>% arrange(weight) %>% select(weight)


#データフレームから行列に変換
x <- as.matrix(x)
y <- as.matrix(y)

#データの数
n <- nrow(x)

#回帰式と推定区間の上限値・下限値をもとめる
Sxx <- sum((x - mean(x))^2)  # 説明変数xの偏差平方和
Sxy <- sum((x - mean(x)) * (y - mean(y)))  # 説明変数xと目的変数yの積和
a <- Sxy/Sxx  # 回帰直線の傾きの推定量
b <- mean(y) - mean(x) * a  # 回帰直線の切片の推定量
estimated_y <- a * x + b  # 推定されたy

residuals <- y - estimated_y  # 残差
sigma <- sqrt(sum((residuals - mean(residuals))^2)/(n - 2))  # 残差標準偏差

significant_level <- 0.05  # 有意水準
t_value <- qt(1 - significant_level/2, n - 2)  # 自由度n-2、有意水準0.05のt値

ci <- t_value * sqrt(1/n + ((x - mean(x))^2)/Sxx) * sigma #区間推定の式
lower <- estimated_y - ci #上限値
upper <- estimated_y + ci #下限値


#ggplotを使う場合
fig_scatter = ggplot(data = input, aes(x = x, y = y), size = 2) + #X軸に身長、Y軸に体重をとる
              stat_smooth(method = lm, se = F, col='blue1') + #回帰直線を引くが、区間はse=Falseで消す
              geom_point(col='green4') +
              geom_line(aes(x = x, y = upper), col='blue1', lty='dashed') + #上限値の直線を引く
              geom_line(aes(x = x, y = lower), col='blue1', lty='dashed') + #下限値の直線を引く
              ylim(30,120) +
              theme_classic() +

              labs(x = TeX('Height ($\\mathrm{cm}$)'),
                  y = TeX('Weight ($\\mathrm{kg}$)')) + #xlabelの設定、TeX()でTeX表記が可能でバックスラッシュ(\)は2つ必要
              annotate('text', x = 190,   y = 40,
                      label = TeX('$R^2$ = $0.689$ (p < $0.001$, n = 50)'), size = 5) + #決定係数、p値の表示(少数点は$$で囲むことで半角のピリオドになる)
              theme(
                    axis.title = element_text(size = 15),
                    rect = element_rect(fill = "transparent"))

#保存
ggsave('C:/Users/user/Desktop/Estimation of Intervals/fig_intervals_ggplot.png', fig_scatter)

ggplot2_lineplot_interval_95

次に、ggplot2を使わない場合です。
やっていることはggplot2を使う場合と同じで、作図方法が異なるだけです。

#描図デバイスを開く(保存)
png('C:/Users/user/Desktop/Estimation of Intervals/test_fig_intervals.png', height=960, width=960, res=144)

plot(x, y,
    xlab= TeX('Height ($\\mathrm{cm}$)'),  ylab=TeX('Weight ($\\mathrm{kg}$)'),
    ylim=c(30,120), col='green4', pch = c(16)) #ggplotと違いfillは使えない。markerstyleの設定(pch)で変更する
lines(x, estimated_y, col='red', lwd = 2) #回帰直線を引く
lines(x, upper, col='red', lty='dashed', lwd = 2) #上限値の直線を引く
lines(x, lower, col='red', lty='dashed', lwd = 2) #下限値の直線を引く
text(x = 190, y = 40, TeX('$R^2$ = $0.689$ (p < $0.001$, n = 50)')) #決定係数、p値の表示(少数点は$$で囲むことで半角のピリオドになる)

#描図デバイスをとじる(保存)
dev.off()

R_lineplot_interval_95

補足:回帰係数の求め方

計算方法は以下のサイト様のコードを参照しています。

簡潔に説明すると、コード上の「a = Sxy/Sxx」の意味は
回帰係数の式(Ⅰ)に、相関係数の式(Ⅱ)を代入するということです。

(1)    \begin{flalign*}   &  \text{(I)} \quad \alpha  =  \mathit{corrcoef}(X, Y) × \frac{SD(Y)}{SD(X)} & \notag \\   &  \text{(I\hspace{-.1em}I)} \quad\mathit{corrcoef}(X, Y)  =  cov(X, Y) × \frac{1}{SD(X)×SD(Y)} \notag & \\ \notag \\   & \text{Substitute Equation (I) for Equation (I\hspace{-.1em}I).} \notag &\\   & \qquad \alpha  =  cov(X, Y) × \frac{1}{SD(X)×SD(Y)} × \frac{SD(Y)}{SD(X)} &  \notag\\   & \quad \qquad =  cov(X, Y) × \frac{1}{SD(X)^2} \notag \end{flalign*}

よって、回帰直線の傾きの推定量には説明変数xと目的変数yの積和と説明変数xの偏差平方和が必要であることがわかります。

注意!プロットするデータは昇順にする必要がある。

プロットするデータによっては、複数線が引かれた図になる場合があります。
これは、Rまたはggplot2がデータが並んでいる順に線を引いているからです。
よって、データを昇順に並べる必要があります。arrange()が地味に重要です。

ちなみに並べないとこんな感じになります・・・・・

fail_lineplot_interval_95


-R
-, ,

Copyright© latotex-blog , 2020 All Rights Reserved Powered by STINGER.