# encording with Shift-JIS

#
# 8章 ダミー変数
#

# C:\EconDat\8_dummy
setwd("C:/EconDat/8_dummy")
getwd()

#
# 例8.1:期末試験得点と性別
#

df_test <- read.csv("../data/testscore/testscore.csv",header = TRUE)
head(df_test)

## 方法1：
# df_test に 女性なら1、男性なら0をとるダミー変数「female」を追加する
df_test[ which(df_test$sex=="女性"), "female" ] <- 1
df_test[ which(df_test$sex=="男性"), "female" ] <- 0
head(df_test)

#「score」列と「female」列を使って単回帰
res1 <- lm(score ~ female, df_test)
summary(res1)

## 方法2：
# 「sex」列のカテゴリー定義順序を確認する
levels( df_test$sex )

# 「sex」列の基準グループを男性に変更
df_test$sex <- relevel(df_test$sex, ref = "男性")
levels(df_test$sex) # 「男性」が最初になったか確認

#「score」列と「sex」列を使って単回帰
res2 <- lm(score ~ sex, df_test)
summary(res2)

#
# 例8.2: 期末試験得点と出席回数・性別
#

res3<- lm(score ~ female+attend_num, df_test)
summary(res3)

#
# 例8.3: 基準グループを設けない回帰モデル
#

# df_test に 男性なら1、女性なら0をとるダミー変数「male」を追加する
df_test[ which(df_test$sex=="女性"), "male" ] <- 0
df_test[ which(df_test$sex=="男性"), "male" ] <- 1
head(df_test)

# 女性ダミー変数と男性ダミー変数の両方を含む回帰モデル
# 定数項も含むモデル
res4<- lm(score ~ female +male +attend_num, df_test)
summary(res4)

# 定数項を除くモデル
res5<- lm(score ~ female +male +attend_num -1, df_test)
summary(res5)

#
# 例8.4:出席回数の期末試験得点への影響に男女差はあるか
#

res6 <- lm(score~ female+attend_num + female:attend_num, df_test)
summary(res6)

#
# 図8.3:（8.20）式が示す回帰直線
#

# グラフ
plot(df_test$attend_num[df_test$female==0],df_test$score[df_test$female==0],col="blue", ylim = c(0,100),xlim = c(0,30),
     xlab="attend_num",ylab = "score")
points(df_test$attend_num[df_test$female==1],df_test$score[df_test$female==1],col="orange", pch=16)
legend(1,95,legend = c("male","female"),col = c("blue","orange"), pch = c(1,16),bty = "n")
abline(a=res6$coefficients[1],b=res6$coefficients[3], col="blue",lwd=2)
abline(a=res6$coefficients[1]+res6$coefficients[2],b=res6$coefficients[3]+res6$coefficients[4], 
       col="orange",lwd=2)

#
# 例8.5: 女性ダミー変数、留学生ダミー変数を含む場合
#

# 留学生ダミー変数の作成
df_test[which(df_test$if_international=="はい"), "international_st"] <-1
df_test[which(df_test$if_international=="いいえ"), "international_st"] <-0
head(df_test)

#回帰式
res7 <- lm(score~ female +international_st +attend_num, df_test)
summary(res7)

#
# 例8.6: 女性ダミー変数と留学生ダミー変数の交差項を含む
#

#回帰式
res8 <- lm(score~ female +international_st +female:international_st +attend_num, df_test)
summary(res8)

#
# 例8.7: 学年の期末試験得点への影響
#

#データの構造
str(df_test)

#回帰式
res9 <- lm(score~ school_year + +attend_num, df_test)
summary(res9)

#
# 例8.8: 学年ダミー変数を用いた期末試験得点の分析
#

# school_yearは数値データなので、factor関数でカテゴリー化し
# levels関数で基準グループを確認する。
levels(factor(df_test$school_year))

res10 <- lm(score~ factor(school_year) +attend_num, df_test)
summary(res10)

