#Tải một số package

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, union
library(ggplot2)

##PHẦN 1: TẢI DỮ LIỆU VÀ PHÂN TÍCH TỔNG QUÁT BỘ DỮ LIỆU

#Lấy dữ liệu từ tập dữ liệu có sẵn carData”

library (carData)

#upload dữ liệu Salaries vào R

data("Salaries")

#Gán dữ liệu là dạng data frame

df <- Salaries

#Hiển thị 6 quan sát đầu tiên

head(Salaries)

#Hiển thị 6 quan sát cuối cùng

tail(Salaries)

#Hiển thị tên các biến

names(Salaries)
## [1] "rank"          "discipline"    "yrs.since.phd" "yrs.service"  
## [5] "sex"           "salary"

#Hiển thị cấu trúc dữ liệu

str(Salaries)
## 'data.frame':    397 obs. of  6 variables:
##  $ rank         : Factor w/ 3 levels "AsstProf","AssocProf",..: 3 3 1 3 3 2 3 3 3 3 ...
##  $ discipline   : Factor w/ 2 levels "A","B": 2 2 2 2 2 2 2 2 2 2 ...
##  $ yrs.since.phd: int  19 20 4 45 40 6 30 45 21 18 ...
##  $ yrs.service  : int  18 16 3 39 41 6 23 45 20 18 ...
##  $ sex          : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 1 ...
##  $ salary       : int  139750 173200 79750 115000 141500 97000 175000 147765 119250 129000 ...

#Hiển thị tóm tắt dữ liệu

summary(Salaries)
##         rank     discipline yrs.since.phd    yrs.service        sex     
##  AsstProf : 67   A:181      Min.   : 1.00   Min.   : 0.00   Female: 39  
##  AssocProf: 64   B:216      1st Qu.:12.00   1st Qu.: 7.00   Male  :358  
##  Prof     :266              Median :21.00   Median :16.00               
##                             Mean   :22.31   Mean   :17.61               
##                             3rd Qu.:32.00   3rd Qu.:27.00               
##                             Max.   :56.00   Max.   :60.00               
##      salary      
##  Min.   : 57800  
##  1st Qu.: 91000  
##  Median :107300  
##  Mean   :113706  
##  3rd Qu.:134185  
##  Max.   :231545

#PHẦN 2: PHÂN TÍCH VÀ VẼ ĐỒ THỊ CÁC BIẾN

#Gán nhãn “salary_prof” cho biến Salary

salary_prof <- Salaries$salary

#Một vài thống kê cho biến Salary

summary(salary_prof)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   57800   91000  107300  113706  134185  231545
mean(salary_prof)
## [1] 113706.5
sd(salary_prof)
## [1] 30289.04
median(salary_prof)
## [1] 107300

#Vẽ biểu đồ histogram cho biến Salary

hist(salary_prof)

hist(salary_prof, col = "blue", main = "Professors Salary")

#Vẽ biểu đồ mật độ cho biến Salary plot(density(salary_prof))

#Vẽ biểu đồ boxplot cho biến Salary

boxplot(salary_prof)

#Vẽ 2 biểu đồ theo dạng 1 hàng 2 cột

par(mfrow=c(1,2))
hist(salary_prof)
plot(density(salary_prof))

#Vẽ 3 biểu đồ theo dạng 1 hàng 3 cột

par(mfrow=c(1,3))
hist(salary_prof)
plot(density(salary_prof))
boxplot(salary_prof)

#Do quan sát trên biểu đồ thấy biến Salary không theo phân phối chuẩn (lệch phải) #Nên ta tạo biên mới là hàm dạng log của Salary, đặt tên là log_salary_prof

log_salary_prof <- log(salary_prof)

#Vẽ biểu đồ histogram cho biến Log_salary_prof #Biểu đồ có dạng gần với phân phối chuẩn

hist(log_salary_prof)

#Vẽ biểu đồ mật độ cho biến Log_salary_prof

plot(density(log_salary_prof))

#Mô tả hàm học vị (rank) trong dữ liệu

summary(Salaries$rank)
##  AsstProf AssocProf      Prof 
##        67        64       266

#Vẽ đồ thị dạng boxplot của lương giáo sư (salary) theo học vị (rank)

boxplot(Salaries$salary ~ Salaries$rank)

boxplot(Salaries$salary ~ Salaries$rank, col="blue", main="Lương theo học vị", ylab="Tổng lương (USD)", xlab="Học vị")

