S&W 第8章サブゼミ
Risako Morita
2018/11/19
1. R Markdownの使い方
基本事項
● R MarkdownとはR（またはR Studio）でドキュメントが作れるという機能のこと
● コードを簡単に再現できるという利点。共有も簡単。
● 実際の流れ
1. RStudioを起動して、画面左上のボタンからR Markdownを選択
2. TitleやAuthorを入力、Output FormatはHTMLがおすすめ
3. ファイルが自動で生成されるので、内容を編集する
4. Knitする。これはレンダリングを実行するコマンドのこと。レンダリングとは：データ記述言語やデータ構造で記述された抽象的で高次の情報から、コンピュータのプログラムを用いて画像・映像・音声などを生成すること（Wikipedia参照）レンダリングが終了すると、自動的に.htmlファイルが作成され、Viewerに表示される。
5.それを確認して、不備がないなら右上のPublishを選択
R Markdown記法
（以下は、四角囲みにコードを、その下に実行例を併記している）
1. 見出し
#レベル１
##レベル2
###レベル3レベル１
レベル２
レベル３
2. 改行と空白行の挿入
文章の最後に半角スペースを2つ以上入れた後にエンターキーで改行することで、強制改行できる。<br/>で空白行が挿入できる。
3. 箇条書き
*の後は必ず半角スペースを入れる。行頭に半角スペースを4つ入れると、記号のレベルが下がる。
* 箇条書き
* 箇条書き
* 箇条書き- 箇条書き
- 箇条書き
- 箇条書き
- 箇条書き
4. 番号付き箇条書き
番号とピリオドの後に半角スペースを入れる。行頭に半角スペースを4つ入れると、レベルが下がる。
1. 箇条書き
1. 箇条書き
2. 箇条書き
1. 箇条書き- 箇条書き
- 箇条書き
- 箇条書き
- 箇条書き
5. 強調文字・斜体文字
**強調文字**
*斜体文字* 強調文字
斜体文字
6. 水平線
3つ以上の*または-で水平線
***
---
7. コードブロック
`この記号を３つ入力したもので内容を挟むと、入力した内容がそのまま表示される。この章の四角囲みもこのコードブロックを使って書いている。（ちなみにコードブロックのコードをコードブロックしたい時は````で挟むといいみたい）
```
（ここに内容を記述）
```
8. Rチャンク（Rのコードを記述したブロック）
先ほどのコードブロックの中にRのコードを書いて、Rで実行した結果とともに表示させる。
```{r}
head(iris)
```こう入力すると、Rのコード（グレーの四角囲み）と、その実行結果が併記されて下のように表示される。
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
9. インラインコード
Rチャンクを作るほどでもない（文章中に簡単な計算結果を埋め込みたいなどの）場合は、 `とrと半角スペースとRのコードと｀
をこの順に入力すれば良い
この相関はr cor(iris[[1]],iris[[2]])`です(はじめのrの前に｀を追加！)この相関は-0.1175698です
R Markdowmの記法はRStudioのHelpバーのMarkdown Quick Referenceにも解説がある。（ただし英語）
以上、参考にしたサイトはこちら：https://kazutan.github.io/kazutanR/Rmd_intro.html#markdown記法
2. S&W第8章
準備
データファイルの読み込み
CPS04 <-read.csv("CPS04.csv")
TeachingRatings <-read.csv("TeachingRatings.csv")
CollegeDistance <-read.csv("CollegeDistance.csv")
Growth <-read.csv("Growth_nm.csv")E8.1
a
res1_1 <- lm(ahe ~age + female + bachelor,
data=CPS04)
summary(res1_1)##
## Call:
## lm(formula = ahe ~ age + female + bachelor, data = CPS04)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.477 -5.098 -1.266 3.431 44.828
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.88380 0.92029 2.047 0.0407 *
## age 0.43920 0.03053 14.387 <2e-16 ***
## female -3.15786 0.18036 -17.508 <2e-16 ***
## bachelor 6.86515 0.17837 38.489 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.884 on 7982 degrees of freedom
## Multiple R-squared: 0.19, Adjusted R-squared: 0.1897
## F-statistic: 624.1 on 3 and 7982 DF, p-value: < 2.2e-16母集団の回帰式は線形であるから、Ageが25歳から26歳になる時も、33歳から34歳になる時も、賃金の増加額は同じ。その金額は、ageの係数を見れば良いから、0.4392042である。つまり時給0.4ドルアップ （参考：p229） ちなみに、この数字を出力するコードはsummary(res1_1)$coefficients[2,1]
しかし、賃金の上昇について分析するにあたって、単位を「ドル」で測ることは適切なのか？「パーセント」で評価する方が適切なのではないか。
理由（参考:p236）
1. そもそもの賃金水準が異なる業種間（または時代間）での比較が可能になる
2. 対数を取ってしまえば、線形の関係になるかもしれない
3. 弾力性についての分析と親和性が高い
対数回帰モデルの復習（板書1参照）
b
まずは賃金の対数を取って、ln_aheという新たな変数を作る
CPS04$ln_ahe <- log(CPS04$ahe)次に回帰
res1_2 <- lm(ln_ahe ~ age + female + bachelor,
data = CPS04 )
summary(res1_2)##
## Call:
## lm(formula = ln_ahe ~ age + female + bachelor, data = CPS04)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1926 -0.2773 0.0141 0.2936 1.5876
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.85646 0.05335 34.80 <2e-16 ***
## age 0.02444 0.00177 13.81 <2e-16 ***
## female -0.18046 0.01046 -17.26 <2e-16 ***
## bachelor 0.40527 0.01034 39.19 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4571 on 7982 degrees of freedom
## Multiple R-squared: 0.1924, Adjusted R-squared: 0.1921
## F-statistic: 633.8 on 3 and 7982 DF, p-value: < 2.2e-16Yは対数を取り、Xは対数を取らないモデル（対数・線形モデル）である。年齢が1歳増加すると、ageの係数に100をかけた%だけ賃金が増加する。つまり25歳から26歳に変化する時も、33歳から34歳に変化する時も、2.4442948%だけ賃金が増加する。
参考：この数字のコードはsummary(res1_2)$coefficients[2,1]*100
c
CPS04$ln_age <- log(CPS04$age)
res1_3 <- lm(ln_ahe ~ ln_age + female + bachelor ,
data = CPS04)
summary(res1_3)##
## Call:
## lm(formula = ln_ahe ~ ln_age + female + bachelor, data = CPS04)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.19384 -0.27859 0.01487 0.29455 1.59413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.12828 0.17662 0.726 0.468
## ln_age 0.72470 0.05204 13.925 <2e-16 ***
## female -0.18030 0.01045 -17.245 <2e-16 ***
## bachelor 0.40523 0.01034 39.195 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.457 on 7982 degrees of freedom
## Multiple R-squared: 0.1927, Adjusted R-squared: 0.1924
## F-statistic: 635 on 3 and 7982 DF, p-value: < 2.2e-16これは対数・対数モデルなので、ln_ageの係数は、年齢に関する賃金の弾力性に値する。年齢が1%上昇すると、賃金は0.7246973%上昇する。つまり25歳から26歳（4％の上昇）は2.8987893%の上昇、33歳から34歳（100/33%の上昇）は2.1960525%の上昇
d
CPS04$age2 <- CPS04$age^2
res1_4 <- lm(ln_ahe ~ age + age2 + female + bachelor,
data = CPS04)
summary(res1_4)##
## Call:
## lm(formula = ln_ahe ~ age + age2 + female + bachelor, data = CPS04)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.19772 -0.27912 0.01245 0.29594 1.61436
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0587335 0.6104789 0.096 0.923357
## age 0.1470452 0.0415124 3.542 0.000399 ***
## age2 -0.0020706 0.0007004 -2.956 0.003125 **
## female -0.1797868 0.0104542 -17.198 < 2e-16 ***
## bachelor 0.4050769 0.0103362 39.190 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4569 on 7981 degrees of freedom
## Multiple R-squared: 0.1933, Adjusted R-squared: 0.1929
## F-statistic: 478 on 4 and 7981 DF, p-value: < 2.2e-16これは2次の回帰モデルであるので、年齢の増加の影響は、変化前の年齢に依存する。
・25歳から26歳
age=25の時のln_aheとage=26の時のln_aheが求まる。
summary(res1_4)$coefficients[2,1]*25
+ summary(res1_4)$coefficients[3,1]*25^2
summary(res1_4)$coefficients[2,1]*26
+ summary(res1_4)$coefficients[3,1]*26^2を計算すると、2.382029 と 2.4234756
・33歳から34歳
summary(res1_4)$coefficients[2,1]*33
+ summary(res1_4)$coefficients[3,1]*33^2
summary(res1_4)$coefficients[2,1]*34
+ summary(res1_4)$coefficients[3,1]*34^2を計算すると2.5976501 と2.6059677
ln_aheの予測値は計算できるが、aheの原数値は直接は求めることが難しい（p243参照)
あえて計算するなら上の対数値の指数を取って
10.8268481 と 11.2850129
13.4321374 と13.5443264
e
bとcの比較：対数・線形モデルと対数・対数モデルの比較には修正済み決定係数の比較が可能である。（参考p241）bは0.1921、cは0.1924であるので、対数・対数モデル(c)の方が当てはまりが良い。
f
bとdの比較：1次の回帰モデルと2次の回帰モデルの比較。もし真の関係が線形であるならば、2次の回帰式のβ_2の係数はゼロである。そこでtテストを行う（参考p229）
β_2=0の帰無仮説。t統計量は
summary(res1_4)$coefficients[3.1]/summary(res1_4)$coefficients[3.2]んー計算できないので手動
-0.0020706/0.0007004-2.9563107 これは絶対値で5％水準の臨界値（1.96)を上回るので、棄却。つまり2次モデル（d）の妥当性が裏付けられた。
g
cとdの比較：ともにln_aheを回帰しているから修正済み決定係数の比較が可能。cは0.1924でdは0.1929であるから、dの方が当てはまりが良いと言える。
h
まずは高卒男性のみを抽出したデータフレームを作る
CPS04_nbm <- subset(CPS04, bachelor == 0 & female == 0)subset関数の注意点：イコールは==（必ず２つ）。条件式で「かつ」を表すのは&
non-bachelor maleの略でnbmです…
ageとln_aheの関係をプロットする
plot(ln_ahe ~ age ,data = CPS04_nbm)
b,c,dの回帰をやり直す（コピペしてnbmを書き加える）
res1_2_nbm <- lm(ln_ahe ~ age + female + bachelor,
data = CPS04_nbm)
res1_3_nbm <- lm(ln_ahe ~ ln_age + female + bachelor ,
data = CPS04_nbm)
res1_4_nbm <- lm(ln_ahe ~ age + age2 + female + bachelor,
data = CPS04_nbm)res1_2_nbm <- lm(ln_ahe ~ age + female + bachelor,
data = CPS04_nbm)
res1_3_nbm <- lm(ln_ahe ~ ln_age + female + bachelor ,
data = CPS04_nbm)
res1_4_nbm <- lm(ln_ahe ~ age + age2 + female + bachelor,
data = CPS04_nbm)summary(res1_2_nbm)##
## Call:
## lm(formula = ln_ahe ~ age + female + bachelor, data = CPS04_nbm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8404 -0.2868 0.0170 0.2952 1.5583
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.982470 0.091605 21.641 < 2e-16 ***
## age 0.020573 0.003069 6.704 2.44e-11 ***
## female NA NA NA NA
## bachelor NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4635 on 2770 degrees of freedom
## Multiple R-squared: 0.01597, Adjusted R-squared: 0.01561
## F-statistic: 44.95 on 1 and 2770 DF, p-value: 2.443e-11summary(res1_3_nbm)##
## Call:
## lm(formula = ln_ahe ~ ln_age + female + bachelor, data = CPS04_nbm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.84320 -0.28698 0.01663 0.29589 1.56384
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.52790 0.30574 1.727 0.0844 .
## ln_age 0.60996 0.09023 6.760 1.68e-11 ***
## female NA NA NA NA
## bachelor NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4634 on 2770 degrees of freedom
## Multiple R-squared: 0.01623, Adjusted R-squared: 0.01587
## F-statistic: 45.69 on 1 and 2770 DF, p-value: 1.68e-11summary(res1_4_nbm)##
## Call:
## lm(formula = ln_ahe ~ age + age2 + female + bachelor, data = CPS04_nbm)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85429 -0.28855 0.01657 0.29984 1.58068
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.487209 1.046617 0.466 0.6416
## age 0.122542 0.071167 1.722 0.0852 .
## age2 -0.001722 0.001201 -1.434 0.1516
## female NA NA NA NA
## bachelor NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4634 on 2769 degrees of freedom
## Multiple R-squared: 0.0167, Adjusted R-squared: 0.01599
## F-statistic: 23.51 on 2 and 2769 DF, p-value: 7.5e-11そして先ほどの散布図に回帰直線を引くんですが、いつも通りabline関数を使えるのはres1_2_nbmの回帰直線のみ。
plot(ln_ahe ~ age ,data = CPS04_nbm)
abline(res1_2_nbm)## Warning in abline(res1_2_nbm): only using the first two of 4 regression
## coefficients
残り二つの非線形回帰を散布図に重ねるにはどうしたらいいのか???わからない
とりあえず大卒女性に対して同じ操作を行う。
CPS04_bf <- subset(CPS04, bachelor == 1 & female == 1)
plot(ln_ahe ~ age ,data = CPS04_bf)
res1_2_bf <- lm(ln_ahe ~ age + female + bachelor,
data = CPS04_bf)
res1_3_bf <- lm(ln_ahe ~ ln_age + female + bachelor ,
data = CPS04_bf)
res1_4_bf <- lm(ln_ahe ~ age + age2 + female + bachelor,
data = CPS04_bf)
summary(res1_2_bf)##
## Call:
## lm(formula = ln_ahe ~ age + female + bachelor, data = CPS04_bf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.90274 -0.24867 0.02397 0.27660 1.15694
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.137898 0.109219 19.574 < 2e-16 ***
## age 0.023119 0.003683 6.277 4.34e-10 ***
## female NA NA NA NA
## bachelor NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4461 on 1737 degrees of freedom
## Multiple R-squared: 0.02218, Adjusted R-squared: 0.02162
## F-statistic: 39.41 on 1 and 1737 DF, p-value: 4.336e-10summary(res1_3_bf)##
## Call:
## lm(formula = ln_ahe ~ ln_age + female + bachelor, data = CPS04_bf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.90110 -0.24949 0.02554 0.27408 1.16212
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5219 0.3649 1.43 0.153
## ln_age 0.6800 0.1079 6.30 3.76e-10 ***
## female NA NA NA NA
## bachelor NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4461 on 1737 degrees of freedom
## Multiple R-squared: 0.02234, Adjusted R-squared: 0.02178
## F-statistic: 39.69 on 1 and 1737 DF, p-value: 3.756e-10summary(res1_4_bf)##
## Call:
## lm(formula = ln_ahe ~ age + age2 + female + bachelor, data = CPS04_bf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.89924 -0.24444 0.02399 0.27426 1.16769
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3561963 1.2514099 1.084 0.279
## age 0.0766237 0.0854075 0.897 0.370
## age2 -0.0009067 0.0014460 -0.627 0.531
## female NA NA NA NA
## bachelor NA NA NA NA
##
## Residual standard error: 0.4462 on 1736 degrees of freedom
## Multiple R-squared: 0.0224, Adjusted R-squared: 0.02128
## F-statistic: 19.89 on 2 and 1736 DF, p-value: 2.871e-09plot(ln_ahe ~ age ,data = CPS04_bf)
abline(res1_2_bf)## Warning in abline(res1_2_bf): only using the first two of 4 regression
## coefficients
lines関数を使ってみたけど、絶対違う…
plot(ln_ahe ~ age ,data = CPS04_bf)
lines(CPS04_bf$age, summary(res1_3_bf)$coef[1,1] + summary(res1_3_bf)$coef[2,1]*log(CPS04_bf$age))
plot(ln_ahe ~ age ,data = CPS04_bf)
lines(CPS04_bf$age, summary(res1_4_bf)$coef[1,1] + summary(res1_4_bf)$coef[2,1]*(CPS04_bf$age) + summary(res1_4_bf)$coef[3,1]*(CPS04_bf$age2))
参考： 綺麗なグラフを描く
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.4.4library(tidyverse)## ─ Attaching packages ───────────────── tidyverse 1.2.1 ─## ✔ tibble 1.4.2 ✔ purrr 0.2.5
## ✔ tidyr 0.8.1 ✔ dplyr 0.7.7
## ✔ readr 1.1.1 ✔ stringr 1.3.1
## ✔ tibble 1.4.2 ✔ forcats 0.3.0## Warning: package 'tibble' was built under R version 3.4.3## Warning: package 'tidyr' was built under R version 3.4.4## Warning: package 'purrr' was built under R version 3.4.4## Warning: package 'dplyr' was built under R version 3.4.4## Warning: package 'stringr' was built under R version 3.4.4## Warning: package 'forcats' was built under R version 3.4.3## ─ Conflicts ─────────────────── tidyverse_conflicts() ─
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()ggplot(CPS04_bf,aes(x=age, y=ln_ahe)) + geom_point() + geom_smooth(method = "lm")