# write.csv関数を用いた回帰分析の出力結果の保存
write.csv( summary(res10)$coef, file = "score_res10.csv")

#
# 例8.9: 学生アルバイトの線形確率モデル
#

# 学生生活アンケート2018
df_survey <- read.csv("../data/student_survey2018/student_survey2018_data.csv", header=TRUE)

## アルバイトをしているか
part_job <- df_survey$Q42

# 性別
sex <- df_survey$Q2

## 留学生かどうか
international_st <- df_survey$Q6

## 学年
school_year <- df_survey$Q4
summary(school_year)

## データの変更
school_year[school_year=='5年'] <-NA
school_year[school_year=='修士課程2年'] <-NA
school_year[school_year=='博士後期課程3年'] <-NA
school_year[school_year=='教員'] <-NA
# school_yearを4つのレベルに再定義
school_year <- factor(school_year)
summary(school_year)

##　アルバイト分析用データフレーム
df_part_job <- data.frame(df_survey$id, part_job, sex, international_st, school_year)
head(df_part_job)

#方法2による分析

##分析用の変数作成

##アルバイトダミー変数
df_part_job[which(df_part_job$part_job=='している'), "part_job_d"]<-1
df_part_job[which(df_part_job$part_job=='していない'), "part_job_d"]<-0
summary(df_part_job$part_job_d)

## level関数によるデータの確認
levels(df_part_job$sex)
levels(df_part_job$school_year)
levels(df_part_job$international_st)

## relevel international_st 
df_part_job$international_st<- relevel(df_part_job$international_st, ref="留学生ではない。")
levels(df_part_job$international_st)
str(df_part_job)

## 回帰分析
part_job_result <- lm(part_job_d ~ sex +international_st +school_year , df_part_job)
summary(part_job_result)


# 参考: 方法1
## ダミー変数を作成し利用した回帰分析

## アルバイトダミー変数 part_time job dummy
df_survey[which(df_survey$Q42=='している'), "part_job_d"]<-1
df_survey[which(df_survey$Q42=='していない'), "part_job_d"]<-0

## 男性ダミー変数male dummy
df_survey[ which(df_survey$Q2=="女性"), "male_d" ] <- 0
df_survey[ which(df_survey$Q2=="男性"), "male_d" ] <- 1

subset(df_survey,  c(Q2, male_d))  

## 女性ダミー変数 female dummy
df_survey[ which(df_survey$Q2=="女性"), "female_d" ] <- 1
df_survey[ which(df_survey$Q2=="男性"), "female_d" ] <- 0

## 留学生ダミー変数
df_survey[ which(df_survey$Q6=='留学生ではない。'), "international_st_d"]<-0
df_survey[ which(df_survey$Q6=='留学生である。'), "international_st_d"]<-1


## アルバイトとfemale dammy, international student dummy, grade
part_job_result1 <- lm(part_job_d ~ male_d+international_st_d+school_year, df_survey)
summary(part_job_result1)

#
# 付録 Appendix
# 8A 8章 ダミー変数を用いた回帰分析・付論: 2項選択モデル
#

#
# 例8.10: アルバイト就業の線形確率モデル
#

# データの読み込みとデータフレームの作成
df_survey <- read.csv("../data/student_survey2018/student_survey2018_data.csv", header=TRUE)
head(df_survey) 　


## アルバイトをしているか
part_job <- df_survey$Q42

# 性別
sex <- df_survey$Q2

## 留学生かどうか
international_st <- df_survey$Q6

## 学年
school_year <- df_survey$Q4
summary(school_year)

## データの変更
school_year[school_year=='5年'] <-NA
school_year[school_year=='修士課程2年'] <-NA
school_year[school_year=='博士後期課程3年'] <-NA
school_year[school_year=='教員'] <-NA
school_year <- factor(school_year) #school_year を四つのlevels に再定義
summary(school_year)

##　アルバイト分析用データフレーム
df_part_job <- data.frame(df_survey$id, part_job, sex, international_st, school_year)
df_part_job <- na.omit(df_part_job)
head(df_part_job)

