読者です 読者をやめる 読者になる 読者になる

Mukku John Blog

取り組んでいること を つらつら と

ggplot2を使って、散布図を作る-4

Rグラフィックスクックブック 散布図

R グラフィックス クックブック 13回目

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. モデルを作る
  2. 予測するデータセットを用意する
  3. モデルから予測値を求める
  4. 求めた予測値を散布図にする
  5. 散布図にモデルをフィットさせた線をプロットする
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

f:id:MukkuJohn:20160831212731p:plain
この辺は、今まで変わりませんね。

5.散布図にモデルをフィットさせた線をプロットする

上で作成した散布図に、モデルの線をプロットします。

線をプロットするには、geom_line()を使います。
geom_line()は、折れ線をプロットするために使いました。
mukkujohn.hatenablog.com

sp + geom_line(data=predicted,size = 1)

f:id:MukkuJohn:20160831213915p:plain
うん、まぁ、こんな感じですね。

色々なモデルで試すことを簡略化するために、上の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)

f:id:MukkuJohn:20160831215600p:plain
線形回帰モデルの線赤の線ですが、最初にプロットした線とは異なります。
今回は、説明変数に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)

f:id:MukkuJohn:20160831221000p:plain



さて、ここからは、グループごとに、モデル構築~プロットを行います。
使っているデータセットには、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)

f:id:MukkuJohn:20160831223653p:plain

ここで気になるのが、それぞれの線の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)

f:id:MukkuJohn:20160831224049p:plain
x軸の範囲が揃いました。


余談ですが、グループごとにプロットする際に
グループをファセットで表現しようと、こんなコードを書いてみました。

上のグラフ <- ggplot(heightweight,aes(x=ageYear,y=heightIn,colour=sex)) +
  geom_point() +
  geom_line(data=predvals)

オブジェクト名が、日本語でもOK!!

そして、上のグラフオブジェクトですが、普通に使えました。

上のグラフ +
  facet_grid(. ~ sex, scales = "free_x", space = "free_x")

f:id:MukkuJohn:20160831223700p:plain
普段は使わないです。絶対。

でも、こうやって、ブログにしていると、
色々と省けるところもあるので、記載内容が前後している時は使おうかなぁ。

次回は、散布図上に、データを説明するための係数の表示や、
ラグ付け/ラベル付けを行っていきます。