Homework#2
Katsuya Ogawa
2022-01-22
#load library, import data, arrangement data
library(SuperLearner)## Warning: パッケージ 'SuperLearner' はバージョン 4.1.2 の R の下で造られました## 要求されたパッケージ nnls をロード中です## 要求されたパッケージ gam をロード中です## Warning: パッケージ 'gam' はバージョン 4.1.2 の R の下で造られました## 要求されたパッケージ splines をロード中です## 要求されたパッケージ foreach をロード中です## Warning: パッケージ 'foreach' はバージョン 4.1.2 の R の下で造られました## Loaded gam 1.20## Super Learner## Version: 2.0-28## Package created on 2021-05-04library(estimatr)## Warning: パッケージ 'estimatr' はバージョン 4.1.2 の R の下で造られましたlibrary(tidyverse)## Warning: パッケージ 'tidyverse' はバージョン 4.1.2 の R の下で造られました## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.1 v forcats 0.5.1## Warning: パッケージ 'ggplot2' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'tibble' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'tidyr' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'readr' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'purrr' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'dplyr' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'stringr' はバージョン 4.1.2 の R の下で造られました## Warning: パッケージ 'forcats' はバージョン 4.1.2 の R の下で造られました## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x purrr::accumulate() masks foreach::accumulate()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x purrr::when() masks foreach::when()library(grf)## Warning: パッケージ 'grf' はバージョン 4.1.2 の R の下で造られましたdata <- read_csv("data(2).csv")## Rows: 2129 Columns: 16## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (7): region, training, health, education, firm_size, marriage, gender
## dbl (9): income, weights, industry, occupation, num_child, working_hour, age...##
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.Y <- data$working_hour
D <- if_else(data$gender == "Female",1, 0)
X <- select(data, age, exp)この課題では予測する値を労働時間、関心のあるパラメータを性別、他のパラメータを年齢、経験年数とした。
#modeling
fit.Y <- CV.SuperLearner(X = X,
Y = Y,
SL.library = c("SL.lm",#ols
"SL.glmnet",#LASSO
"SL.ranger")#random forest
)## 要求されたパッケージ glmnet をロード中です## 要求されたパッケージ ranger をロード中ですfit.D <- CV.SuperLearner(X = X,
Y = D,
SL.library = c("SL.lm",
"SL.glmnet",
"SL.ranger")
)
summary(fit.Y)##
## Call:
## CV.SuperLearner(Y = Y, X = X, SL.library = c("SL.lm", "SL.glmnet", "SL.ranger"))
##
##
## Risk is based on: Mean Squared Error
##
## All risk estimates are based on V = 10
##
## Algorithm Ave se Min Max
## Super Learner 87.928 4.9610 52.255 113.22
## Discrete SL 87.879 4.9644 52.049 113.34
## SL.lm_All 87.922 4.9643 52.062 113.84
## SL.glmnet_All 87.905 4.9676 52.049 113.34
## SL.ranger_All 93.780 4.9155 61.755 119.40summary(fit.D)##
## Call:
## CV.SuperLearner(Y = D, X = X, SL.library = c("SL.lm", "SL.glmnet", "SL.ranger"))
##
##
## Risk is based on: Mean Squared Error
##
## All risk estimates are based on V = 10
##
## Algorithm Ave se Min Max
## Super Learner 0.21287 0.0035678 0.19764 0.22445
## Discrete SL 0.21601 0.0035132 0.19653 0.22749
## SL.lm_All 0.21613 0.0035357 0.19639 0.22776
## SL.glmnet_All 0.21587 0.0035070 0.19653 0.22749
## SL.ranger_All 0.22150 0.0044559 0.20036 0.23911fit.Yで結果変数の予測モデルを作る。このモデルではolsの値が最も小さいためパフォーマンスが一番優れていると考えられる。
同様にfit.Dで原因変数の予測モデルを作る。最もパフォーマンスが良いのはSuper Learnerである。
#parameter estimation
X.matrix <- model.matrix(~ age + exp ,
data = X)
X.matrix <- X.matrix[,-1]
fit.tau <- causal_forest(X = X.matrix,
Y = Y,
W = D,
Y.hat = fit.Y$SL.predict,
W.hat = fit.D$SL.predict)
hist(fit.tau$predictions)
best_linear_projection(fit.tau,scale(X.matrix))##
## Best linear projection of the conditional average treatment effect.
## Confidence intervals are cluster- and heteroskedasticity-robust (HC3):
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.6215199 0.4116030 -11.2281 < 2.2e-16 ***
## age -0.0013235 0.5078944 -0.0026 0.997921
## exp 1.5748800 0.5134204 3.0674 0.002186 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ここでは効果の異質性を推定する。性別のパラメータtauがどれだけ労働時間に影響を与えるかをcausal_forestをつかって推定する。
ヒストグラムを見ると、因果効果の値であるtauが-6から-4のサンプル数が多くなっているが、-5あたりを中心として裾野の広い分布をしているため、因果効果の異質性が高いと考えられる。つまり、女性であることによる労働時間の変化は-6時間から-4時間が最も多いが、ばらつきは大きいという解釈ができる。
また、層別の効果の推定にはbest_linear_projection を使う。この推定の結果によると経験年数が大きいほうが性別による効果の異質性が大きく、p値も0.0027と小さい値を出しているため信頼できる値であると考えられる。