線形回帰の95%CI:ポリゴンver.

線形モデルの95%CIをグレイのシェードで描きたいときは、polygon()を使います。ポアソンモデルを例に、グレイシェイドの95%CIを描きます(下図)。

SPONSOR LINK

#1. 統計モデルの構築
model <- glm(y~x, family=poisson)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- exp(pred$fit) #予測値の平均値
up <- exp(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- exp(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogなので、predict()の返値をlogの逆関数exp()で変換する

#4. polygon()のx,yの値を設定する
X <- c(new$x, rev(new$x))
Y <- c(low, rev(up))
# rev()はベクターの並び順を逆転させる

#5. プロット
plot(y~x, type="n", ylim=c(min(Y,y),max(Y,y)),
     ann=F, axes=F)
box(bty="l")
polygon(X,Y,border = F, col=rgb(0,0,0,alpha=0.3))
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish abundance", side = 2, line = 3)

polygon()を後のほうに描いてしまうと、先に描いたポイントやラインも塗りつぶしてしまうので、先に描いておきます。また、プロットのy軸の範囲をylim=c(min(Y,y),max(Y,y))とすることで、データポイントおよび95%CIの最大・最小値をカバーできる描画領域を指定しています。

※polygon()補足
polygon()は、グラフ上のXY座標を指定し、その座標に囲まれた範囲を塗りつぶします。例えば、次のようなデータフレームがあるとします。

#ポリゴンの座標
X <- c(1,2,2,1); Y <- c(1,1,2,2)
data.frame(X,Y)
  X Y
1 1 1
2 2 1
3 2 2
4 1 2

これは、座標を(X, Y)と表記したとき、データ1が(1, 1)を示し、データ4が(1, 2)を示す、という意味なります。これをpolygon()に書かせてみます(右図)。polygon()では、第一引数でX座標を指定し、第二引数でY座標を指定します。

#ポリゴンを描画するプロット領域
plot(0, type="n", xlim=c(0,3),ylim=c(0,3),
     xlab="X", ylab="Y")
polygon(X,Y,col="black")
# colは塗りつぶしの色指定;borderは縁どりの有無を指定

線形回帰の95%CI:二項分布

二項分布は、「複数回の試行のうち、何回成功したか」のようなデータ(上限値が分かっている非負の離散値)に対して適用されることが多い確率分布です(ベルヌーイ分布の一般化)。生態学のデータでは、種の在/不在データに適用されることが多いと思います。二項分布のパラメーターは、ある事象Xが生じる確率pであり、pは0-1の範囲に収まる必要があります。そのため、リンク関数にはlogit関数(logit(p) = log(p/(1-p)))が用いられます。二項分布モデルの場合、predict()の返値をlogitスケールになっていますので、予測値をプロットするためには逆ロジット関数(inverse logit: ilogit(x) = exp(x)/(1+exp(x))で変換する必要があります。その他、注意点などについてはポアソン分布の場合と同様です。サンプルデータは記事末尾にあります。

SPONSOR LINK

#1. 統計モデルの構築
model <- glm(y~x, family=binomial)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

#4. inverse logit functionを作る
ilogit <- function(x){ exp(x)/(1+exp(x)) }
# Rのデフォルト関数にilogit関数はない

#5. predict()をつかって予測値を導く
y_pred <- ilogit(pred$fit) #予測値の平均値
up <- ilogit(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- ilogit(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogitなので、predict()の返値をlogitの逆関数ilogit()で変換する

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Probability of fish present", side = 2, line = 3)
ilogit <- function(x){ exp(x)/(1+exp(x)) }
N <- 50
a <- 0.4; b <- 1.5
x <- rnorm(N,0,1)

p <- y <- NULL
for(i in 1:N){
  p[i] <- ilogit(a + b*x[i])
  y[i] <- rbinom(1,1,p[i])
}

線形回帰の95%CI:ポアソン分布

ポアソン分布は、非負の離散値(カウントデータなど)に対して適用されることの多い確率分布です。ポアソン分布を仮定した線形モデルでは、母数である平均ラムダが負の値をとれないため、リンク関数にlogを使います。そのため、predict()の返値は対数スケールの値になっており、予測値曲線をプロットする際にはもとのスケールに戻す必要があります。記事末尾のサンプルデータを使い、ポアソンモデルの予測曲線をプロットします(右図)。

#1. 統計モデルの構築
model <- glm(y~x, family=poisson)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- exp(pred$fit) #予測値の平均値
up <- exp(pred$fit + 1.96*pred$se.fit) #Upper 95%CI:平均予測値に1.96*SEを加える
low <- exp(pred$fit - 1.96*pred$se.fit) #Lower 95%CI
# リンク関数がlogなので、predict()の返値をlogの逆関数exp()で変換する
# SEの演算は指数変換を施す前に行う
# ダメな例:exp(pred$fit) + exp(1.96*pred$se.fit)

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish abundance", side = 2, line = 3)

基本的なスクリプトは正規分布の場合と変わりませんが、辺値を変換するステップが必須となります。また、予測値の95%CIを出す際に、exp()で変換後に演算すると負の値が現れたりしますので注意が必要です。

※注 predict(model, new, se.fit=T, type=”response”)とすると、実測値スケールで予測値を返してくれますが、このSEをもとに四則演算を行うと上記の「ダメな例」と同様のことをすることになります。

N <- 50; a <- 1.4; b <- 0.7
x <- rnorm(N,0,1)

lambda <- y <- NULL
for(i in 1:N){
  lambda[i] <- exp(a + b*x[i])
  y[i] <- rpois(1,lambda[i])
}

線形回帰の95%CI:正規分布

lmやglmの予測値をプロットし、それらの信頼区間を合わせて表示したいケースは多いかと思います。正規分布のモデルであればpredict()の引数を調整することで簡単に描けますが、それ以外の確率分布となったときに案外描けなかったりします。以下のスクリプトでは正規分布を例に示していますが、この手順を覚えればほかの確率分布のモデルにも適用できます(多少の変更は必要)。プロットの基本的な調整でつくったサンプルデータを使ってスクリプトを書きます。以下のスクリプトから右図が描けます。

#1. 統計モデルの構築
model <- lm(y~x)

#2. モデルの予測値を示すための説明変数のデータフレームを作る
new <- data.frame(x = seq(min(x),max(x),length=100))
#予測用データフレーム内の説明変数の名称はモデルの変数名と一致する必要あり
#seq(min, max, length=サンプル数)となっており、最小値から最大値までの範囲で等間隔のベクターを生成

#3. predict()をつかって予測値を導く
pred <- predict(model,new,se.fit=T)

y_pred <- pred$fit #予測値の平均値
up <- pred$fit + 1.96*pred$se.fit #Upper 95%CI:平均予測値に1.96*SEを加える
low <- pred$fit - 1.96*pred$se.fit #Lower 95%CI

#4. プロット
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
lines(new$x, y_pred) # 平均予測値のラインを描く
lines(new$x, up, lty=2); lines(new$x, low, lty=2)# 95%CIを描く
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

おまけ:matplot()を使った方法

Y <- cbind(up,y_pred,low)
matplot(new$x,Y, type="l", lty=c(2,1,2), col=rgb(0,0,0), ann=F, axes=F)
box(bty="l")
points(x,y,cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3))
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

なお、複数の説明変数がある重回帰モデルの場合、注目する説明変数以外は平均値(もしくは中央値)に固定したデータセットを作る必要があります(そうしない限り、ほかの変数の影響を受けた予測値になる)。以下にサンプルを示します。

model <- lm(y~x1+x2+x3)
new_x1 <- seq(min(x1),max(x1),length=100)
new <- data.frame(x1 = new_x1, x2 = median(x2), x3 = median(x3))

半透明のシンボル

半透明のプロットを描く場合、シンボルの中身の色を調整できるシンボルを指定します。ここでは、シンボルの引数pchは21とし(円の中身の色指定ができるシンボル)、bgで中身の色を指定します。pch=21の場合、colはシンボルの縁どりの色を指定します。

rgb()はred,green,blueの三原色の割合を指定することで色調整する関数です。最初の三つの引数がそれぞれRGBに対応しています。すべて0とすると黒、1とすると白になります。alphaが透明度を指定する引数で、0(完全な透明)から1(完全な不透明)の値をとります。

#points()内の引数の値を調整
plot(y~x, type="n", ann=F, axes=F)
box(bty="l")
points(x, y, cex=1.5, pch=21, col=NA, bg=rgb(0,0,0,alpha=0.3)
axis(1); axis(2, las=2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)

 

プロットの調整

下のplotを使い、見た目の調整方法についてまとめます。私個人の意見では、一度まっさらなプロット領域を用意し、好みの設定でデータや軸を追加していくのがのちのち微調整しやすいスクリプトになると思います。

調整前
#1. バックグラウンドとなるプロットを作る
plot(y~x, type="n", ann=F, axes=F)
#type="n"でプロット領域に何も書かない
#ann=F で両軸ラベルの消去
#axes=F で両軸の消去

#2. 外枠の追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
#bty="l"でL字型の外枠

#3. データポイントの追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex=1.5)
#cexでポイントサイズの調整:2とするとデフォルトの二倍の大きさ

#4. 軸の追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex = 1.5)
axis(1); axis(2, las = 2)
#axisの第一引数は軸の位置を指定:1=下横軸、2=左縦軸、3=上横軸、4=右横軸
#las=2とすることで軸の値を横表記

#5. 軸ラベルの追加
plot(y~x, type="n", ann=F, axes=F)
box(bty = "l")
points(x, y, cex=1.5)
axis(1); axis(2, las = 2)
mtext("Standardized lake area", side = 1, line = 3)
mtext("Fish density", side = 2, line = 3)
#sideはラベル位置の指定:1=下、2=左、3=上、4=右
#lineはラベルの軸からの距離(行数)

※ステップ1~6まで継ぎ足しで書いてます

調整後

今回使ったサンプルデータの生成

N <- 50; a <- 3.4; b <- 0.2
x <- rnorm(N,0,1)
sigma <- 0.3

mu <- y <- NULL
for(i in 1:N){
  mu[i] <- a + b*x[i]
  y[i] <- rnorm(1,mu[i],sigma)
}