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

次元の海で溺れる

Rとデータ解析と統計手法たちとわたし

ggplot2で信頼区間の描画をして遊びたい 練習編

面白いデータセットを探し続けているものの
なかなか見つからないので歌詞でも収集して形態素に切ってしまおうプロジェクトが
水面下で進行中です。

ネットワークでも描こうか、、とも思ったけど卒論で死ぬ思いした残像がまだ、、、


###

SappoRo.R前にggplot2で遊ぼうシリーズ第2弾、
信頼図を描こう。です。

テーマが地味。

テレビとかに出てくるお医者さんやら何かしらの専門家の人が、
「信頼区間95%」をうまーく表現丸めて喋ってるのを見ると尊敬します。

『○○にならない可能性が5%、つまりかなり高い確率で〇〇になるってことですね~』

絶妙なライン攻めてるかんじ。


ggplot2初心者なので、
相変わらず教科書はこれです。

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

いつもお世話になっております。

ではさっそく。

#---信頼区間を図示してみたい--------------------

#使用するデータセットについて
#gcookbookパッケージよりclimateデータセットをお借りします

install.packages("gcookbook")
library(gcookbook)
library(ggplot2)


> #データ構造

> head(climate)
    Source Year Anomaly1y Anomaly5y Anomaly10y Unc10y
1 Berkeley 1800        NA        NA     -0.435  0.505
2 Berkeley 1801        NA        NA     -0.453  0.493
3 Berkeley 1802        NA        NA     -0.460  0.486
4 Berkeley 1803        NA        NA     -0.493  0.489
5 Berkeley 1804        NA        NA     -0.536  0.483
6 Berkeley 1805        NA        NA     -0.541  0.475

> str(climate)
'data.frame':	499 obs. of  6 variables:
 $ Source    : chr  "Berkeley" "Berkeley" "Berkeley" "Berkeley" ...
 $ Year      : num  1800 1801 1802 1803 1804 ...
 $ Anomaly1y : num  NA NA NA NA NA NA NA NA NA NA ...
 $ Anomaly5y : num  NA NA NA NA NA NA NA NA NA NA ...
 $ Anomaly10y: num  -0.435 -0.453 -0.46 -0.493 -0.536 -0.541 -0.59 -0.695 -0.763 -0.818 ...
 $ Unc10y    : num  0.505 0.493 0.486 0.489 0.483 0.475 0.468 0.461 0.453 0.451 ...
> summary(climate)
    Source               Year        Anomaly1y          Anomaly5y         Anomaly10y      
 Length:499         Min.   :1800   Min.   :-0.60070   Min.   :-0.4995   Min.   :-1.01500  
 Class :character   1st Qu.:1884   1st Qu.:-0.21629   1st Qu.:-0.1053   1st Qu.:-0.28350  
 Mode  :character   Median :1926   Median :-0.02797   Median :-0.0042   Median :-0.07328  
                    Mean   :1923   Mean   : 0.01277   Mean   : 0.0555   Mean   :-0.07869  
                    3rd Qu.:1968   3rd Qu.: 0.15155   3rd Qu.: 0.1620   3rd Qu.: 0.05065  
                    Max.   :2011   Max.   : 0.96354   Max.   : 0.8953   Max.   : 0.88400  
                                   NA's   :207        NA's   :373       NA's   :20        
     Unc10y      
 Min.   :0.0110  
 1st Qu.:0.0430  
 Median :0.1040  
 Mean   :0.1452  
 3rd Qu.:0.2220  
 Max.   :0.5050  
 NA's   :294     

#変数は6つ。
#Anomaly10yは1950年~1980年までの平均気温からの偏差の10年移動平均
#Unc10yは95%信頼区間
#Sourceは出典というか、調査機関って感じですかね。NASAとか。
#Anomaly1yと5yは1年移動平均と5年移動平均ぽいですがSourceによって欠損値がバラバラなので今回は不使用

#使うデータ部分を切り出してきて格納
#Berkely調査の平均気温偏差10年移動平均、その信頼区間を使用

clim=subset(climate,climate$Source=="Berkeley",select=c("Year","Anomaly10y","Unc10y"))

> clim
    Year Anomaly10y Unc10y
