pacman::p_load(tidyverse, ggplot2, car)
dta <- read.csv("data0505.csv", header = T)
options(digits = 3)
dta <- dta %>% mutate( Gender = relevel(Gender, ref = "女"),
Sector = relevel(Sector, ref = "私立"),
Field = relevel(Field, ref = "大眾傳播學群"),
EduLv = factor(EduLv, levels=c("博士","碩士","普通大學","科技大學",
"技術學院","五專","三專",
"二專","高中","高職","國中")),
EduLv = relevel(EduLv, ref = "技術學院"),
Region = factor(Region, levels =c("宜花東離島","北北基","桃竹苗",
"中彰投","雲嘉南","高屏澎")),
Age = as.numeric(Age),
J_year = as.numeric(J_year),
JobZone = as.numeric(JobZone),
EduZone = as.numeric(EduZone),
JobZone_D = as.numeric(EduZone-JobZone),
Salary = as.numeric(Salary),
SubEduOver = relevel(SubEduOver, ref="符合工作要求"),
Core = recode_factor(as.factor(JobCor), "1" = "無關聯",
"2" = "部分關聯",
"3" = "核心關聯"),
SubEduOver = factor(SubEduOver,levels =c("符合工作要求","高於工作要求","低於工作要求")))
# data construction
glimpse(dta)
## Observations: 1,571
## Variables: 26
## $ SID <fctr> A1, A10, A100, A102, A103, A104, A105, A106, A107...
## $ Gender <fctr> 女, 女, 女, 女, 男, 女, 女, 女, 女, 女, 女, 男, 男, 男, 女, 男, 男...
## $ Sector <fctr> 國立（公立）, 國立（公立）, 私立, 國立（公立）, 國立（公立）, 國立（公立）, 國立（公立...
## $ EduLv <fctr> 碩士, 碩士, 普通大學, 普通大學, 高職, 普通大學, 普通大學, 普通大學, 普通大學, 普...
## $ SubEduOver <fctr> 符合工作要求, 高於工作要求, 符合工作要求, 符合工作要求, 符合工作要求, 符合工作要求, 符...
## $ Require <fctr> 碩士, 高中/高職, 普通大學, 普通大學, 高中/高職, 普通大學, 普通大學, 普通大學, 普...
## $ Field <fctr> 教育學群, 資訊學群, 外語學群, 教育學群, 工程學群, 文史哲學群, 文史哲學群, 大眾傳播學...
## $ City <fctr> 臺南市, 高雄市, 苗栗縣, 新北市, 高雄市, 南投縣, 嘉義市, 臺北市, 臺北市, 南投縣,...
## $ Category <fctr> 受雇於公營機關, 受雇於公營機關, 受雇於公營機關, 受雇於公營機關, 受雇者於私營企業, 受雇於...
## $ Staff <fctr> 10-29人, 50-99人, 50-99人, 10-29人, 2-9人, 100-199人, 1...
## $ Hours <int> 48, 40, 70, 50, 57, 51, 64, 50, 50, 47, 50, 60, 45...
## $ J_year <dbl> 2, 8, 4, 1, 21, 1, 6, 0, 1, 1, 17, 7, 3, 23, 1, 2,...
## $ J_total <dbl> 2, 8, 4, 1, 30, 1, 6, 0, 2, 2, 28, 7, 30, 26, 1, 1...
## $ income <fctr> 2萬以下, 2-3萬以下, 3-4萬以下, 4-5萬以下, 3-4萬以下, 4-5萬以下, 2萬以...
## $ SubMismatch <int> 4, 2, 3, 4, 5, 4, 5, 4, 3, 3, 4, 4, 4, 4, 5, 5, 3,...
## $ JobSat <int> 6, 4, 3, 5, 5, 6, 7, 5, 3, 6, 3, 5, 4, 7, 3, 4, 4,...
## $ EduZone <dbl> 5, 5, 4, 4, 2, 4, 4, 4, 4, 4, 5, 5, 3, 5, 4, 4, 4,...
## $ Region <fctr> 雲嘉南, 高屏澎, 桃竹苗, 北北基, 高屏澎, 中彰投, 雲嘉南, 北北基, 北北基, 中彰投,...
## $ Salary <dbl> 20000, 25000, 35000, 45000, 35000, 45000, 20000, 3...
## $ Age <dbl> 26, 34, 30, 25, 62, 25, 21, 24, 25, 26, 57, 35, 54...
## $ JobZone <dbl> 4, 3, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4, 4, 5, 4, 4, 3,...
## $ JobCor <int> 3, 1, 2, 3, 1, 2, 1, 3, 1, 1, 1, 3, 1, 1, 3, 1, 3,...
## $ Core <fctr> 核心關聯, 無關聯, 部分關聯, 核心關聯, 無關聯, 部分關聯, 無關聯, 核心關聯, 無關聯,...
## $ ObjOver <fctr> over, over, adequate, adequate, under, adequate, ...
## $ X <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA...
## $ JobZone_D <dbl> 1, 2, 0, 0, -1, 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, ...
# check and pick out
lapply(dta[,c("Sector", "Field", "City", "Region","EduLv", "SubEduOver", "ObjOver")], levels)
## $Sector
## [1] "私立" "國外學校" "國立（公立）"
##
## $Field
## [1] "大眾傳播學群" "工程學群" "文史哲學群" "外語學群"
## [5] "生命科學學群" "生物資源學群" "地球與環境學群" "法政學群"
## [9] "社會與心理學群" "建築與設計學群" "財經學群" "教育學群"
## [13] "資訊學群" "遊憩與運動學群" "管理學群" "數理化學群"
## [17] "醫藥衛生學群" "藝術學群"
##
## $City
## [1] "宜蘭縣" "花蓮縣" "金門縣" "南投縣" "屏東縣" "苗栗縣" "桃園市"
## [8] "高雄市" "基隆市" "雲林縣" "新北市" "新竹市" "新竹縣" "嘉義市"
## [15] "嘉義縣" "彰化縣" "臺中市" "臺北市" "臺東縣" "臺南市" "澎湖縣"
##
## $Region
## [1] "宜花東離島" "北北基" "桃竹苗" "中彰投" "雲嘉南"
## [6] "高屏澎"
##
## $EduLv
## [1] "技術學院" "博士" "碩士" "普通大學" "科技大學" "五專"
## [7] "三專" "二專" "高中" "高職" "國中"
##
## $SubEduOver
## [1] "符合工作要求" "高於工作要求" "低於工作要求"
##
## $ObjOver
## [1] "adequate" "over" "under"
## [1] "SID" "Gender" "Sector" "EduLv" "SubEduOver"
## [6] "Require" "Field" "City" "Category" "Staff"
## [11] "Hours" "J_year" "J_total" "income" "SubMismatch"
## [16] "JobSat" "EduZone" "Region" "Salary" "Age"
## [21] "JobZone" "JobCor" "Core" "ObjOver" "X"
## [26] "JobZone_D"
p <- dplyr::select(dta, -City, -income, -JobSat, -X)%>%
filter(Age <= 65 & Age >= 20 & Category != "自營者" )
p <- p[p$EduLv %in% c("博士","碩士", "普通大學", "科技大學", "技術學院"),]
p <- p %>% mutate( EduLv = factor(EduLv, levels=c("博士","碩士","普通大學","科技大學",
"技術學院")),
EduLv = relevel(EduLv, ref = "技術學院"),
Category = factor(Category,levels=c("受雇者於私營企業","受雇於公營機關")))
#
# Fit the full model
full<- lm(Salary ~ Gender+ Sector+ Field+ Region+ EduLv+ SubEduOver +SubMismatch +ObjOver +Core + J_year+ EduZone+ Hours, data = p)
# Fit the null model
null <- lm(Salary ~1, data= p)
#
# 最終模型設定
# 薪資對性別
a1 <- update(null,~.+Gender)
# 薪資對性別 / 公私立.教育程度.學群
a2 <- update(null,~.+Gender+Sector+ EduLv+ Field)
# 薪資對性別 / 公私立.教育程度.學群 / 地區.年資.工時
a3 <- update(null,~.+Gender+Sector+ EduLv+ Field+ Region+ J_year+ Hours)
# 薪資對性別 / 公私立.教育程度.學群 / 地區.年資.工時 / 自評過量
a41 <- update(null,~.+Gender+Sector+ EduLv+ Field+ Region+ J_year+ Hours+ SubEduOver)
# 薪資對性別 / 公私立.教育程度.學群 / 地區.年資.工時 / 自評過量.自評關聯
a42 <- update(null,~.+Gender+Sector+ EduLv+ Field+ Region+ J_year+ Hours+ SubEduOver+SubMismatch)
# 薪資對性別 / 公私立.教育程度.學群 / 地區.年資.工時 / 客評過量
a51 <- update(null,~.+Gender+Sector+ EduLv+ Field+ Region+ J_year+ Hours+ ObjOver)
# 薪資對性別 / 公私立.教育程度.學群 / 地區.年資.工時 / 客評過量.客評關聯
a52 <- update(null,~.+Gender+Sector+ EduLv+ Field+ Region+ J_year+ Hours+ ObjOver+ Core)
# 或許拿掉公司人數?
# 薪資對性別 / 公私立.教育程度.學群 / 地區.年資.工時 / 自評過量.自評關聯
a53 <- update(null,~.+Gender+Sector+ EduLv+ Field+ Region+ J_year+ Hours+ ObjOver+ Core)
model1
# 殘差圖
par(mfrow=c(2,2))
plot(a2)
## GVIF Df GVIF^(1/(2*Df))
## Gender 1.33 1 1.15
## Sector 1.46 2 1.10
## EduLv 1.50 4 1.05
## Field 1.93 17 1.02
#誤差獨立性
# P值不顯著，說明無相關
durbinWatsonTest(a2)
## lag Autocorrelation D-W Statistic p-value
## 1 0.041 1.91 0.108
## Alternative hypothesis: rho != 0
model2
# 殘差圖
par(mfrow=c(2,2))
plot(a3)
## GVIF Df GVIF^(1/(2*Df))
## Gender 1.35 1 1.16
## Sector 1.47 2 1.10
## EduLv 1.62 4 1.06
## Field 2.24 17 1.02
## Region 1.22 5 1.02
## J_year 1.08 1 1.04
## Hours 1.03 1 1.02
#誤差獨立性
# P值不顯著，說明無相關
durbinWatsonTest(a3)
## lag Autocorrelation D-W Statistic p-value
## 1 0.0326 1.93 0.228
## Alternative hypothesis: rho != 0
model3
# 殘差圖
par(mfrow=c(2,2))
plot(a42)
## GVIF Df GVIF^(1/(2*Df))
## Gender 1.36 1 1.17
## Sector 1.50 2 1.11
## EduLv 1.73 4 1.07
## Field 2.48 17 1.03
## Region 1.23 5 1.02
## J_year 1.09 1 1.05
## Hours 1.04 1 1.02
## SubEduOver 1.28 2 1.06
## SubMismatch 1.25 1 1.12
#誤差獨立性
# P值不顯著，說明無相關
durbinWatsonTest(a42)
## lag Autocorrelation D-W Statistic p-value
## 1 0.0514 1.9 0.05
## Alternative hypothesis: rho != 0
model4
# 殘差圖
par(mfrow=c(2,2))
plot(a52)
## GVIF Df GVIF^(1/(2*Df))
## Gender 1.35 1 1.16
## Sector 1.49 2 1.11
## EduLv 2.01 4 1.09
## Field 3.04 17 1.03
## Region 1.23 5 1.02
## J_year 1.09 1 1.04
## Hours 1.03 1 1.02
## ObjOver 1.60 2 1.13
## Core 1.32 2 1.07
#誤差獨立性
# P值不顯著，說明無相關
durbinWatsonTest(a52)
## lag Autocorrelation D-W Statistic p-value
## 1 0.0388 1.92 0.118
## Alternative hypothesis: rho != 0