#方法2による分析

##分析用の変数作成

##アルバイトダミー変数
df_part_job[which(df_part_job$part_job=='している'), "part_job_d"]<-1
df_part_job[which(df_part_job$part_job=='していない'), "part_job_d"]<-0
summary(df_part_job$part_job_d)

## level関数によるデータの確認
levels(df_part_job$sex)
levels(df_part_job$school_year)
levels(df_part_job$international_st)

## relevel international_st 
df_part_job$international_st<- relevel(df_part_job$international_st, ref="留学生ではない。")
levels(df_part_job$international_st)


## 回帰分析
part_job_result <- lm(part_job_d ~ sex +international_st +school_year , df_part_job)
summary(part_job_result)


#残差の保存
part_job_resid <- residuals(part_job_result)
#予測値の保存
part_job_fit <- fitted(part_job_result)


# 実際と予測値のクロス集計表
fit_order <- cut(part_job_fit, breaks=c(-0.5,0, 0.5,1.0, 1.5), right= FALSE, include.lowest = TRUE)
actual_fit_ctab <- table(df_part_job$part_job_d, fit_order)
actual_fit_ctab_m <- addmargins(actual_fit_ctab)
actual_fit_ctab_m

#
# 例8.11：アルバイト就業のロジット・モデルとプロビット・モデル
#

## ロジットモデルの推定
part_job_d_3_logit <- glm(part_job_d ~ sex +international_st +school_year, family=binomial(link = "logit"), data=df_part_job)
summary(part_job_d_3_logit)

#残差の保存
part_job_d_3_logit_resid <- residuals(part_job_d_3_logit)
#予測値の保存
part_job_d_3_logit_fit <- fitted(part_job_d_3_logit)
summary(part_job_d_3_logit_fit)

# 実際と予測値のクロス集計表
fit_order <- cut(part_job_d_3_logit_fit, breaks=c(-0.5, 0,0.5,1.0, 1.5), right= FALSE, include.lowest = TRUE)
#breaks=c(0.1,0.2,0.3, 0.4,0.5,0.6, 0.7, 0.8,0.9,1.0)
actual_fit_ctab <- table(df_part_job$part_job_d, fit_order)
actual_fit_ctab_m <- addmargins(actual_fit_ctab)
actual_fit_ctab_m

# 比率
actual_fit_ctab_p <- prop.table(actual_fit_ctab)
actual_fit_ctab_pm <- addmargins(actual_fit_ctab_p)
actual_fit_ctab_pm

## プロビット
part_job_d_3_probit <- glm(part_job_d ~ sex +international_st +school_year, family=binomial(link = "probit"), data=df_part_job)
summary(part_job_d_3_probit)

#残差の保存
part_job_d_3_probit_resid <- residuals(part_job_d_3_probit)
#予測値の保存
part_job_d_3_probit_fit <- fitted(part_job_d_3_probit)
summary(part_job_d_3_probit_fit)

# 予測値と実際のクロス集計表
fit_order <- cut(part_job_d_3_probit_fit,  breaks=c(-0.5,0,0.5,1.0, 1.5), right= FALSE, include.lowest = TRUE)
#breaks=c(0.1,0.2,0.3, 0.4,0.5,0.6, 0.7, 0.8,0.9,1.0)
actual_fit_ctab <- table(df_part_job$part_job_d, fit_order)
actual_fit_ctab_m <- addmargins(actual_fit_ctab)
actual_fit_ctab_m

#
# 例8.12: 不均一分散一致標準誤差を利用したアルバイト就業の確率線形モデル
#

#不均一分散一致標準誤差を求めるためのパッケージ

#パッケージのインストール< 1度実行すればよい>
# install.packages("estimatr")

#パッケージの読み込み
library(estimatr)

#不均一分散一致推定量
library(estimatr)
part_job_d_rob <- lm_robust(part_job_d ~ sex +international_st +school_year ,se_type="HC0",  df_part_job)
summary(part_job_d_rob)