1   1800     -0.435  0.505
2   1801     -0.453  0.493
3   1802     -0.460  0.486
4   1803     -0.493  0.489
5   1804     -0.536  0.483
6   1805     -0.541  0.475
7   1806     -0.590  0.468
8   1807     -0.695  0.461
9   1808     -0.763  0.453
10  1809     -0.818  0.451
11  1810     -0.842  0.437
12  1811     -0.912  0.406
13  1812     -0.963  0.395
14  1813     -1.015  0.389
15  1814     -0.949  0.369
16  1815     -0.935  0.368
17  1816     -0.918  0.361
18  1817     -0.802  0.330
19  1818     -0.743  0.349
20  1819     -0.703  0.382
21  1820     -0.625  0.395
22  1821     -0.527  0.363
23  1822     -0.392  0.339
24  1823     -0.332  0.351
25  1824     -0.321  0.343
26  1825     -0.278  0.344
27  1826     -0.270  0.295
28  1827     -0.378  0.284
29  1828     -0.398  0.254
30  1829     -0.388  0.243
31  1830     -0.429  0.247
32  1831     -0.506  0.247
33  1832     -0.625  0.257
34  1833     -0.689  0.246
35  1834     -0.667  0.224
36  1835     -0.689  0.220
37  1836     -0.727  0.221
38  1837     -0.662  0.218
39  1838     -0.612  0.222
40  1839     -0.613  0.256
41  1840     -0.638  0.274
42  1841     -0.542  0.267
43  1842     -0.472  0.282
44  1843     -0.407  0.269
45  1844     -0.402  0.285
46  1845     -0.390  0.285
47  1846     -0.351  0.284
48  1847     -0.346  0.280
49  1848     -0.362  0.273
50  1849     -0.357  0.260
51  1850     -0.290  0.244
52  1851     -0.321  0.239
53  1852     -0.372  0.228
54  1853     -0.363  0.226
55  1854     -0.351  0.228
56  1855     -0.348  0.216
57  1856     -0.374  0.203
58  1857     -0.418  0.187
59  1858     -0.420  0.186
60  1859     -0.427  0.178
61  1860     -0.473  0.172
62  1861     -0.458  0.154
63  1862     -0.412  0.149
64  1863     -0.393  0.145
65  1864     -0.373  0.141
66  1865     -0.344  0.148
67  1866     -0.317  0.147
68  1867     -0.282  0.145
69  1868     -0.269  0.145
70  1869     -0.241  0.149
71  1870     -0.233  0.149
72  1871     -0.255  0.146
73  1872     -0.271  0.142
74  1873     -0.247  0.140
75  1874     -0.248  0.132
76  1875     -0.271  0.128
77  1876     -0.252  0.124
78  1877     -0.235  0.120
79  1878     -0.257  0.115
80  1879     -0.297  0.119
81  1880     -0.321  0.124
82  1881     -0.314  0.123
83  1882     -0.325  0.117
84  1883     -0.376  0.108
85  1884     -0.373  0.111
86  1885     -0.368  0.110
87  1886     -0.386  0.109
88  1887     -0.397  0.107
89  1888     -0.394  0.108
90  1889     -0.364  0.104
91  1890     -0.342  0.102
92  1891     -0.323  0.099
93  1892     -0.298  0.100
94  1893     -0.273  0.107
95  1894     -0.286  0.107
96  1895     -0.251  0.106
97  1896     -0.201  0.104
98  1897     -0.176  0.104
99  1898     -0.153  0.104
100 1899     -0.174  0.106
101 1900     -0.171  0.108
102 1901     -0.162  0.109
103 1902     -0.177  0.108
104 1903     -0.199  0.104
105 1904     -0.223  0.105
106 1905     -0.241  0.107
107 1906     -0.294  0.106
108 1907     -0.312  0.105
109 1908     -0.328  0.103
110 1909     -0.281  0.101
111 1910     -0.247  0.099
112 1911     -0.243  0.097
113 1912     -0.257  0.100
114 1913     -0.268  0.100
115 1914     -0.257  0.097
116 1915     -0.249  0.095
117 1916     -0.214  0.096
118 1917     -0.201  0.096
119 1918     -0.176  0.096
120 1919     -0.182  0.097
121 1920     -0.193  0.097
122 1921     -0.167  0.098
123 1922     -0.128  0.096
124 1923     -0.075  0.097
125 1924     -0.064  0.098
126 1925     -0.065  0.100
127 1926     -0.050  0.100
128 1927     -0.020  0.099
129 1928     -0.018  0.099
130 1929     -0.026  0.100
131 1930     -0.014  0.101
132 1931     -0.047  0.098
133 1932     -0.035  0.096
134 1933     -0.017  0.093
135 1934      0.020  0.092
136 1935      0.053  0.089
137 1936      0.063  0.085
138 1937      0.048  0.081
139 1938      0.073  0.079
140 1939      0.113  0.076
141 1940      0.113  0.072
142 1941      0.134  0.071
143 1942      0.134  0.069
144 1943      0.127  0.070
145 1944      0.111  0.068
146 1945      0.072  0.066
147 1946      0.035  0.066
148 1947      0.042  0.064
149 1948      0.045  0.063
150 1949      0.013  0.062
151 1950      0.010  0.058
152 1951     -0.017  0.054
153 1952     -0.040  0.047
154 1953     -0.040  0.043
155 1954     -0.032  0.038
156 1955     -0.022  0.035
157 1956      0.012  0.031
158 1957      0.007  0.028
159 1958      0.002  0.027
160 1959      0.002  0.026
161 1960     -0.019  0.026
162 1961     -0.001  0.021
163 1962      0.017  0.018
164 1963      0.004  0.016
165 1964     -0.028  0.018
166 1965     -0.006  0.017
167 1966     -0.024  0.017
168 1967     -0.041  0.019
169 1968     -0.025  0.020
170 1969     -0.019  0.024
171 1970      0.010  0.026
172 1971      0.007  0.022
173 1972      0.015  0.015
174 1973      0.028  0.012
175 1974      0.049  0.014
176 1975      0.068  0.012
177 1976      0.128  0.011
178 1977      0.158  0.012
179 1978      0.167  0.013
180 1979      0.193  0.012
181 1980      0.186  0.016
182 1981      0.217  0.016
183 1982      0.235  0.014
184 1983      0.270  0.014
185 1984      0.318  0.014
186 1985      0.344  0.013
187 1986      0.352  0.012
188 1987      0.380  0.011
189 1988      0.370  0.013
190 1989      0.366  0.017
191 1990      0.433  0.019
192 1991      0.467  0.018
193 1992      0.496  0.017
194 1993      0.526  0.019
195 1994      0.554  0.020
196 1995      0.563  0.019
197 1996      0.565  0.022
198 1997      0.618  0.022
199 1998      0.680  0.023
200 1999      0.734  0.025
201 2000      0.748  0.026
202 2001      0.793  0.027
203 2002      0.856  0.028
204 2003      0.869  0.028
205 2004      0.884  0.029