#Cách 2: Vẽ đồ thị dạng boxplot của lương (salary) theo học hàm (rank)

ggplot(Salaries) + geom_boxplot(aes(x=rank, y=salary))

#Vẽ đồ thị dạng plot của tiền lương (salary) theo học vị (rank) và lĩnh vực nghiên cứu (displine)

ggplot(Salaries) + geom_boxplot(aes(x=rank, y=salary))

ggplot(Salaries) + geom_boxplot(aes(x=discipline, y=salary))

salary_discipline_rank <-  Salaries %>% ggplot() +
  geom_jitter(aes(rank,salary, colour=discipline)) + 
  geom_smooth(aes(rank, salary, colour=discipline), method=lm, se=FALSE) +
  labs(x = "Rank", y = "Tổng lương (USD",
       title = "Tiền lương với học hàm và lĩnh vực nghiên cứu")
salary_discipline_rank
## `geom_smooth()` using formula 'y ~ x'

#Mô tả số năm tính từ khi có học vị PhD trong tập dữ liệu

summary(Salaries$yrs.since.phd)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   12.00   21.00   22.31   32.00   56.00

#Vẽ đồ thị dạng plot giữa lương (salary) và số năm có PhD (yrs.since.phd)

plot(Salaries$salary ~ Salaries$yrs.since.phd)

plot(Salaries$salary ~ Salaries$yrs.since.phd, col="blue", main = "Lương theo năm có PhD", ylab = "Tổng lương (USD)", xlab = "Số năm có PhD")

#Mô tả lĩnh vực nghiên cứu (discipline) trong tập dữ liệu

summary(Salaries$discipline)
##   A   B 
## 181 216

#Vẽ đồ thị dạng boxplot giữa lương (salary) và lĩnh vực (discipline)

boxplot(Salaries$salary ~ Salaries$discipline)

boxplot(Salaries$salary ~ Salaries$discipline, col="blue", main="Lương theo lĩnh vực", ylab="Tổng lương (USD)", xlab="Lĩnh vực")

#Mô tả về giới tính giáo sư (sex) trong bảng dữ liệu

summary(Salaries$sex)
## Female   Male 
##     39    358

#Vẽ đồ thị dạng boxplot giữa lương (salary) và giới tính (sex)

boxplot(Salaries$salary ~ Salaries$sex)

boxplot(Salaries$salary ~ Salaries$sex, col="blue", main="Lương theo giới tính", ylab="Tổng lương (USD)", xlab="Giới tính")

#Mô tả về số năm làm việc (yrs.service) trong bảng dữ liệu

summary(Salaries$yrs.service)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    7.00   16.00   17.61   27.00   60.00

#Vẽ đồ thị dạng plot giữa lương (salary) và số năm làm việc (yrs.service)

plot(Salaries$salary ~ Salaries$yrs.service)

plot(Salaries$salary ~ Salaries$yrs.service, col="blue", main = "Lương theo năm làm việc", ylab = "Tổng lương (USD)", xlab = "Số năm làm việc")

#Vẽ 4 biểu đồ mức lương (salary) theo 4 biến (rank, yrs.since.phd, displine, sex) thành dạng 2 hàng 2 cột

par(mfrow=c(2,2))
boxplot(Salaries$salary ~ Salaries$rank, col="blue", main="Lương theo học vị", ylab="Tổng lương (USD)", xlab="Học vị")
plot(Salaries$salary ~ Salaries$yrs.since.phd, col="blue", main = "Lương theo năm có PhD", ylab = "Tổng lương (USD)", xlab = "Số năm có PhD")
boxplot(Salaries$salary ~ Salaries$discipline, col="blue", main="Lương theo lĩnh vực", ylab="Tổng lương (USD)", xlab="Lĩnh vực")
boxplot(Salaries$salary ~ Salaries$sex, col="blue", main="Lương theo giới tính", ylab="Tổng lương (USD)", xlab="Giới tính")

##PHẦN 3: PHÂN TÍCH HỒI QUY TUYẾN TÍNH

#Tính hệ số tương quan giữa mức lương (salary) và số năm có PhD (yrs.since.phd)

cor(Salaries$salary, Salaries$yrs.since.phd, use="everything", method=c("pearson"))
## [1] 0.4192311

#Vẽ biểu đồ ma trận tương quan cặp giữa các biến

pairs(Salaries)

##PHẦN 3.1: MÔ HÌNH HỒI QUY TUYẾN TÍNH ĐƠN

