ggplot2を使って、散布図を作る-4
R グラフィックス クックブック １３回目
ggplot2パッケージを利用して、散布図を作成していきます？
散布図にカテゴリしてますが、モデル作って、フィットする線をプロットします。
あるデータに、フィットするモデルの線を重ねて表示する場合は、
stat_smooth()を使うのが手っ取り早いです。
mukkujohn.hatenablog.com
使うデータ
今までも使ってきた、こちらのデータを使います。
> str(heightweight) 'data.frame': 236 obs. of 5 variables: $ sex : Factor w/ 2 levels "f","m": 1 1 1 1 1 1 1 1 1 1 ... $ ageYear : num 11.9 12.9 12.8 13.4 15.9 ... $ ageMonth: int 143 155 153 161 191 171 185 142 160 140 ... $ heightIn: num 56.3 62.3 63.3 59 62.5 62.5 59 56.5 62 53.8 ... $ weightLb: num 85 105 108 92 112 ...
モデルを作る
こんな手順で進めていきます。
- モデルを作る
- 予測するデータセットを用意する
- モデルから予測値を求める
- 求めた予測値を散布図にする
- 散布図にモデルをフィットさせた線をプロットする
1.モデルを作る
lm()で、線形モデルが作成できます。
#heightIn ~ ageYear + I(ageYear^2)：モデル式 # 目的変数:heightIn # 説明変数:ageYear 、 I(ageYear^2) model <- lm(heightIn ~ ageYear + I(ageYear^2), heightweight)
2.予測するデータセットを用意する
最小値と最大値を決めて、その範囲内で指定した個数のデータを用意します。
#最小値を取る xmin <- min(heightweight$ageYear) #最大値を取る xmax <- max(heightweight$ageYear) #最小値～最大値まで、等間隔のデータを100個用意する predicted <- data.frame(ageYear=seq(xmin, xmax, length.out = 100))
predictedオブジェクトは、こんなデータセットです。
> str(predicted) 'data.frame': 100 obs. of 1 variable: $ ageYear: num 11.6 11.6 11.7 11.8 11.8 ...
3.モデルから予測値を求める
モデルから予測値を求めるためには、predict()を使います。
#モデルから予測された値を、heightIn列に入れる predicted$heightIn <- predict(model, predicted)
上記のコードを実行した後の、
predictedオブジェクトは、このような形式になります。
> str(predicted) 'data.frame': 100 obs. of 2 variables: $ ageYear : num 11.6 11.6 11.7 11.8 11.8 ... $ heightIn: num 56.8 57 57.2 57.3 57.5 ...
4.求めた予測値を散布図にする
散布図にするのは、今まで扱ってきた通り
geom_point()を使います。
#後で使うため、オブジェクトspに格納しておきます。 sp <- ggplot(heightweight, aes(x=ageYear,y=heightIn)) + geom_point(colour = "grey40")
図を表示します。
sp
この辺は、今まで変わりませんね。
5.散布図にモデルをフィットさせた線をプロットする
上で作成した散布図に、モデルの線をプロットします。
線をプロットするには、geom_line()を使います。
geom_line()は、折れ線をプロットするために使いました。
mukkujohn.hatenablog.com
sp + geom_line(data=predicted,size = 1)
うん、まぁ、こんな感じですね。
色々なモデルで試すことを簡略化するために、上の2.と3.を関数にします。
#関数の引数は、モデル、説明変数、目的変数、範囲、個数です #範囲、個数はデフォルト値がついています。 predictvals <- function(model, xvar, yvar, xrange=NULL, samples=100, ...){ #説明変数をモデルから取る方法が、モデルによって違います #範囲が指定されていない場合は、モデルから値を取って範囲に使います if(is.null(xrange)){ if(any(class(model) %in% c("lm","glm"))) xrange <- range(model$model[[xvar]]) else if (any(class(model) %in% "loess")) xrange <- range(model$x) } #予測するためのデータを作る。2.で行っていることです newdata <- data.frame(x = seq(xrange[1], xrange[2], length.out = samples)) #列名を設定します names(newdata) <- xvar #予測値を求めます newdata[[yvar]] <- predict(model, newdata=newdata, ...) newdata }
上の関数がある前提で、1.～3.の手順がこうなります。
#線形回帰モデルを作ります modlinear <- lm(heightIn ~ ageYear, heightweight) #loessモデルを作ります modloess <- loess(heightIn ~ ageYear, heightweight)
#線形回帰モデルで予測されたデータセットを作ります lm_predicted <- predictvals(modlinear,"ageYear","heightIn") #loessモデルで予測されたデータセットを作ります loess_predicted <- predictvals(modloess,"ageYear","heightIn")
散布図に、線形回帰モデルの線と、loessモデルの線をプロットします
sp + geom_line(data=lm_predicted, colour="red", size=.8) + geom_line(data=loess_predicted, colour="blue",size=.8)
線形回帰モデルの線（赤の線）ですが、最初にプロットした線とは異なります。
今回は、説明変数にageYearしか使っていませんね。
最初にモデルを作成した時は、ageYearの2乗値も説明変数としています。
ロジスティック回帰についても確認してみます。
ロジスティック回帰 - Wikipedia
この記事と同様に、データを用意します。
mukkujohn.hatenablog.com
library(MASS) > str(biopsy) 'data.frame': 699 obs. of 11 variables: $ ID : chr "1000025" "1002945" "1015425" "1016277" ... $ V1 : int 5 5 3 6 4 8 1 2 2 4 ... $ V2 : int 1 4 1 8 1 10 1 1 1 2 ... $ V3 : int 1 4 1 8 1 10 1 2 1 1 ... $ V4 : int 1 5 1 1 3 8 1 1 1 1 ... $ V5 : int 2 7 2 3 2 7 2 2 2 2 ... $ V6 : int 1 10 2 4 1 10 10 1 1 1 ... $ V7 : int 3 3 3 3 3 9 3 3 1 2 ... $ V8 : int 1 2 1 7 1 7 1 1 1 1 ... $ V9 : int 1 1 1 1 1 1 1 1 5 1 ... $ class: Factor w/ 2 levels "benign","malignant": 1 1 1 1 1 2 1 1 1 1 ...
V1列と
class列をファクタから変えたclassn列の値を使います。
b <- biopsy b$classn[b$class=="benign"] <- 0 b$classn[b$class=="malignant"] <- 1 > str(b[,c("V1","classn")]) 'data.frame': 699 obs. of 2 variables: $ V1 : int 5 5 3 6 4 8 1 2 2 4 ... $ classn: num 0 0 0 0 0 1 0 0 0 0 ...
モデルを作って、予測値を求めます。
#classn列値を、V1値から求めるモデルを作ります。 fitlogistic <- glm(classn ~ V1, b, family=binomial()) #predictvals関数を使って、予測値を求めます #ロジスティック回帰は、predict関数にrespoinseを指定する必要があります glm_predicted <- predictvals(fitlogistic,"V1","classn",type="response")
過去記事と同様に、同じポイントにデータが多くオーバープロットになっています。
ジッターをかけて点をプロットし、予測値の線をプロットします。
ggplot(b, aes(x=V1,y=classn)) + geom_point(position=position_jitter(width = .3,height=.08),alpha=0.4,shape=21,size=1.5)+ geom_line(data=glm_predicted,colour="#1177FF",size=1)
さて、ここからは、グループごとに、モデル構築～プロットを行います。
使っているデータセットには、sex列があるので、この値でグループ化します。
data.frameに対する操作を行う必要がありますので、plyrパッケージを使います。
library(plyr)
さらに、グループごとにモデルを作りやすくするために、関数を作っておきます。
#データセットが固定で、線形回帰モデル固定だけど… make_model <- function(data){ lm(heightIn ~ ageYear, data) }
heightweightデータセットを、sex値で分割して、
make_model関数を適用した結果のリストを作ります。
models <- dlply(heightweight, "sex", .fun= make_model) > class(models) [1] "list"
リストから、predictvals関数を通して、予測値のdata.frameにします
predvals <- ldply(models, .fun=predictvals, xvar="ageYear", yvar="heightIn")
上と同様にプロットします。
ggplot(heightweight,aes(x=ageYear,y=heightIn,colour=sex)) + geom_point() + geom_line(data=predvals)
ここで気になるのが、それぞれの線のx軸範囲が異なることですね。
予測値のdata.frameを作成する際に、x軸範囲を指定をしないで、
predictvals関数を使ったため、グループごとに範囲が決定しています。
x軸範囲を揃えるためには、predictvals関数を使う際に、
指定すれば良いだけです。
predvals <- ldply(models, .fun=predictvals, xvar="ageYear",yvar="heightIn", xrange=range(heightweight$ageYear)) ggplot(heightweight,aes(x=ageYear,y=heightIn,colour=sex)) + geom_point() + geom_line(data=predvals)
x軸の範囲が揃いました。
余談ですが、グループごとにプロットする際に
グループをファセットで表現しようと、こんなコードを書いてみました。
上のグラフ <- ggplot(heightweight,aes(x=ageYear,y=heightIn,colour=sex)) + geom_point() + geom_line(data=predvals)
オブジェクト名が、日本語でもＯＫ！！
そして、上のグラフオブジェクトですが、普通に使えました。
上のグラフ + facet_grid(. ~ sex, scales = "free_x", space = "free_x")
普段は使わないです。絶対。
でも、こうやって、ブログにしていると、
色々と省けるところもあるので、記載内容が前後している時は使おうかなぁ。
次回は、散布図上に、データを説明するための係数の表示や、
ラグ付け／ラベル付けを行っていきます。