これで準備はOKなはず。
コメントにあれこれ書いてますが、
Climateデータセットの一部を切り出してきて使用したってことです。

#いざ描画

#信頼区間は
#Anomaly10yにUnc10yを加算したものを上限(ymax)
#Anomaly10yからUnc10yを減算したものを下限(ymin)として設定
#網掛けの透過度は80%(alpha=0.2)

ggplot(clim,aes(x=Year,y=Anomaly10y))+
  geom_ribbon(aes(ymin=Anomaly10y-Unc10y,ymax=Anomaly10y+Unc10y),alpha=0.2)+
  geom_line()

f:id:WAFkw:20150523161942p:plain

できたー
イメージどおりー

geom_ribbon(網掛け)の上に、geom_lineを重ねてるので
こういう順番でのコード記述になる。ふーん。

てことは

順番を逆にしたら網掛けが上に来てしまって見えなくなる・・・?

ためしてみた

#あえて網掛け部分と実線を逆に描画してみる
#しかも網掛けの透過度を0%にしてみる(alpha=1)
ggplot(clim,aes(x=Year,y=Anomaly10y))+
  geom_line()+
 geom_ribbon(aes(ymin=Anomaly10y-Unc10y,ymax=Anomaly10y+Unc10y),alpha=1)

#なるほど見えない

f:id:WAFkw:20150523162208p:plain


コメントにも書いてますが、なるほど見えない。

一回透過80%でやってみたら実線が透けて見えちゃってたので、
最終的には黒く塗りつぶしてみた。


せっかくなので街で良く見る点線パターンでの描画もしてみる

#点線バージョン

ggplot(clim,aes(x=Year,y=Anomaly10y))+
  geom_line(aes(y=Anomaly10y-Unc10y),colour="blue",linetype="dashed")+
  geom_line(aes(y=Anomaly10y+Unc10y),colour="red",linetype="dashed")+
  geom_line()

f:id:WAFkw:20150523162422p:plain

いろいろ遊んではみたけど、
linetypeはdottedよりdashedの方が見やすいと思う派でした。
colourも上限下限で色変えた方が自分のコードの反映状況が分かりやすいかなーとおもって。

またいろいろ遊んでみよう。


整体のじかんだー