Rmarkdown
fumita
2019/12/16
https://www.slideshare.net/KojiKosugi/ss-50740386?next_slideshow=1
library(readr)## Warning: package 'readr' was built under R version 3.6.2income = read_csv("income.csv")## Parsed with column specification:
## cols(
## id = col_character(),
## sex = col_character(),
## age = col_double(),
## height = col_double(),
## weight = col_double(),
## income = col_number(),
## generation = col_character()
## )head(income)## # A tibble: 6 x 7
## id sex age height weight income generation
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 AU male 70 160. 58.3 201 elder
## 2 AY female 70 156. 44 487 elder
## 3 AB male 69 173. 75.7 424 elder
## 4 AM male 67 166. 69.3 1735 elder
## 5 M male 66 171. 76.5 929 elder
## 6 CM female 66 164. 67.3 397 elderdim(income)## [1] 100 7summary(income)## id sex age height
## Length:100 Length:100 Min. :20.00 Min. :148.0
## Class :character Class :character 1st Qu.:36.00 1st Qu.:158.1
## Mode :character Mode :character Median :45.00 Median :162.9
## Mean :45.96 Mean :163.7
## 3rd Qu.:57.25 3rd Qu.:170.2
## Max. :70.00 Max. :180.5
## weight income generation
## Min. :28.30 Min. : 24.0 Length:100
## 1st Qu.:48.95 1st Qu.: 134.8 Class :character
## Median :59.95 Median : 298.5 Mode :character
## Mean :59.18 Mean : 434.4
## 3rd Qu.:67.33 3rd Qu.: 607.2
## Max. :85.60 Max. :2351.0library(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerstargazer(as.data.frame(income),
type = "html")| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| age | 100 | 45.960 | 13.331 | 20 | 36 | 57.2 | 70 |
| height | 100 | 163.746 | 7.692 | 148.000 | 158.100 | 170.175 | 180.500 |
| weight | 100 | 59.179 | 12.647 | 28.300 | 48.950 | 67.325 | 85.600 |
| income | 100 | 434.400 | 445.777 | 24 | 134.8 | 607.2 | 2,351 |
hist(income$height,freq = FALSE,
xlab = "身長",
ylab = "確率密度",
main = "身長の分布")
library(ggplot2) #ggplotパッケージをロードggplot(income,aes(height))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
theme_classic()
ggplot(income,aes(height))+
geom_histogram(aes(y = ..density..),
bins = 10,
colour = "gray",
fill = "blue")+
labs(title = "身長の分布",
x="身長",
y="確率密度")+
theme_classic()
stem(income$height)#チャンクオプションにcomment=""を指定##
## The decimal point is 1 digit(s) to the right of the |
##
## 14 | 88
## 15 | 001112334
## 15 | 555667777788888899999999
## 16 | 00001112222233333444
## 16 | 556666777778999999
## 17 | 0001111222233333344
## 17 | 5555889
## 18 | 1stem(income$height,scale=0.5)#チャンクオプションにcomment=""を指定##
## The decimal point is 1 digit(s) to the right of the |
##
## 14 | 88
## 15 | 001112334555667777788888899999999
## 16 | 00001112222233333444556666777778999999
## 17 | 00011112222333333445555889
## 18 | 1ggplot(income,aes(x = sex,y = height)) +
geom_boxplot() +
labs(x="性別",y="身長")+
scale_x_discrete(labels=c("女性","男性"))+
theme_classic()
ggplot(income,aes(x = sex,y = height)) +
geom_boxplot() +
coord_flip()+
labs(y="身長",x="性別")+
scale_x_discrete(labels=c("女性","男性"))+
theme_classic()
##ｐ113
ggplot(income,aes(x = sex,y = height)) +
geom_violin() +
geom_point()+
labs(x="性別",y="身長")+
scale_x_discrete(labels=c("女性","男性"))+
theme_classic()
##ｐ114
ggplot(income,aes(height,weight)) +
geom_point()+
labs(x="身長",y="体重")+
theme_classic()
##ｐ115間違いstat_smoothnの引数を入れないとこうなる
ggplot(income,aes(height,weight)) +
geom_point()+
stat_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
labs(x="身長",y="体重") +
theme_classic()## NULL##ｐ115-116 身長と体重の相関図 ##ggplot2のstat_smooth()のmethod引数で直線回帰を指定。 #＃https://rpubs.com/uri-sy/ggplot2_stat_smooth参照 ##エラー箇所を= lmを= “lm”に修正したら回帰直線がでるようになった。
ggplot(income,aes(height,weight)) +
geom_point()+
stat_smooth(method = "lm", se = FALSE, colour = "black", size = 1) +
labs(x="身長",y="体重") +
theme_classic()
#se = TRUEにするこで初期設定で95%信頼区間も表示される。
ggplot(income,aes(height,weight)) +
geom_point()+
stat_smooth(method = "lm", se = TRUE, colour = "black", size = 1) +
labs(x="身長",y="体重") +
theme_classic()
##ｐ117 身長と体重の観測点95％信頼区間有
ggplot(income,aes(height,weight)) +
geom_point()+
stat_smooth(method = "lm", se = TRUE, colour = "black", size = 1) +
geom_text(aes(y=weight+1,label=id),size=3,vjust=0)+
labs(x="身長（身長と体重の観測点95％信頼区間有）",y="体重") +
theme_classic()
##ｐ117 身長と体重の観測点p2
ggplot(income,aes(height,weight)) +
geom_point()+
stat_smooth(method = "lm", se = FALSE, colour = "black", size = 1) +
geom_text(aes(y=weight+1,label=id),size=3,vjust=0)+
labs(x="身長",y="体重") +
theme_classic()
##ｐ117 身長と体重の観測点 間違いサンプル
ggplot(income,aes(height,weight)) +
geom_point(aes(color=sex))+
geom_smooth(se=FALSE,aes(color=sex))+
scale_color_grey()+
labs(x="身長",y="体重") +
theme_classic()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
library("dplyr")##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary("ggplot2")library("DT")## Warning: package 'DT' was built under R version 3.6.3library("broom")##stat_smoothwをgeom_smoothnにしても同じ直線が引ける ##なお、ハッチのかかっている範囲は「母集団から標本を取ってきて、その平均から95%信頼区間を求める、という作業を100回やったときに、95回はその区間の中に母平均が含まれる」という意味です。 https://bellcurve.jp/statistics/course/8891.html
ggplot(income,aes(height,weight)) +
geom_point()+
geom_smooth(method = "lm", se = TRUE, colour = "black", size = 1) +
labs(x="身長",y="体重") +
theme_classic()
##ｐ118図4.46 身長と体重の散布図（男女色別） scale_color_grey()+
ggplot(income,aes(x = height,y = weight,shape = sex)) +
geom_point(aes(color=sex))+
stat_smooth(method = "lm", se = FALSE,aes(color=sex))+
labs(x="身長",y="体重") +
theme_classic()
##ｐ119 図4.46 身長と体重の散布図（男女ドットの濃淡比） scale_color_grey()+で色を消す
ggplot(income,aes(height,weight)) +
geom_point(aes(color=sex))+
stat_smooth(method = "lm", se = FALSE,aes(color=sex))+
scale_color_grey()+
labs(x="身長",y="体重") +
theme_classic()
##ｐ120 折れ線グラフ
library(gapminder)data("gapminder")glimpse(gapminder)## Observations: 1,704
## Variables: 6
## $ country <fct> Afghanistan, Afghanistan, Afghanistan, Afghanistan, Afgha...
## $ continent <fct> Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asi...
## $ year <int> 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 199...
## $ lifeExp <dbl> 28.801, 30.332, 31.997, 34.020, 36.088, 38.438, 39.854, 4...
## $ pop <int> 8425333, 9240934, 10267083, 11537966, 13079460, 14880372,...
## $ gdpPercap <dbl> 779.4453, 820.8530, 853.1007, 836.1971, 739.9811, 786.113...head(gapminder)## # A tibble: 6 x 6
## country continent year lifeExp pop gdpPercap
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Afghanistan Asia 1952 28.8 8425333 779.
## 2 Afghanistan Asia 1957 30.3 9240934 821.
## 3 Afghanistan Asia 1962 32.0 10267083 853.
## 4 Afghanistan Asia 1967 34.0 11537966 836.
## 5 Afghanistan Asia 1972 36.1 13079460 740.
## 6 Afghanistan Asia 1977 38.4 14880372 786.tail(gapminder)## # A tibble: 6 x 6
## country continent year lifeExp pop gdpPercap
## <fct> <fct> <int> <dbl> <int> <dbl>
## 1 Zimbabwe Africa 1982 60.4 7636524 789.
## 2 Zimbabwe Africa 1987 62.4 9216418 706.
## 3 Zimbabwe Africa 1992 60.4 10704340 693.
## 4 Zimbabwe Africa 1997 46.8 11404948 792.
## 5 Zimbabwe Africa 2002 40.0 11926563 672.
## 6 Zimbabwe Africa 2007 43.5 12311143 470.summary(gapminder)## country continent year lifeExp
## Afghanistan: 12 Africa :624 Min. :1952 Min. :23.60
## Albania : 12 Americas:300 1st Qu.:1966 1st Qu.:48.20
## Algeria : 12 Asia :396 Median :1980 Median :60.71
## Angola : 12 Europe :360 Mean :1980 Mean :59.47
## Argentina : 12 Oceania : 24 3rd Qu.:1993 3rd Qu.:70.85
## Australia : 12 Max. :2007 Max. :82.60
## (Other) :1632
## pop gdpPercap
## Min. :6.001e+04 Min. : 241.2
## 1st Qu.:2.794e+06 1st Qu.: 1202.1
## Median :7.024e+06 Median : 3531.8
## Mean :2.960e+07 Mean : 7215.3
## 3rd Qu.:1.959e+07 3rd Qu.: 9325.5
## Max. :1.319e+09 Max. :113523.1
## ##p123 図4.50
table(gapminder$country)##
## Afghanistan Albania Algeria
## 12 12 12
## Angola Argentina Australia
## 12 12 12
## Austria Bahrain Bangladesh
## 12 12 12
## Belgium Benin Bolivia
## 12 12 12
## Bosnia and Herzegovina Botswana Brazil
## 12 12 12
## Bulgaria Burkina Faso Burundi
## 12 12 12
## Cambodia Cameroon Canada
## 12 12 12
## Central African Republic Chad Chile
## 12 12 12
## China Colombia Comoros
## 12 12 12
## Congo, Dem. Rep. Congo, Rep. Costa Rica
## 12 12 12
## Cote d'Ivoire Croatia Cuba
## 12 12 12
## Czech Republic Denmark Djibouti
## 12 12 12
## Dominican Republic Ecuador Egypt
## 12 12 12
## El Salvador Equatorial Guinea Eritrea
## 12 12 12
## Ethiopia Finland France
## 12 12 12
## Gabon Gambia Germany
## 12 12 12
## Ghana Greece Guatemala
## 12 12 12
## Guinea Guinea-Bissau Haiti
## 12 12 12
## Honduras Hong Kong, China Hungary
## 12 12 12
## Iceland India Indonesia
## 12 12 12
## Iran Iraq Ireland
## 12 12 12
## Israel Italy Jamaica
## 12 12 12
## Japan Jordan Kenya
## 12 12 12
## Korea, Dem. Rep. Korea, Rep. Kuwait
## 12 12 12
## Lebanon Lesotho Liberia
## 12 12 12
## Libya Madagascar Malawi
## 12 12 12
## Malaysia Mali Mauritania
## 12 12 12
## Mauritius Mexico Mongolia
## 12 12 12
## Montenegro Morocco Mozambique
## 12 12 12
## Myanmar Namibia Nepal
## 12 12 12
## Netherlands New Zealand Nicaragua
## 12 12 12
## Niger Nigeria Norway
## 12 12 12
## Oman Pakistan Panama
## 12 12 12
## Paraguay Peru Philippines
## 12 12 12
## Poland Portugal Puerto Rico
## 12 12 12
## Reunion Romania Rwanda
## 12 12 12
## Sao Tome and Principe Saudi Arabia Senegal
## 12 12 12
## Serbia Sierra Leone Singapore
## 12 12 12
## Slovak Republic Slovenia Somalia
## 12 12 12
## South Africa Spain Sri Lanka
## 12 12 12
## Sudan Swaziland Sweden
## 12 12 12
## Switzerland Syria Taiwan
## 12 12 12
## Tanzania Thailand Togo
## 12 12 12
## Trinidad and Tobago Tunisia Turkey
## 12 12 12
## Uganda United Kingdom United States
## 12 12 12
## Uruguay Venezuela Vietnam
## 12 12 12
## West Bank and Gaza Yemen, Rep. Zambia
## 12 12 12
## Zimbabwe
## 12library(dplyr)Japan <- gapminder %>%
filter(country=="Japan") %>%
select(year,lifeExp) #yearとlifeexpだけを抜き出す##p126
ggplot(Japan,aes(x=year,y=lifeExp))+
geom_point()+
geom_line()+
ggtitle("日本人の平均寿命")+
labs(x="年",y="平均寿命")+
theme_classic()
##ｐ128 filter(country==“Japan” | ##country==“chaina”)の中の|はorを表す：https://kazutan.github.io/kazutanR/hands_on_170730/filter.html
jpn.chi <- gapminder %>%
filter(country=="China" | country=="Japan") %>%
select(year,country,lifeExp) #year,country,lifeExpだけを抜き出すjpn.chi## # A tibble: 24 x 3
## year country lifeExp
## <int> <fct> <dbl>
## 1 1952 China 44
## 2 1957 China 50.5
## 3 1962 China 44.5
## 4 1967 China 58.4
## 5 1972 China 63.1
## 6 1977 China 64.0
## 7 1982 China 65.5
## 8 1987 China 67.3
## 9 1992 China 68.7
## 10 1997 China 70.4
## # ... with 14 more rowstail(jpn.chi)## # A tibble: 6 x 3
## year country lifeExp
## <int> <fct> <dbl>
## 1 1982 Japan 77.1
## 2 1987 Japan 78.7
## 3 1992 Japan 79.4
## 4 1997 Japan 80.7
## 5 2002 Japan 82
## 6 2007 Japan 82.6##p129 図4.57日本人と中国人の平均寿命の折れ線グラフ
ggplot(jpn.chi, aes(x=year, y=lifeExp, group=country)) +
geom_point() +
geom_line(aes(linetype=country)) +
scale_linetype_manual(values=c("twodash","dottend")) +
ggtitle("平均寿命（日本と中国)1952-2007") +
labs(x="年", y="平均寿命") +
scale_linetype_discrete(name="国名",labels = c("日本","中国"))+
theme_classic()## Scale for 'linetype' is already present. Adding another scale for 'linetype',
## which will replace the existing scale.
##ｐ131 図4.58衆院選データの最初の6行
library(readr)hr <- read_csv("hr96-17.csv",na = ".")## Parsed with column specification:
## cols(
## year = col_double(),
## ku = col_character(),
## kun = col_double(),
## party = col_character(),
## party_code = col_double(),
## name = col_character(),
## age = col_double(),
## nocand = col_double(),
## vote = col_double(),
## voteshare = col_double(),
## eligible = col_double()
## )head(hr)## # A tibble: 6 x 11
## year ku kun party party_code name age nocand vote voteshare eligible
## <dbl> <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1996 aichi 1 NFP 8 KAWA~ 47 7 66876 40 346774
## 2 1996 aichi 1 LDP 1 IMAE~ 72 7 42969 25.7 346774
## 3 1996 aichi 1 DPJ 3 SATO~ 53 7 33503 20.1 346774
## 4 1996 aichi 1 JCP 2 IWAN~ 43 7 22209 13.3 346774
## 5 1996 aichi 1 othe~ 100 ITO,~ 51 7 616 0.4 346774
## 6 1996 aichi 1 koku~ 22 YAMA~ 51 7 566 0.3 346774library(dplyr)
shinzo <- hr %>%
filter(name=="ABE, SHINZO") %>%
select(year,ku,kun,party,age,nocand,vote,voteshare)shinzo## # A tibble: 8 x 8
## year ku kun party age nocand vote voteshare
## <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 1996 yamaguchi 4 LDP 42 3 93459 54.3
## 2 2000 yamaguchi 4 LDP 45 2 121835 71.7
## 3 2003 yamaguchi 4 LDP 49 3 140347 79.7
## 4 2005 yamaguchi 4 LDP 50 3 137701 73.6
## 5 2009 yamaguchi 4 LDP 54 3 121365 64.3
## 6 2012 yamaguchi 4 LDP 58 3 118696 78.2
## 7 2014 yamaguchi 4 LDP 60 3 100829 76.3
## 8 2017 yamaguchi 4 LDP 63 5 104825 72.6ここまでで横浜市立大学の講義は終了以降は、http://rmecab.jp/ranko/chap2.htmlを参考図書として実施
-統計の基礎 -正規分布のグラフ作成 もっとも有名な確率分布です．身長体重，学校の成績など，現実の多くのデータを確率的に説明するのに使われる分布です． 下の図のような曲線になるのが特徴です． curve (dnorm, from = -4, to = 4)
curve (dnorm, from = -4, to = 4)
curve (dnorm, from = -10, to = 10)
- ヒストグラム
ヒストグラムは，データを区間ごとに分割して，それぞれ区間の属するデータの数をY軸にとったグラフです．X軸は，区間を表します．
結果として，棒グラフが並んだグラフとなります．このとき，中央の棒が大きく，左右対象に棒が低くなっている場合，データを正規分布にあてはめて考えることができます．
ここで，正規屋さんの一ヶ月の売上をヒストグラムにしてみます．ヒストグラムにしたとき，正規分布の曲線におおむね一致しているようなら，確率分布として正規分布を想定して分析を行うことができます．
install.packages (“Ranko”, repos = “http://rmecab.jp/R”) ; demo(startup)
『とある弁当屋の統計技師（データサイエンティスト）』
付録パッケージ：バージョン2.1
第1章前半のデモを実行するには demo(chap1) と入力してEnterを押します
第2章であれば demo(chap2) と数値部分を変えて入力します
```r
library(Ranko)##
## 『とある弁当屋の統計技師（データサイエンティスト）』
## 付録パッケージ：バージョン2.1
##
##
## 第1章前半のデモを実行するには demo(chap1) と入力してEnterを押します
## 第2章であれば demo(chap2) と数値部分を変えて入力します##
## Attaching package: 'Ranko'## The following object is masked from 'package:utils':
##
## demodemo(chap1)##
##
## demo(chap1)
## ---- ~~~~~
## デモを途中でやめる場合はEscキーを押します
##
## > # 文太：これからRで本書の内容をおさらいします
## > # 1章から6章までの内容を，個別に確認します．
## > # demo() の丸括弧の中に chap1 とかを入力して
## > # 閉じ括弧の後ろでEnterキーを押していただければいいのです
## >
## >
## > # 乱子：このコンソールとか言うWindowの右端でEnter押すのね
## > # 文太：そうです．
## > # 2章以降は demo (chap2); demo(chap3) などとなります
## > #
## > # では，1週間分のお弁当データの確認から始めますね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太: Rは複数の数値をひとかたまりのデータとして扱います
## >
## > data(bento7); bento7
## 月 火 水 木 金 土 日
## 181 194 265 206 208 272 221
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：なに，この bento7 というのは？
## > # 文太：Rでは，データに名前をつけて扱います．この場合は
## > # 181 194 206 208 221 265 272
## > # の7つの数字をまとめた塊に bento7 という名前をつけたのです
## > # 乱子：相変わらずセンス無いわね
## > # 文太：XとかYよりは分かりやすいと思ったんですけどね．．．
## > # で，統計解析ではこうしたデータの塊につけた名前を
## > # 変数あるいはオブジェクトといいます
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：いま，このデータ(変数名bento7)をある週における
## > # 1日のお弁当売上個数だと考えましょう
## > # 平均を求めてみます．以下の式を計算するわけです
## >
## > 弁当平均の計算式()
## ( 181 + 194 + 265 + 206 + 208 + 272 + 221 ) / 7
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：こんな長ったらしい式を入力するわけ？
## > # 文太：いや，まさか
## > # Rを使えば，もっと簡単に平均値を求められます．
## > # 平均値のことを英語で mean といいます．
## > # Rで平均値を求めるにはmean関数を使います．
## > # 乱子：関数ってなによ
## > # 文太：関数というのは，hogehoge() みたいな記号です
## > # コードともいいます
## > # mean()は平均値を求める関数です
## > # 丸括弧内にデータの塊を表す変数(bento7)を指定して
## > # こんなふうに実行します
## >
## > mean (bento7)
## [1] 221
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：このように関数に続く丸括弧の内部に，
## > # データのかたまりを表す変数(bento7)を指定して
## > # Enterキーを押すと次の行に結果が表示されます
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：次に中央値です．
## > # 中央値は並び替えて真ん中に位置するデータです
## > # まず並びかえてみます
## >
## > sort (bento7)
## 月 火 木 金 日 水 土
## 181 194 206 208 221 265 272
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：Rではsort関数に続く丸括弧内に変数を指定すると
## > # 並びかえた結果を表示してくれます
## > # 中央に位置している値が中央値です
## > #
## > # 乱子：この変数の場合は7つしか数値がないから良いけど
## > # たくさんデータがあったら，
## > # どこが真ん中だかわからないけど？
## > #
## > # 文太：もちろん，中央値を求める関数もありますって
## > # 中央値のことを英語で median といいます
## > # Rではmedianという関数を使うと中央値が表示されます
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > median (bento7)
## [1] 208
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ここで土曜日の売上を272個から500個に変更します
## > # Rではデータのかたまり変数から一部を指定できます
## >
## > bento7 [6] <- 500
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：なによ，この矢印 (<-) は？
## > # 文太：代入記号といいます．矢印 (<-)の左にあるものと
## > # 右にあるものをひもづけるわけです
## > # たとえば x <- 8 とすると，x に 8 が代入されて，以降
## > # x は 8 を指していることになります
## >
## > x <- 8
##
## > x
## [1] 8
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：bento7変数の後ろにある "[6]" ってのは？
## > # 文太：bento7 [6] とすると，変数bento7の6番目の要素
## > # つまり土曜日のデータを指定したことになります
## >
## > bento7 [6]
## 土
## 500
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：bento7 [6] <- 500 というのは，変数 bento7 の
## > # 6番目の要素を 500に置き換えるという命令（コード）です
## > # ではあらためて平均値と中央値を求めてみます
## >
## > mean (bento7)
## [1] 253.5714
##
## > median (bento7)
## [1] 208
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：なるほど，平均値は前と比べて大きくなったけど
## > # 中央値は変わらないわけね
## > # 文太：そうなんです．
## > # 乱子：これでおしまい？
## > # 文太：いえいえ，これからが本番です
## > # それではバラツキについて復習しましょう
## > # グラフを使ってお話ししたいと思います
## > # その前にグラフ作成機能を拡張しましょう
## > # 乱子：拡張？
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > # 文太：はい．Rでデフォルトで入っているグラフも
## > # 悪くないんですが，もう少し見栄えのいい
## > # グラフを作成したいと思います
## > # そのための「ggplot2パッケージ」を導入します
## >
## > # 乱子：ggplot2？
## > # 文太：Rのグラフィックス機能を拡張する有名なパッケージです
## > # Rではこういう拡張機能をパッケージと言って
## > # 後からいくつでも自由に無料で追加できます
## > # 一度インストールすれば，以降は何度でも自由に使えます
## > # 使う直前に一回だけこんな1行を実行しておきます
## > # library (ggplot2)
## > # 乱子：library()ってのも関数？
## > # 文太：そうです
## > # いま実行したのは，これからggplot2パッケージを使うかんな！
## > # って宣言です
## > # では本文でも紹介した箱ヒゲ図を作ってみましょうか
## > # 正規屋さんとガンマ屋さんそれぞれのお弁当の売上1ヶ月分でした
## > # データは bento という名前だとします
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > data(bento)
##
## > ## if (.Platform$pkgType == "mac.binary"){
## > ## quartzFonts(HiraMaru=quartzFont(rep(c("HiraMaruProN-W3","HiraMaruProN-W6"), 2)))
## > ## quartzFonts(HiraKaku=quartzFont(rep(c("HiraKakuProN-W3","HiraKakuProN-W6"), 2)))
## > ## theme_set(theme_grey(base_family = "HiraKaku"))
## > ## }else if () (.Platform$pkgType == "win.binary"){
## > ## windowsFonts(MEI = windowsFont("Meiryo"))
## > ## theme_set(theme_grey(base_family = "MEI"))
## > ## }else {
## > ## theme_set(theme_gray(base_size = 12,base_family="Japan1"))
## > ## }
## >
## >
## > m <- ggplot(bento, aes(factor (お店), 売上個数))
##
## > m + geom_boxplot(aes(fill = factor(お店))) + xlab("とあるお弁当屋さん２軒") + ylab ("１週間の売上")
##
## > ##############################
## > # 乱子：グラフはたしかにできたけど，この上に書いてある命令，
## > # コード？まったく意味わからないけど
## > # 文太：ええ．ggplot2パッケージの使い方にはコツがあるんです
## > # もっと簡単に箱ヒゲ図を作る方法もあります
## > # boxplot() 関数を使うんです
## > #
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > boxplot(売上個数 ~ お店, bento)
##
## > ##############################
## > # 乱子：地味ね．．．
## > # でもboxplot関数の丸括弧内にデータ名 bento を指定しただけ
## > # これなら想像つくかな
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それぞれのお弁当屋さん売上個数を確認しておきましょうか
## > # まず正規屋さんの先月の売上データ
## > bento [ bento$お店 == "正規屋", ]
## 売上個数 お店
## 1 181 正規屋
## 2 194 正規屋
## 3 265 正規屋
## 4 206 正規屋
## 5 208 正規屋
## 6 272 正規屋
## 7 221 正規屋
## 8 152 正規屋
## 9 176 正規屋
## 10 185 正規屋
## 11 252 正規屋
## 12 217 正規屋
## 13 219 正規屋
## 14 207 正規屋
## 15 181 正規屋
## 16 274 正規屋
## 17 223 正規屋
## 18 124 正規屋
## 19 231 正規屋
## 20 184 正規屋
## 21 160 正規屋
## 22 194 正規屋
## 23 162 正規屋
## 24 174 正規屋
## 25 178 正規屋
## 26 136 正規屋
## 27 237 正規屋
## 28 209 正規屋
## 29 157 正規屋
## 30 253 正規屋
##
## > # 文太：次に（架空の）がんま亭の先月の売上データ
## > head (bento [ bento$お店 == "がんま亭", ], 30)
## 売上個数 お店
## 31 207 がんま亭
## 32 204 がんま亭
## 33 209 がんま亭
## 34 209 がんま亭
## 35 209 がんま亭
## 36 208 がんま亭
## 37 208 がんま亭
## 38 205 がんま亭
## 39 203 がんま亭
## 40 203 がんま亭
## 41 202 がんま亭
## 42 204 がんま亭
## 43 199 がんま亭
## 44 216 がんま亭
## 45 211 がんま亭
## 46 199 がんま亭
## 47 203 がんま亭
## 48 203 がんま亭
## 49 209 がんま亭
## 50 205 がんま亭
## 51 206 がんま亭
## 52 205 がんま亭
## 53 205 がんま亭
## 54 212 がんま亭
## 55 204 がんま亭
## 56 213 がんま亭
## 57 197 がんま亭
## 58 208 がんま亭
## 59 206 がんま亭
## 60 206 がんま亭
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それぞれのお店の中央値や分位数です
## >
## > by(bento$売上個数, bento$お店, summary)
## bento$お店: がんま亭
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 197.0 203.2 205.5 205.9 208.8 216.0
## ------------------------------------------------------------
## bento$お店: 正規屋
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 124.0 176.5 200.0 201.1 222.5 274.0
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：関数の方はともかく，表示内容を説明して？
## > # 文太：そうですね (^^;
## > # Rとその関数や，使い方などについては
## > # 石田基広『Rで学ぶデータ・プログラミング入門』
## > # 共立出版 ISBN 4320110293
## > # をぜひ参照して下さい
## > #
## > # 乱子：出力の説明！
## > # 文太：あ，はい！
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > # まず「bento$お店: がんま亭」は，そのまんまbentoデータから
## > # お店ががんま亭のデータを対象にした，ということです
## > #
## > # 「Min. 1st Qu. Median Mean 3rd Qu. Max. 」は左から
## > # 「最小値 第1四分位数 中央値 平均値 第3四分位数 最大数」
## > # に対応します
## > #
## > # 乱子：我が正規屋の場合，一番売れ行きの悪かった日は 124個
## > # 一番良かった日は 274個，中央値は 200 個
## > # 第1四分位数は 176.5 で第3四分位数は 222.5 ってわけね
## > # 文太：はい，そうです
## > # じゃあ，今度は正規屋さんのオムライスとオムライス弁当
## > # それぞれの売上個数を箱ヒゲ図にしてみましょう
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > data (omu)
##
## > m <- ggplot(omu, aes(factor (弁当名), 売上個数))
##
## > m + geom_boxplot(aes(fill = factor(弁当名))) + xlab("オムライス2 種の売上数") + ylab ("1 ヶ月の売上")
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太，それぞれのお弁当の要約も求めてみましょう
## >
## > by(omu$売上個数, omu$弁当名, summary)
## omu$弁当名: オムライス
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 7.00 8.00 8.80 9.75 18.00
## ------------------------------------------------------------
## omu$弁当名: オムライス弁当
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 11.00 12.50 13.07 14.75 21.00
##
## > browseURL("http://rmecab.jp/ranko/chap1.html")
##
## > # ここで一旦休憩しましょう！
## > # 続けて2章後半のデモをみる場合は
## > # demo (chap2)
## > # と入力してEnterを押してください
## > # なおRを終了させようとすると
## > # 「作業スペースを保存しますか？」
## > # と尋ねてきますが，これには
## > # 「いいえ」を選んで構いません（「はい」でも問題はありません）．
## > # Rをあらためて起動して，本書のデモを実行する際は，最初に
## > # library (Ranko)
## > # を実行してください
## >demo(chap2)##
##
## demo(chap2)
## ---- ~~~~~
## デモを途中でやめる場合はEscキーを押します
##
## > options(scipen=1)
##
## > library (grid)
##
## > ##############################
## > # 文太：第2章ではデータ分析の王道を学びます！
## > # 乱子：王道ってのが，なんか怪しいけど．．．
## >
## >
## >
## > # 文太：それでは正規屋さんの先月の売上を表示してみます
## >
## > data (bento30); bento30
## [1] 90500 97000 132500 103000 104000 136000 110500 76000 88000 92500
## [11] 126000 108500 109500 103500 90500 137000 111500 62000 115500 92000
## [21] 80000 97000 81000 87000 89000 68000 118500 104500 78500 126500
##
## > # 乱子：端数は省略されているのね
## > # 文太：はい
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > options(digits=8)
##
## > ##############################
## > # 文太：さてRのmean関数を使って平均値と中央値を求めます
## >
## > mean (bento30)
## [1] 100533.33
##
## > # 中央値
## > median (bento30)
## [1] 100000
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > # options (digits=7)
## > ##############################
## > # 文太：次に平均値に対するバラツキを考えてみます
## > # 乱子：中央値の場合の四分位範囲みたいなやつね
## > # 文太：そうです
## > # まずデータのそれぞれから平均値を引いてみます
## >
## > bento30 - mean (bento30)
## [1] -10033.3333 -3533.3333 31966.6667 2466.6667 3466.6667 35466.6667
## [7] 9966.6667 -24533.3333 -12533.3333 -8033.3333 25466.6667 7966.6667
## [13] 8966.6667 2966.6667 -10033.3333 36466.6667 10966.6667 -38533.3333
## [19] 14966.6667 -8533.3333 -20533.3333 -3533.3333 -19533.3333 -13533.3333
## [25] -11533.3333 -32533.3333 17966.6667 3966.6667 -22033.3333 25966.6667
##
## > # 文太：引き算の結果を自乗します
## > (bento30 - mean (bento30))^2
## [1] 100667777.8 12484444.4 1021867777.8 6084444.4 12017777.8
## [6] 1257884444.4 99334444.4 601884444.4 157084444.4 64534444.4
## [11] 648551111.1 63467777.8 80401111.1 8801111.1 100667777.8
## [16] 1329817777.8 120267777.8 1484817777.8 224001111.1 72817777.8
## [21] 421617777.8 12484444.4 381551111.1 183151111.1 133017777.8
## [26] 1058417777.8 322801111.1 15734444.4 485467777.8 674267777.8
##
## > # 文太：その合計を求めます
## > # Rで合計は sum()関数で求めます
## >
## > sum( (bento30 - mean (bento30))^2)
## [1] 11155966667
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：分散はこの合計をデータ数で割ります
## >
## >
## > 弁当の分散 <- sum( (bento30 - mean (bento30))^2) / 30
##
## > 弁当の分散
## [1] 371865556
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：その平方根が「標準偏差」です
## > # 平方根はsqrt()関数で求めます
## >
## > sqrt (弁当の分散)
## [1] 19283.816
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ただしですね，本書では説明しませんでしたが
## > # 不偏分散という考え方があります
## > # 本書あとがきで説明していますが
## > # データ数ではなく，データ数-1 で割った数値です
## > # 多くの統計ソフトで分散といえば不偏分散を意味します
## > # そして不偏分散を求めるための関数が用意されています
## > # Rだと var() 関数です．標準偏差は sd() 関数です
## >
## > var (bento30)
## [1] 384688506
##
## > sd (bento30)
## [1] 19613.478
##
## > # 乱子：普通はvar() 関数と sd() 関数を使って構わないのね
## > # 文太：そうですが，一応，どちらの分散が出力されているのかは
## > # 注意してください
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：お弁当売上のヒストグラムを描いてみます
## > # データを整理しなおして bentoDF と名前を変えました
## >
## > data (bentoDF)
##
## > m <- ggplot(bentoDF, aes( x = 売上))
##
## > m + geom_histogram(colour = "black", fill = "white", binwidth = 10000) +
## + labs(title = "1日の平均売上数のヒストグラム") +
## + scale_x_continuous(breaks=seq(50000,150000,by=10000))
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：実際にお弁当売上について本当の平均値範囲を
## > # 計算してみましょう
## >
## >
## > c (100533.3 - 1.96 * (19283.82) / sqrt (30),
## + 100533.3 + 1.96 * (19283.82) / sqrt (30) )
## [1] 93632.674 107433.926
##
## > # 乱子：この範囲に，本当の平均値が含まれているのね
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：仮想の受験者たちの偏差値データを
## > # ヒストグラムで表してみます
## >
## > data (hensati)
##
## > ggplot(hensati, aes(x = 偏差値)) +
## + geom_histogram(fill="white", colour="black", binwidth=5) +
## + geom_density(aes(y=6*..count..), colour="black", adjust=4) +
## + labs (title="偏差値は平均が50：標準偏差が10") +
## + theme(panel.grid.major = element_blank ()) +
## + xlab("偏差値") + ylab("人数") +
## + scale_x_continuous(breaks=seq(10,90,by=10))
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：平均が0で標準偏差が１の正規分布を表すデータを
## > # グラフにしてみます
## >
## > data (normalDist)
##
## > p <- ggplot(normalDist, aes(x = 正規分布, y = 確率密度)) +
## + geom_line(colour="black", size = 1)
##
## > p + labs(title="平均が0で標準偏差が１の正規分布")
##
## > # 乱子：曲線の下の面積を考えた場合，
## > # x軸で-1.96 から 1.96 の範囲の面積が 0.95，
## > # つまり95%なのね
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：そうです．95%範囲を塗りつぶしてみます
## >
## > DF2 <- normalDist[normalDist$確率密度 > 0.058,] #0.05844
##
## > p + geom_ribbon(data = DF2 , fill = "skyblue", aes ( ymin= 0, ymax=確率密度)) + labs(title="平均が0で標準偏差が１の正規分布") + annotate("text", label = "この着色した面積が全体の95%", x = 0, y = 0.05, family = "JP1") # + annotate("text", label = "-1.96", x = -1.96, y = .0008,family = "JP1") + annotate("text", label = "1.96", x = 1.96, y = .0008,family = "JP1") + geom_segment(aes(x = -1.9, y = 0.01, xend = 1.9, yend = 0.01), arrow = arrow(ends = "both"))
##
## > # , length = unit(0.5, "cm")))
## >
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：偏差値の理論グラフを書いてみます
## >
## > data (hensati2)
##
## > p <- ggplot(hensati2, aes(x = 偏差値, y = 確率密度)) +
## + geom_line(colour="black", size = 1)
##
## > p + labs(title="偏差値でデータの広がりは±標準偏差") +
## + xlab("偏差値") + ylab("")
##
## > # 乱子：受験生の95%は偏差値が30から70の間にいるというわけね
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：同じく，この偏差値についても95%範囲を
## > # 塗りつぶしてみますね
## >
## > DF2 <- hensati2 [ hensati2$確率密度 > 0.0053,]
##
## > p + labs(title = "偏差値でデータの広がりは±標準偏差") + xlab("偏差値")+ ylab("") + geom_ribbon(data = DF2 ,fill="greenyellow", aes ( ymin= 0, ymax=確率密度)) + annotate("text", label = "この着色した面積が全体の95%", x = 50, y = 0.005, family = "JP1") + annotate("text", label = "30", x = 30, y = .0008,family = "JP1") + annotate("text", label = "70", x = 70, y = .0008,family = "JP1") + geom_segment(aes(x = 30, y = 0.003, xend = 70, yend = 0.003), arrow = arrow(ends = "both", length = unit(0.5, "cm")))
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それでは正規屋さんのお弁当売上高について
## > # 平均値と標準偏差を使って理論的な曲線グラフを
## > # 作成してみましょう
## >
## > data (obento)
##
## > p <- ggplot(obento, aes(x = お弁当売上, y = 確率密度)) + geom_line(colour="black", size = 1)
##
## > DF2 <- obento[obento$お弁当売上 > 93633 & obento$お弁当売上 < 107434, ]
##
## > p + labs(title="正規屋さんのお弁当売上の分布") + ylab("") + geom_ribbon(data = DF2 , fill="greenyellow", aes ( ymin= 0, ymax=確率密度)) + annotate("text", label = "この着色した面積が全体の95%", x = 100000, y = 20, family = "JP1") + annotate("text", label = "93633", x = 93633, y = 1, family = "JP1") + annotate("text", label = "107434", x = 107434, y = 1, family = "JP1") + geom_segment(aes(x = 93633, y = 10, xend = 107434, yend = 10), arrow = arrow(ends = "both"))# , length = unit(0.5, "cm")
##
## > # 乱子：まあ，だいたい93600円から107400円の間に
## > # 月売上の本当の平均値があるわけね
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：サザエさんの視聴率について標本平均値が 19.3%
## > # だったとして本当の視聴率の範囲を求めてみましょう
## >
## >
## > c(
## + 0.193 - 1.96 * sqrt ( (0.193 * (1 - 0.193)) / 600),
## + 0.193 + 1.96 * sqrt ( (0.193 * (1 - 0.193)) / 600)
## + )
## [1] 0.16142118 0.22457882
##
## > # 乱子：20%を超えている可能性があるわけだから
## > # サザエさん制作スタッフは左遷されないですむわけね
## >
## >
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：二項分布にもとづいてコインの表が出る枚数
## > # についてそれぞれの確率をグラフにしてみます
## >
## > data (coin)
##
## > p <- ggplot (coin, aes(x= 枚数, y = 確率))# + geom_bar(stat="bin")
##
## > p + geom_point(size=5) + geom_step() + labs(title="表が出る確率") +
## + xlab ("表の枚数") + ylab("確率") +
## + scale_x_discrete(breaks = seq(0, 10, by = 1))
##
## > # 文太：ちなみに10個から3個を選ぶ場合の数は
## > # Rではchoose()関数を使って求められます
## >
## > choose (10, 3)
## [1] 120
##
## > browseURL("http://rmecab.jp/ranko/chap2.html")demo(chap3)##
##
## demo(chap3)
## ---- ~~~~~
## デモを途中でやめる場合はEscキーを押します
##
## > ##############################
## > # 文太：第3章は相関と回帰分析について説明します
## > # とビールの消費量と気温の関係を例に説明しましょう
## >
## > data (beer)
##
## > beer[, 1:3]
## 月 気温 消費量
## 1 1 5.5 2.38
## 2 2 6.6 3.85
## 3 3 8.1 4.41
## 4 4 15.8 5.67
## 5 5 19.5 5.44
## 6 6 22.4 6.03
## 7 7 28.3 8.15
## 8 8 28.9 8.23
## 9 9 27.8 7.35
## 10 10 18.2 5.50
## 11 11 13.7 4.90
## 12 12 8.7 4.60
##
## > # 乱子：みんなビール好きよね．．．
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：消費量と気温という2つの変数からなるデータがある場合
## > # 散布図を描いてみると，それぞれの関係が直感的に分かります
## >
## >
## > p1 <- ggplot (beer, aes(x = 気温, y = 消費量))
##
## > p1 + geom_point(size = 5)+ labs(title = "ビール消費量と気温の関係図")
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ちょっと逆の関係を見てみましょう
## > # 駅からの距離と地価の関係です
## >
## > data (tika)
##
## > p <- ggplot(tika, aes(距離 , 地価) )
##
## > p + geom_point(size = 5) + labs(title = "駅からの距離と地価の相関")
##
## > # 乱子：こっちは距離が大きくなると地価は下がっていくのね
## > # 文太：ええ，負の相関があると表現します
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：グラフだけで判断していると
## > # こういうデータの場合
## > # 相関があるかどうか，はっきりとはわかりません
## >
## > data (anime)
##
## > p <- ggplot(anime, aes(アニメ , 勉強時間) )
##
## > p + geom_point(size = 5) +
## + labs(title = "アニメ視聴時間と勉強時間相関")
##
## > # 乱子：アニメ視聴時間が長くなれば，当然，
## > # 勉強時間は減ると思うけどねぇ．．．
## > # 文太：それは偏見でしょう
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：すこし面倒な計算を続けます
## > # まず偏差という数値です
## > # 個々のデータから平均値を引きます
## >
## > beer$消費量 - mean (beer$消費量)
## [1] -3.1625 -1.6925 -1.1325 0.1275 -0.1025 0.4875 2.6075 2.6875 1.8075
## [10] -0.0425 -0.6425 -0.9425
##
## > beer$気温 - mean (beer$気温)
## [1] -11.4583333 -10.3583333 -8.8583333 -1.1583333 2.5416667 5.4416667
## [7] 11.3416667 11.9416667 10.8416667 1.2416667 -3.2583333 -8.2583333
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：beer$消費量っていうので
## > # 間にあるドルマーク($)は何？
## > # 文太：これはbeerという名前で登録したデータ
## > # の消費量という列，という意味です
## > # beer というデータには，こういう列があります
## >
## > head (beer[, 1:3])
## 月 気温 消費量
## 1 1 5.5 2.38
## 2 2 6.6 3.85
## 3 3 8.1 4.41
## 4 4 15.8 5.67
## 5 5 19.5 5.44
## 6 6 22.4 6.03
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それぞれの偏差どうしを掛け算してみます
## > # まず消費量の偏差と，気温の偏差を
## > # ペアごとに掛け算してみます
## >
## > (beer$消費量 - mean (beer$消費量)) * ( beer$気温 - mean (beer$気温))
## [1] 36.236979167 17.531479167 10.032062500 -0.147687500 -0.260520833
## [6] 2.652812500 29.573395833 32.093229167 19.596312500 -0.052770833
## [11] 2.093479167 7.783479167
##
## > # 文太：次にミリリットルを単位とした偏差の積
## >
## > (beer$ミリ - mean (beer$ミリ)) * ( beer$気温 - mean (beer$気温))
## [1] 36236.979167 17531.479167 10032.062500 -147.687500 -260.520833
## [6] 2652.812500 29573.395833 32093.229167 19596.312500 -52.770833
## [11] 2093.479167 7783.479167
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：いまの計算をもう一度繰り返しますが
## > # 今度は合計します
## >
## > # リットル単位での偏差積和
## > sum ( (beer$消費量 - mean (beer$消費量)) * ( beer$気温 - mean (beer$気温)) )
## [1] 157.13225
##
## > # ミリリットル単位での偏差積和
## > sum ( (beer$ミリ - mean (beer$ミリ)) * ( beer$気温 - mean (beer$気温)) )
## [1] 157132.25
##
## > # リットル単位での消費量と気温の共分散
## > sum ( (beer$消費量 - mean (beer$消費量)) * ( beer$気温 - mean (beer$気温)) ) / 12
## [1] 13.094354
##
## > # ミリリットル単位での消費量と気温の共分散
## > sum ( (beer$ミリ - mean (beer$ミリ)) * ( beer$気温 - mean (beer$気温)) ) / 12
## [1] 13094.354
##
## > # 乱子：なんだか，ややっこしいわね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：さらにややっこしくて恐縮ですが
## > # 分散と標準偏差はこうやって計算します
## >
## > # リットル単位での消費量の分散
## > sum ( (beer$消費量 - mean (beer$消費量))^2) / 12
## [1] 2.7503854
##
## > # リットル単位での消費量の標準偏差
## > sqrt (sum ( (beer$消費量 - mean (beer$消費量))^2) / 12)
## [1] 1.6584286
##
## > # ミリリットル単位での消費量の分散
## > sum ( (beer$ミリ - mean (beer$ミリ))^2) / 12
## [1] 2750385.4
##
## > # ミリリットル単位での消費量の標準偏差
## > sqrt (sum ( (beer$ミリ - mean (beer$ミリ))^2) / 12)
## [1] 1658.4286
##
## > # 気温の分散
## > sum ( (beer$気温 - mean (beer$気温))^2) / 12
## [1] 68.634097
##
## > # リットル単位での消費量の標準偏差
## > sqrt (sum ( (beer$消費量 - mean (beer$消費量))^2) / 12)
## [1] 1.6584286
##
## > # 気温の標準偏差
## > sqrt (sum ( (beer$気温 - mean (beer$気温))^2) / 12)
## [1] 8.2845698
##
## > # 乱子：何しようとしているんだっけ？
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：いや，消費量と気温の相関を数字で表現しようと
## > # しているんでした
## > # ただ，データの単位に影響されないよう無名化した数値を
## > # 使いたかったので，ちょっと，ややっこしい計算を
## > # してもらったんです
## > # ここまで数値を使うと，以下のように計算できます
## >
## >
## > 13.1 / ( 1.66 * 8.28)
## [1] 0.95308771
##
## > 13094 / (1658.43 * 8.28)
## [1] 0.95355304
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：実は，Rではですね，相関係数を求める関数も
## > # ちゃんとあります．．．痛い！
## > # 乱子：それならそうと，早く言ってよ
## > # 文太：まあ，計算式を一度は確認しておくのもいいかな
## > # と思いまして．．．
## > # cor()関数というのに，コンマをはさんで2つのデータ
## > # を指定します
## >
## > cor (beer$消費量 , beer$気温)
## [1] 0.95305364
##
## > cor (beer$ミリ , beer$気温)
## [1] 0.95305364
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：無名化に関連して，データの標準化についても
## > # 確認します
## > # たとえばこんな小さなデータがあったとします
## >
## > (dat1 <- c(1, 2, 3, 4, 5))
## [1] 1 2 3 4 5
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：まず，それぞれから全体平均を引きます
## >
## > (dat1 - mean (dat1))
## [1] -2 -1 0 1 2
##
## > # 文太：次に標準偏差ですが，ここでは不偏分散にもとづく
## > # 標準偏差を使います
## > # 標準偏差は sd() 関数で求められます
## >
## > (dat1 - mean (dat1)) / sd (dat1)
## [1] -1.26491106 -0.63245553 0.00000000 0.63245553 1.26491106
##
## > # 文太：dat1 を 100倍しただけのデータを用意します
## >
## > (dat2 <- dat1 * 100)
## [1] 100 200 300 400 500
##
## > # 文太：標準化してみます
## >
## > (dat2 - mean (dat2)) / sd (dat2)
## [1] -1.26491106 -0.63245553 0.00000000 0.63245553 1.26491106
##
## > # 乱子：元のデータを100倍したところで
## > # 標準化してみると同じなのね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ビールの消費量と気温の関係に完全な相関が
## > # ある場合はこんなグラフになります
## >
## > beer.lm <- lm (消費量 ~ 気温, data = beer)
##
## > coefs <- coef(beer.lm)
##
## > p1 <- ggplot (beer, aes(x = 気温, y = 消費量の予測)) +
## + labs(title = "完全な相関！！")#
##
## > p1 + geom_point(size = 5) +
## + geom_abline(intercept = coefs[1], slope = coefs[2])
##
## > # 乱子：この直線は？
## > # 文太：回帰直線です．つまり 売上 = a * 気温 + b を表します
## > # 乱子：直線上にデータがのっているのね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：でも現実のデータの場合は，こんな散布図になります
## >
## >
## > p1 <- ggplot (beer, aes(x = 気温, y = 消費量))
##
## > p1 + geom_point(size = 5) +
## + geom_abline(intercept = coefs[1], slope = coefs[2]) +
## + labs(title="回帰直線を重ねる")
##
## > # 乱子：ようするに，現実のデータは数式では完全に
## > # 表現できないってことね
## > # 文太：ええ，でも，近い数式を導くことはできます
## > # それが回帰式です
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ついでに，観測データと回帰直線をY軸に平行な直線で結んでみます．
## > # （初版で「垂線」と誤って表記されていますが，違います）
## > # （垂線とは，回帰直線に対して垂直な線のことです．．．m(..)m ）
## >
## > p1 <- ggplot (beer, aes(x = 気温, y = 消費量))
##
## > p1 + geom_point(size = 5)+ labs(title = "Y軸と平行にデータと回帰直線を結ぶ") +
## + geom_segment(aes(x = beer$気温, y = beer$消費量, yend = beer$消費量の予測, xend = beer$気温 ))+
## + geom_abline(intercept = coefs[1], slope = coefs[2])
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：回帰分析は直線でデータを表す方法なので
## > # こういうデータには使えません
## >
## > data (df)
##
## > p <- ggplot(df, aes(あるデータ , もう一方のデータ) )
##
## > p + geom_point(size = 5) + labs(title = "相関がある？ない？") +
## + theme(panel.grid.major = element_blank ())
##
## > # 乱子：これはこれでキレイな関係がありそうだけどね
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それでは実際にRで回帰分析をおこなってみましょう
## > # lm() という関数を使います
## >
## > beer.lm <- lm (消費量 ~ 気温, data = beer)
##
## > # 文太：結果を，ここではbeer.lm という名前で保存して
## > # その概要を summary()関数で確認します
## >
## > summary(beer.lm)
##
## Call:
## lm(formula = 消費量 ~ 気温, data = beer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.97642 -0.34722 0.13143 0.41783 0.63307
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.307105 0.361779 6.3771 8.066e-05 ***
## 気温 0.190785 0.019168 9.9531 1.659e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5501 on 10 degrees of freedom
## Multiple R-squared: 0.90831, Adjusted R-squared: 0.89914
## F-statistic: 99.065 on 1 and 10 DF, p-value: 1.6595e-06
##
##
## > # 乱子：「消費量 ~ 気温」ってあるけど，この波みたいな記号の意味が
## > # わからないけど
## > # 文太：それ「チルダ」っていいます
## > # 「消費量 ~ 気温」は，消費量の大きさを気温で回帰する
## > # っていいます
## > # 乱子：回帰するって？
## > # 文太：ようするに消費量を左辺とし，右辺に気温を入れた回帰式を
## > # 求めなさい，って意味です
## > # ちなみに回帰係数だけを取り出して確認したい場合はこうします
## >
## > (coefs <- coef(beer.lm))
## (Intercept) 気温
## 2.30710489 0.19078497
##
## > # 乱子：気温の係数 0.19 と切片の 2.3 ね
## >
## >
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：本書では言及していませんが，もうひとつ回帰分析の例を
## > # 実行してみましょう
## > # Rには車の速度（マイル）と停止までに要した距離（フィート）
## > # のデータがあります
## >
## > cars
## speed dist
## 1 4 2
## 2 4 10
## 3 7 4
## 4 7 22
## 5 8 16
## 6 9 10
## 7 10 18
## 8 10 26
## 9 10 34
## 10 11 17
## 11 11 28
## 12 12 14
## 13 12 20
## 14 12 24
## 15 12 28
## 16 13 26
## 17 13 34
## 18 13 34
## 19 13 46
## 20 14 26
## 21 14 36
## 22 14 60
## 23 14 80
## 24 15 20
## 25 15 26
## 26 15 54
## 27 16 32
## 28 16 40
## 29 17 32
## 30 17 40
## 31 17 50
## 32 18 42
## 33 18 56
## 34 18 76
## 35 18 84
## 36 19 36
## 37 19 46
## 38 19 68
## 39 20 32
## 40 20 48
## 41 20 52
## 42 20 56
## 43 20 64
## 44 22 66
## 45 23 54
## 46 24 70
## 47 24 92
## 48 24 93
## 49 24 120
## 50 25 85
##
## > # 文太：このデータには変数が2つあります
## > # その列名を確認してみます
## >
## > colnames(cars)
## [1] "speed" "dist"
##
## > # 文太：距離とスピードを回帰分析してみましょう
## > # 乱子さん，できそうですか？
## > # 乱子：こんな感じ？
## >
## > lm(dist ~ speed, cars)
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Coefficients:
## (Intercept) speed
## -17.5791 3.9324
##
##
## > # 乱子：Coefficients って？
## > # 文太：係数のことです
## > # 切片が?17.58，speedの係数が3.93ってことですね
## > # 乱子：ビールの消費量と気温とは出力が違うわね
## > # 文太：分析結果をいったん保存して，保存結果にsummary()
## > # を適用すると同じになりますよ
## >
## > cars.lm <-lm(dist ~ speed, cars)
##
## > summary(cars.lm)
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.0691 -9.5253 -2.2719 9.2147 43.2013
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.57909 6.75844 -2.6011 0.01232 *
## speed 3.93241 0.41551 9.4640 1.49e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared: 0.65108, Adjusted R-squared: 0.64381
## F-statistic: 89.567 on 1 and 48 DF, p-value: 1.4898e-12
##
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：気温が20度の場合の予想値を求めてみます
## >
## > 0.19 * 20 + 2.3
## [1] 6.1
##
## > # 文太：ただしRでは以下のように予測値を計算できます
## >
## > beer2 <- data.frame (気温 = c(19, 20, 21))
##
## > predict(beer.lm, newdata = data.frame(気温 = 20))
## 1
## 6.1228043
##
## > # 乱子：関数の名前がpredict()だから，
## > # 予測しようってだけは分かるわ
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ｔ値のグラフを描きます
## >
## > data (tDist)
##
## > p <- ggplot(tDist, aes(x = t値, y = 確率密度)) + geom_line(colour="black", size = 1)
##
## > p + labs(title=" 自由度10のt分布")
##
## > # 乱子：9.953なんて，はるか右のほうね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：最後にビールの実際の消費量と
## > # 回帰モデルが予測した消費量をならべてみましょう
## >
## > beer [, c(5,3)]
## 消費量の予測 消費量
## 1 3.36 2.38
## 2 3.57 3.85
## 3 3.85 4.41
## 4 5.32 5.67
## 5 6.03 5.44
## 6 6.58 6.03
## 7 7.71 8.15
## 8 7.82 8.23
## 9 7.61 7.35
## 10 5.78 5.50
## 11 4.92 4.90
## 12 3.97 4.60
##
## > # 乱子：まあ，だいたいあっているって，感じ？
## > # それよりお腹減ったわね．．．
## >
## > browseURL("http://rmecab.jp/ranko/chap3.html")demo(chap4)##
##
## demo(chap4)
## ---- ~~~~~
## デモを途中でやめる場合はEscキーを押します
##
## > # 文太：第4章は重回帰分析です
## > # 乱子：複回帰分析と呼ぶべきじゃない？
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：説明変数が2つ以上の回帰分析を
## > # 重回帰分析といいます
## > # 売上を気温と曜日から予測する？モデル
## >
## > data(bentoC)
##
## > bento.lm <- lm (売上 ~ 気温 + 曜日, data = bentoC)
##
## > # 文太：「売上 ~ 気温 + 曜日」は
## > # 売上を，気温と曜日から説明する，って意味です
## > # lm 関数は linear model 「線形モデル」の略みたいですね
## > # 乱子：でも，結果が表示されないじゃない
## > # 文太：結果はbento.lmという名前で保存しましたので
## > # サマリーを出します
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > summary (bento.lm)
##
## Call:
## lm(formula = 売上 ~ 気温 + 曜日, data = bentoC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.4209 -5.3360 1.1366 7.7990 19.5337
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.5296 19.0258 1.9726 0.061244 .
## 気温 9.5424 1.4118 6.7589 8.594e-07 ***
## 曜日2火 5.4443 7.8470 0.6938 0.495061
## 曜日3水 12.7038 7.5892 1.6739 0.108302
## 曜日4木 7.5393 7.4016 1.0186 0.319465
## 曜日5金 29.6144 8.7817 3.3723 0.002747 **
## 曜日6土 36.4980 10.1858 3.5832 0.001657 **
## 曜日7日 44.8063 12.6341 3.5464 0.001810 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.424 on 22 degrees of freedom
## Multiple R-squared: 0.93566, Adjusted R-squared: 0.91519
## F-statistic: 45.704 on 7 and 22 DF, p-value: 1.1698e-11
##
##
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：回帰式を使って気温が 16 度で金曜日
## > # の売上個数を予測してみましょう
## >
## > 9.542 * 16.0 + 5.444 * 0 + 12.704 * 0 + 7.539 * 0 + 29.614 * 1 + 36.498 * 0 + 44.806 * 0 + 37.53
## [1] 219.816
##
## > # 乱子：金曜日だけフラグ 1 がたっているのね
## >
## >
## > pause ()
## この行の右端でエンターキーを押して下さい:
##
## > # 文太：実は，Rには指定してモデルから予測行う関数
## > # predictが用意されています
## >
## > predict (bento.lm, newdata = data.frame (気温 = 16, 曜日 = "5金"))
## 1
## 219.82261
##
## > # 文太：わかりにくと思いますが，ようはですね
## > # 新しいデータとして「気温」16度，曜日は5番目の金曜に
## > # してみたってことです
## >
## > browseURL("http://rmecab.jp/ranko/chap4.html")demo(chap5)##
##
## demo(chap5)
## ---- ~~~~~
## デモを途中でやめる場合はEscキーを押します
##
## > # 文太：本章では息抜きに乱子さんがくるくる回る
## > ## アニメgif画像を造りましょう
## > # 乱子：なによそれ？
## > # 文太：詳細はhttp://rmecab.jp/ranko/chap5.html
## > ## を参照してください
## >
## >
## >
## > browseURL("http://rmecab.jp/ranko/chap5.html")
##
## > ## install.packages ("rgl") # rglパッケージのインストール
## >
## > ## library (rgl) # rgl ライブラリを使うという宣言
## >
## > ## open3d() # 3Dグラフィックのウィンドウのみを生成
## >
## > ## とりあえず単純な球体に乱子の画像を貼り付ける
## >
## > ## spheres3d (c(0,0,0),
## > ## texture= "ranko_mini.png", radius=1, color="white",
## > ### alpha=0.8)
## >
## > ## この状態でマウスクリックで画像を主導で回転できます
## > ## ホイールを使うと拡大もできます．
## >
## > ## 自動的に回転させてみる(この場合ImageMagickを
## > ## インストールしている必要はありません)
## > ## if (!rgl.useNULL()) play3d(spin3d(axis=c(0,0,1), rpm=8),
## > ### duration=20 )
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ## C:\temp というフォルダを用意して，アニメgifを作成
## > ### （ImageMagickが必要です）
## >
## > ## Mac の場合は
## > ### Sys.setenv(PATH=paste("/opt/local/bin",
## > ### Sys.getenv("PATH"), sep=":")
## > ### としてImageMagickのconvertプログラムの場所を指定する
## > ### 必要があります
## >
## > ## movie3d( spin3d(), , fps = 16, duration=5 ,
## > ### movie = "C:/temp/anime")
## >
## > ## 作成されたanime.gifはブラウザで開きます．
## > ### 画像をInternetExplorerやFirefoxなどの
## > ### アイコンにドラッグ＆ドロップしてください．demo(chap6)##
##
## demo(chap6)
## ---- ~~~~~
## デモを途中でやめる場合はEscキーを押します
##
## > # 文太：第6章はもりだくさんです
## > # ロジスティック回帰分析，決定木，カイ自乗検定
## > # を説明しますよ！
## >
## > # 乱子：はしゃぎすぎていると，お母さんがやってくるわよ
## > # 文太：うっ．．．
## >
## >
## >
## > ##############################
## > # 文太：それでは乱子さんの失敗モデルを実行してみましょう
## > # 乱子：なんか嫌味たっぷりね！
## > # 文太：いえいえ，だめもとで色々試してみるのは悪いことでは
## > # ないですよ．ええっと，赤牛つけ麺弁当とハッカソンスープを
## > # セットで買うかどうか，でしたね
## >
## >
## > data (aka)
##
## > aka.lm <- lm ( 買上 ~ 気温 + 天気 + 曜日 + IT, data = aka)
##
## > summary (aka.lm)
##
## Call:
## lm(formula = 買上 ~ 気温 + 天気 + 曜日 + IT, data = aka)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.647036 -0.118106 -0.016356 0.200914 0.467415
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.049310 0.787319 -0.0626 0.95072
## 気温 -0.022134 0.057763 -0.3832 0.70584
## 天気2雨 0.207544 0.279056 0.7437 0.46614
## 天気3晴 0.023619 0.287806 0.0821 0.93545
## 曜日2火 0.186393 0.276339 0.6745 0.50811
## 曜日3水 0.626960 0.264270 2.3724 0.02838 *
## 曜日4木 0.191065 0.265925 0.7185 0.48120
## 曜日5金 0.488452 0.301320 1.6210 0.12149
## 曜日6土 0.520510 0.350173 1.4864 0.15357
## 曜日7日 0.469631 0.449737 1.0442 0.30948
## IT 0.530265 0.194301 2.7291 0.01332 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.37334 on 19 degrees of freedom
## Multiple R-squared: 0.5486, Adjusted R-squared: 0.31102
## F-statistic: 2.3091 on 10 and 19 DF, p-value: 0.05589
##
##
## > # 乱子：私が実行したのと同じ結果ね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：このモデルで予測してみると，こんなふうになっちゃいます
## >
## > head (data.frame (日付= 1:30, 予測値 = round (predict (aka.lm), 2)), 10)
## 日付 予測値
## 1 1 0.03
## 2 2 1.00
## 3 3 0.03
## 4 4 0.26
## 5 5 0.74
## 6 6 -0.03
## 7 7 0.60
## 8 8 -0.15
## 9 9 0.56
## 10 10 0.18
##
## > # 乱子：私は，1か0のどっちかになると期待していたんだけど
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ロジスティック回帰分析でも使われる
## > # ロジット曲線ってのはこんなです
## >
## >
## > x <- seq(-4, 4, length.out = 1000)
##
## > DF <- data.frame(ロジット = x , 確率 = plogis(x))
##
## > p <- ggplot(DF, aes(x = ロジット, y = 確率)) + geom_line(colour="black", size = 1 )
##
## > p + labs(title="ロジットと確率の関係") + ylab("確率") + xlab ("ロジットの値")##
## > # 乱子：Y軸が0.0から1.0の間におさまるのがミソなのね
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それでは乱子さんのモデルに
## > # ロジスティック回帰分析を適用してみましょう
## >
## > aka.glm <- suppressWarnings (glm (買上 ~ 気温 + 天気 + 曜日 + IT, data = aka, family=binomial))
##
## > summary (aka.glm)
##
## Call:
## glm(formula = 買上 ~ 気温 + 天気 + 曜日 + IT, family = binomial,
## data = aka)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.635102 -0.000046 0.000000 0.000000 1.635102
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -46.09537 13685.21580 -0.0034 0.9973
## 気温 -1.07833 0.85525 -1.2608 0.2074
## 天気2雨 -1.09710 11381.28385 -0.0001 0.9999
## 天気3晴 -17.32735 8952.18274 -0.0019 0.9985
## 曜日2火 23.66925 15012.91719 0.0016 0.9987
## 曜日3水 61.22836 14765.99410 0.0041 0.9967
## 曜日4木 40.50613 19147.44304 0.0021 0.9983
## 曜日5金 41.56332 13266.33432 0.0031 0.9975
## 曜日6土 43.69431 13266.33470 0.0033 0.9974
## 曜日7日 48.03602 19764.46350 0.0024 0.9981
## IT 39.74752 9561.95593 0.0042 0.9967
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 34.7949 on 29 degrees of freedom
## Residual deviance: 10.1580 on 19 degrees of freedom
## AIC: 32.158
##
## Number of Fisher Scoring iterations: 20
##
##
## > # predict (aka.glm)
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：これも本書には書いていませんが
## > # Rでロジスティック回帰分析は次のように実行します
## > # ここではCOUNTというパッケージに付属する
## > # タイタニック号データを例にとります
## > # 乱子：タイタニック号？
## > # 文太：はい，沈没の際に生き残ることに年齢，性別が影響しているか
## > # どうかを調べるんです
## > ## install.packages("COUNT")
## > ## library (COUNT)
## >
## > data (titanic)
##
## > # 文太；データの最初の部分を確認してみます
## >
## > head (titanic)
## survived age sex class
## 1 1 1 1 1
## 2 1 1 1 1
## 3 1 1 1 1
## 4 1 1 1 1
## 5 1 1 1 1
## 6 1 1 1 1
##
## > # 乱子：survivedって列は二値変数で，生き残ったら 1
## > # 死んでしまったら 0 ってこと？
## > # 文太：はい．他にageは大人なら1，子供なら 0
## > # sexは男性は1，女性は0
## > # classは1ならば1等客室，2は2等客室，3は3等客室です
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ロジスティック回帰分析を実行してみます
## >
## > titanic.glm <- glm(survived ~ age + sex + class, family = binomial,
## + data=titanic)
##
## > # 乱子：glm()関数ってのは？
## > # 文太：一般化線形モデルを実行する関数です
## > # 乱子：一般化線形モデル？
## > # 文太：要するに線形回帰分析やロジスティック回帰分析とか
## > # 汎用的に実行できる関数です
## > # survived ~ age + sex + class は生存したかどうかを，
## > # 年齢と性別，客室等級 で説明するモデルです
## > # 乱子：family=binomialは？
## > # 文太：この部分は，要するに目的変数(survived)が 1 か 0
## > # であることを指示しています
## > # 乱子：で結果は？
## > # 文太：titanic.glm という名前で保存しましたので
## > # summary()を出します
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > summary(titanic.glm)
##
## Call:
## glm(formula = survived ~ age + sex + class, family = binomial,
## data = titanic)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.06517 -0.67181 -0.47400 0.79299 2.11750
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.06188 0.29800 10.2746 < 2.2e-16 ***
## age -1.05561 0.24266 -4.3502 1.360e-05 ***
## sex -2.36946 0.14525 -16.3130 < 2.2e-16 ***
## class2 -1.01056 0.19493 -5.1841 2.171e-07 ***
## class3 -1.76637 0.17072 -10.3468 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1746.76 on 1315 degrees of freedom
## Residual deviance: 1276.20 on 1311 degrees of freedom
## AIC: 1286.2
##
## Number of Fisher Scoring iterations: 4
##
##
## > # 乱子：age sex class のいずれも生存に影響しているってことね
## > # 文太：ええ，実はタイタニック号のデータは女性や子供
## > # そして1等乗客かどうかが生存に影響していることで
## > # 有名なデータなんです
## > # 乱子：要するに女性や子供は率先して救出されたってこと？
## > # 文太：貴族なんかもいた1,2等客室に関してはそうです
## > # 乱子：3等は？
## > # 文太：性別とかあまり関係無いみたいです
## > # 乱子：なんだか，突っ込みどころの多いデータみたいね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ではハッカソンスープと一緒に売れたお弁当の個数を
## > # 確認してみましょう
## >
## > data (hackathon)
##
## > xtabs(~ 性別 + 弁当 , data = hackathon)
## 弁当
## 性別 A B C D E 果 野
## 女 2 1 1 2 1 14 11
##
## > # 乱子：果は果物こんやく丼，野は野菜寿司ね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ハッカソンスープを買ったお客さんがお弁当を選ぶ
## > # 条件を決定木で表してみましょう
## > # 決定木を実行するには，いろいろ外部パッケージを導入します
## > library(rpart)
##
## > #library(rattle)
## > library(rpart.plot)## Warning: package 'rpart.plot' was built under R version 3.6.2
##
## > library(RColorBrewer)
##
## > hackathon.tree <- rpart (弁当~ ., data = hackathon)
##
## > hackathon.tree
## n= 32
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 32 18 果 (0.062 0.031 0.031 0.062 0.031 0.44 0.34)
## 2) 時間=昼 20 11 野 (0.1 0.05 0.05 0.1 0.05 0.2 0.45)
## 4) 年齢< 25 8 5 果 (0.25 0 0.12 0 0 0.38 0.25) *
## 5) 年齢>=25 12 5 野 (0 0.083 0 0.17 0.083 0.083 0.58) *
## 3) 時間=晩 12 2 果 (0 0 0 0 0 0.83 0.17) *
##
## > if (grepl("mac", .Platform$pkgType)){
## + suppressWarnings (prp(hackathon.tree, type=2, extra=101, nn=TRUE, fallen.leaves=TRUE, faclen=0, varlen=0, shadow.col="grey", branch.lty=3, cex = 1.2, split.cex=1.2, under.cex = 1.2, family = "JP1"))
## + } else {
## + suppressWarnings (prp(hackathon.tree, type=2, extra=101, nn=TRUE, fallen.leaves=TRUE, faclen=0, varlen=0, shadow.col="grey", branch.lty=3, cex = 1.2, split.cex=1.2, under.cex = 1.2))
## + }
##
## > # 乱子：グラフを見ると上から，まず時間帯が夜かどうかが重要で
## > # 昼の場合は年齢が25歳未満かどうかが重要と判断されるわけね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：野菜寿司と果物こんやく丼の売上個数（あるいは購買客数）について
## > # 昼と晩それぞれにわけてみましょう．
## >
## >
## > 表1 <- hackathon[hackathon$ 弁当 %in% c("果","野"), ]
##
## > 表1 <- droplevels(表1)
##
## > (表1x <- xtabs(~ 弁当 + 時間, data = 表1))
## 時間
## 弁当 昼 晩
## 果 4 10
## 野 9 2
##
## > # 乱子：分割表ってやつね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：独立性の検定（カイ自乗検定）ではカイ自乗分布という
## > # 検定統計量を使います
## > # カイ自乗値の分布はこんなグラフになります
## >
## > DF <- data.frame(カイ自乗値 = seq(0,10,.1), 確率密度 = dchisq(seq(0,10,.1), 1))
##
## > p <- ggplot(DF, aes(x = カイ自乗値, y = 確率密度)) +
## + geom_line(colour="black", size = 1)
##
## > p + labs(title=" 自由度１のカイ自乗分布")
##
## > # 乱子：正規分布やｔ分布と違って左右対称じゃないのね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：さっきの分割表に行と列それぞれの合計を追加してみます
## >
## > addmargins (表1x)
## 時間
## 弁当 昼 晩 Sum
## 果 4 10 14
## 野 9 2 11
## Sum 13 12 25
##
## > # 乱子：分かりにくいけど Sum ってとこが列ないし行の合計なのね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：列合計と行合計をつかって各セルの割合を求めてみましょう
## >
## > matrix (c(
## + 14/25 * 13/25, 14/25 * 12/25,
## + 11/25 * 13/25, 11/25 * 12/25
## + ), byrow =T, ncol = 2)
## [,1] [,2]
## [1,] 0.2912 0.2688
## [2,] 0.2288 0.2112
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：それぞれの割合がわかれば，それぞれに総数を掛け算
## > # すれば，期待値が求められます
## >
## > (期待値 <- matrix (c( (14/25 * 13/25) * 25, (14/25 * 12/25) * 25, (11/25 * 13/25) * 25, (11/25 * 12/25) * 25), byrow =T, ncol = 2))
## [,1] [,2]
## [1,] 7.28 6.72
## [2,] 5.72 5.28
##
## > # 乱子：割合からしたら，こんな数値になっているはずだっていう
## > # 理想的な値ね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：この期待値と，現実のデータの差異が大きいかどうか
## > # を調べます
## > # そのために期待値と現実の数値を引き算します
## >
## > 表1x - 期待値
## 時間
## 弁当 昼 晩
## 果 -3.28 3.28
## 野 3.28 -3.28
##
## > # 文太：ただし，ここでも自乗します
## >
## > (表1x - 期待値)^2
## 時間
## 弁当 昼 晩
## 果 10.7584 10.7584
## 野 10.7584 10.7584
##
## > # 文太：自乗した結果は期待値で割っておきます
## >
## > (表1x - 期待値)^2 / 期待値
## 時間
## 弁当 昼 晩
## 果 1.4778022 1.6009524
## 野 1.8808392 2.0375758
##
## > # 文太：その合計です．これがカイ自乗統計量です
## >
## > sum ((表1x - 期待値)^2 / 期待値)
## [1] 6.9971695
##
## > # 乱子：まあ，だいたい 7 ってことね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：また怒られそうですが，上の計算は Rではこうやっても
## > # 求められます
## >
## > (xx <- chisq.test(表1x, correct =F))
##
## Pearson's Chi-squared test
##
## data: 表1x
## X-squared = 6.99717, df = 1, p-value = 0.0081639
##
##
## > # 乱子：chisq.test()？
## > # 文太：カイ自乗検定を行うRの関数です
## > # 丸括弧の中に分割表を表す変数を指定すれば
## > # カイ自乗値を計算し，さらに検定結果，つまりは
## > # p値を表示してくれます
## >
## > # 乱子：p-value ってのがp値ね
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：いま検定結果をxxという名前で保存しましたけど
## > # 次のようにすると，期待値を確認できます
## >
## > xx$expected
## 時間
## 弁当 昼 晩
## 果 7.28 6.72
## 野 5.72 5.28
##
## > # 文太；統計検定量（カイ自乗値）をあえて計算してみると
## > # こうですね
## >
## > sum ((表1x - xx$expected)^2 / xx$expected)
## [1] 6.9971695
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：ただ1つ注意があります
## > # 乱子： 注意？
## >
## > # 文太：はい，うえで使ったデータのように2行2列の分割表では
## > # 多くの統計ソフトではイェーツの補正を行います
## > # 乱子：なにそれ
## > # 文太：本書のあとがきにも書いてありますけど
## > # 分析対象のもとの数値は離散値なんですが，検定で使う
## > # カイ自乗値は連続量なんです．そこで，ちょっと調整するんです
## > # 乱子：どうするの？
## > # 文太：いえ，さっきのchisq.test(表x1, correct = F) から
## > # , correct = F を削除するだけです
## >
## > suppressWarnings ( chisq.test(表1x) )
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: 表1x
## X-squared = 5.02648, df = 1, p-value = 0.024963
##
##
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 乱子：p値が変わったけど
## > # 文太：はい．検定の判断そのものが
## > # 代わってしまうこともあります．
## > # 利用しているソフトウェアが
## > # イェーツの補正を行なっているのかどうかを
## > # 必ず確認してください
## >
## >
## > pause()
## この行の右端でエンターキーを押して下さい:
##
## > ##############################
## > # 文太：お弁当と年齢の関係についても
## > # 分割表を作成して，検定してみます
## >
## > 表2 <- xtabs(~ 弁当 + 年齢, data = 表1)
##
## > dimnames (表2) <- list (弁当 = c("果" , "野"), 年齢 = c("25歳未満", "25歳上"))
##
## > 表2
## 年齢
## 弁当 25歳未満 25歳上
## 果 11 3
## 野 4 7
##
## > suppressWarnings ( chisq.test(表2, correct = F) )
##
## Pearson's Chi-squared test
##
## data: 表2
## X-squared = 4.57251, df = 1, p-value = 0.032489
##
##
## > # 乱子：これの帰無仮説は？
## > # 文太：果物こんやく丼を選ぶか野菜寿司を選ぶかに
## > # 年齢は関係ない，つまり年齢とは独立である
## > # 乱子：対立仮説は？
## > # 文太：どっちのお弁当を選ぶかは年齢が影響する
## > # つまり年齢は独立ではない
## > # 乱子：カイ自乗値が4.57でp値が0.03だから帰無仮説は棄却されるのね
## > # 文太：はい．つまり2つのお弁当の選択には年齢が影響していることに
## > # なります．
## >
## > # 乱子：あれ？お母さんが来た
## > # 文太：え？
## >
## >
## >
## > browseURL("http://rmecab.jp/ranko/chap6.html")