#Mô hình hồi quy tuyến tính đơn giữa lương (salary) và số năm có PhD (yrs.since.phd)

salary_phd <- lm(Salaries$salary ~ Salaries$yrs.since.phd)
summary(salary_phd)
## 
## Call:
## lm(formula = Salaries$salary ~ Salaries$yrs.since.phd)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -84171 -19432  -2858  16086 102383 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             91718.7     2765.8  33.162   <2e-16 ***
## Salaries$yrs.since.phd    985.3      107.4   9.177   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 27530 on 395 degrees of freedom
## Multiple R-squared:  0.1758, Adjusted R-squared:  0.1737 
## F-statistic: 84.23 on 1 and 395 DF,  p-value: < 2.2e-16

#Kết quả cho thấy #Ước lượng của hệ số tương quan và hệ số tự do đều có ý nghĩa thống kê (p-value < 0.05) #Dữ liệu giải thích được 17.37% mô hình thực nghiệm (Adjusted R-squared = 0.1737) #Hàm thực nghiệm là: salary = 91718.7 + 985.3 * yrs.since.phd

#Vẽ đồ thị mối quan hệ tuyến tính giữa mức lương (salary) và số năm có PhD (yrs.since.phd)

plot(Salaries$yrs.since.phd, Salaries$salary, col = "blue", main = "Mức lương và Số năm có PhD", ylab = "Tổng lương (USD)", xlab = "Số năm có PhD")
abline(salary_phd, col=2, lwd=2)

#Chẩn đoán mô hình thực nghiệm thông qua 4 loại các biểu

par(mfrow=c(2,2))
plot(salary_phd)

##PHẦN 3.2: MÔ HÌNH HỒI QUY TUYẾN TÍNH ĐA BIẾN

#Mô hình hồi quy tuyến tính của lương (salary) theo 4 biến độc lập là: rank, displine, sex, yrs.since.phd

full_model <- lm(salary~rank+discipline+yrs.since.phd+yrs.service+sex, data = df)

#Mô hình trên có thể viết tắt dưới dạng

full_model <- lm(salary~., data = df)
#Chạy kết quả phân tích hồi quy đa biến
summary(full_model)
## 
## Call:
## lm(formula = salary ~ ., data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65248 -13211  -1775  10384  99592 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    65955.2     4588.6  14.374  < 2e-16 ***
## rankAssocProf  12907.6     4145.3   3.114  0.00198 ** 
## rankProf       45066.0     4237.5  10.635  < 2e-16 ***
## disciplineB    14417.6     2342.9   6.154 1.88e-09 ***
## yrs.since.phd    535.1      241.0   2.220  0.02698 *  
## yrs.service     -489.5      211.9  -2.310  0.02143 *  
## sexMale         4783.5     3858.7   1.240  0.21584    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22540 on 390 degrees of freedom
## Multiple R-squared:  0.4547, Adjusted R-squared:  0.4463 
## F-statistic:  54.2 on 6 and 390 DF,  p-value: < 2.2e-16

#Đồ thị mô hình tổng thể của full_model

f1<- full_model$fitted.values
plot(f1,Salaries$salary,main="Mô hình hồi quy tuyến tính đa biến f1")
abline(lm(Salaries$salary~f1),col=2,lwd=2)

#Kết quả cho thấy #Câu hỏi 1: Tại sao lại chỉ xuất hiện một số biến là: rankAssocProf, rankProf, disciplineB, sexMale #Câu hỏi 2: Các biến khác tại sao ko có trong mô hình, ví dụ rankAsstProf, disciplineA, sexFemale #Câu hỏi 3: Có cần chuyển các biến phân loại (categorical) như biến rank (AsstProf, AssocProf, Prof) sang định dạng là số (0, 1, 2)

#Học hàm AssocProf, Prof tương quan dương với mức lương (salary) #Lĩnh vực B tương quan dương với mức lương (salary) #Thời gian nhận học vị PhD (yrs.since.phd) tương quan dương với mức lương (salary) #Số năm làm việc (yrs.service) tương quan âm với mức lương (salary). Đây là một hiện tượng bất thường. #Gới tính nam (sexMale) không có tương quan với mức lương (salary) #Số liệu giải thích 44.63% mô hình (Adjusted R-squared = 0.4463)

#Từ kết quả trên ta thực hiện 2 bước tiếp theo là #Bước 1: Tính hệ số tương quan giữa số năm nhận học vị PhD (yrs.since.phd) và thời gian làm việc (yrs.service) #Bước 2: Loại bỏ biến giới tính (sex) ra khỏi mô hình

#Hệ số tương quan là 0.91 là rất lớn nên có thể xảy ra hiện tượng đa cộng tuyến

cor(Salaries$yrs.since.phd, Salaries$yrs.service)
## [1] 0.9096491

#Do vậy ta phải loại bỏ 1 trong 2 biến trên khỏi mô hình. #Ta lựa chọn bỏ biến nào? Loại bỏ biến yrs.service vì có hệ số tương quan thấp hơn (r=0.334) biến yrs.since.phd (r=0.419)

cor(Salaries$salary, Salaries$yrs.since.phd)
## [1] 0.4192311
cor(Salaries$salary, Salaries$yrs.service)
## [1] 0.3347447

#Như vậy mô hình mới gồm các biến độc lập là: rank, displine và yrs.since.phd

new_model <- lm(salary~rank+discipline+yrs.since.phd, data = df)

#Chạy kết quả phân tích hồi quy đa biến cho mô hình mới

summary(new_model)
## 
## Call:
## lm(formula = salary ~ rank + discipline + yrs.since.phd, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -67395 -13480  -1536  10416  97166 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   71405.40    3278.32  21.781  < 2e-16 ***
## rankAssocProf 13030.16    4168.17   3.126   0.0019 ** 
## rankProf      46211.57    4238.52  10.903  < 2e-16 ***
## disciplineB   14028.68    2345.90   5.980 5.03e-09 ***
## yrs.since.phd    71.92     126.68   0.568   0.5706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22670 on 392 degrees of freedom
## Multiple R-squared:  0.4454, Adjusted R-squared:  0.4398 
## F-statistic: 78.72 on 4 and 392 DF,  p-value: < 2.2e-16

#Kết quả của mô hình mới #Học hàm AssocProf và Prof tương quan dương với mức lương (salary) #Lĩnh vực B tương quan dương với mức lương (salary) #Thời gian nhận học vị PhD (yrs.since.phd) không tương quan dương với mức lương (salary) #Số liệu giải thích 43.98% mô hình (Adjusted R-squared = 0.4398) #Ta tiếp tục loại biến yrs.since.phd ra khỏi mô hình

#Đồ thị mô hình tổng thể của new_model

f2<- new_model$fitted.values
plot(f2,Salaries$salary,main="Mô hình hồi quy tuyến tính đa biến f2")
abline(lm(Salaries$salary~f2),col=2,lwd=2)

#Chạy mô hình mới

new_model2 <- lm(salary~rank+discipline, data=df)
summary(new_model2)
## 
## Call:
## lm(formula = salary ~ rank + discipline, data = df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -65990 -14049  -1288  10760  97996 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      71944       3135  22.948  < 2e-16 ***
## rankAssocProf    13762       3961   3.475 0.000569 ***
## rankProf         47844       3112  15.376  < 2e-16 ***
## disciplineB      13761       2296   5.993 4.65e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22650 on 393 degrees of freedom
## Multiple R-squared:  0.445,  Adjusted R-squared:  0.4407 
## F-statistic:   105 on 3 and 393 DF,  p-value: < 2.2e-16

#Kết quả #Hàm thực nghiệm là salary = 71944 + (1-X1)47844 + X113762 + (1-X2)*13761 #Với X1=1 khi học vị là AssocProf, X1=0 khi học vị là Prof
#Với X2 = 1 khi lĩnh vực A, X2=0 khi lĩnh vực B

#Đồ thị mô hình tổng thể

f3<- new_model2$fitted.values
plot(f3,Salaries$salary,main="Mô hình hồi quy tuyến tính đa biến f3")
abline(lm(Salaries$salary~f3),col=2,lwd=2)

#Đồ thị so sánh 3 mô hình tổng thể f1, f2 và f3

par(mfrow=c(1,3))
f1<- full_model$fitted.values
plot(f1,Salaries$salary,main="Mô hình f1")
abline(lm(Salaries$salary~f1),col=2,lwd=2)
f2<- new_model$fitted.values
plot(f2,Salaries$salary,main="Mô hình f2")
abline(lm(Salaries$salary~f2),col=2,lwd=2)
f3<- new_model2$fitted.values
plot(f3,Salaries$salary,main="Mô hình f3")
abline(lm(Salaries$salary~f3),col=2,lwd=2)