library(dplyr)
library(tidyr)
library(readxl)
library(dichromat)
library(tidyverse)
library(skedastic)
library(car)
library(lmtest)
library(olsrr)Input Data
df <- readxl::read_excel("EPL 2020-2021.xlsx")
str(df)## tibble [293 × 12] (S3: tbl_df/tbl/data.frame)
## $ MarketValue : num [1:293] 1304 1130 956 695 782 ...
## $ Age : num [1:293] 21 24 29 28 26 21 21 27 22 19 ...
## $ Starts : num [1:293] 32 29 24 23 21 18 18 15 12 10 ...
## $ Mins : num [1:293] 2890 2602 2146 2010 1815 ...
## $ Goals : num [1:293] 6 6 0 7 0 4 4 2 6 2 ...
## $ Assists : num [1:293] 5 8 2 1 1 2 3 3 1 3 ...
## $ Passes_Attempted: num [1:293] 1881 826 1504 1739 1737 ...
## $ xG : num [1:293] 0.21 0.41 0.04 0.31 0.05 0.28 0.37 0.15 0.56 0.12 ...
## $ xA : num [1:293] 0.24 0.21 0.05 0.09 0.09 0.14 0.09 0.28 0.07 0.26 ...
## $ Yellow_Cards : num [1:293] 2 2 7 2 4 2 2 3 0 0 ...
## $ PlayerStatus : num [1:293] 0 1 1 1 1 1 1 1 0 0 ...
## $ TeamStatus : num [1:293] 1 1 1 1 1 1 1 1 1 1 ...Exploratory Data Analysis
ggplot(df) +
aes(x ="", y = MarketValue, fill = as.factor(TeamStatus)) +
geom_boxplot() +
scale_fill_hue(direction = 1) +
coord_flip() +
theme_minimal()
ggplot(df) +
aes(x ="", y = MarketValue, fill = as.factor(PlayerStatus)) +
geom_boxplot() +
scale_fill_hue(direction = 1) +
coord_flip() +
theme_minimal()
GGally::ggpairs(df[1:10])
Terlihat bahwa peubah yang memiliki korelasi paling tinggi dengan Y adalah Assists, lalu diikuti oleh Goals. Sebaran peubah respon terlihat menjulur ke kanan yang menunjukkan bahwa sebagian besar pemain liga inggris memiliki Market Value yang rendah.
Histogram
par(mfrow=c(3,3))
for(i in c(1:9)){
hist(df[[i]], col="green", main = paste("Histogram of",colnames(df)[i]), xlab=colnames(df)[i])
}
Barplot
par(mfrow=c(1,3))
green_palette <- colorRampPalette(c("lightgreen", "darkgreen"))
bar_colors <- green_palette(length(table(df$Yellow_Cards)))
barplot(table(df$Yellow_Cards), col=bar_colors, xlab="Category", ylab="Frequency",main = "Barplot of Yellow Cards")
barplot(table(df$TeamStatus), xlab="Category", ylab="Frequency", main = "Barplot of Team Status", col=c("green","darkgreen"))
barplot(table(df$PlayerStatus), col=c("green","darkgreen"), xlab="Category", ylab="Frequency",main = "Barplot of Player Status")
Boxplot
par(mfrow=c(1,3))
for(i in c(1:9)){
boxplot(df[i], main = paste("Boxplot of",colnames(df)[i]), col = "green")
}
Regresi Linear Berganda (OLS)
ols <- lm(MarketValue ~ ., data = df)
summary(ols)##
## Call:
## lm(formula = MarketValue ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -577.54 -131.35 -14.22 126.06 757.49
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 458.62942 79.28373 5.785 1.93e-08 ***
## Age -19.50909 3.17648 -6.142 2.78e-09 ***
## Starts -4.68537 12.07767 -0.388 0.6984
## Mins 0.01054 0.14408 0.073 0.9418
## Goals 28.74046 5.85908 4.905 1.58e-06 ***
## Assists 42.62719 8.62492 4.942 1.33e-06 ***
## Passes_Attempted 0.23938 0.05296 4.520 9.12e-06 ***
## xG 58.67790 105.32516 0.557 0.5779
## xA 32.05856 147.73562 0.217 0.8284
## Yellow_Cards 3.56849 8.12389 0.439 0.6608
## PlayerStatus 51.05874 28.39792 1.798 0.0733 .
## TeamStatus 266.16171 36.17773 7.357 2.07e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 218.6 on 281 degrees of freedom
## Multiple R-squared: 0.6613, Adjusted R-squared: 0.6481
## F-statistic: 49.88 on 11 and 281 DF, p-value: < 2.2e-16Cek multikolinearitas
vif(ols)## Age Starts Mins Goals
## 1.109478 116.682090 123.753891 3.582797
## Assists Passes_Attempted xG xA
## 2.613650 5.694110 1.980655 1.381770
## Yellow_Cards PlayerStatus TeamStatus
## 2.269907 1.105888 1.273737Starts dan Mins berkorelasi kuat, sehingga perlu eliminasi peubah.
df <- as.data.frame(df)
df <- df %>% dplyr::select(-Starts)
str(df)## 'data.frame': 293 obs. of 11 variables:
## $ MarketValue : num 1304 1130 956 695 782 ...
## $ Age : num 21 24 29 28 26 21 21 27 22 19 ...
## $ Mins : num 2890 2602 2146 2010 1815 ...
## $ Goals : num 6 6 0 7 0 4 4 2 6 2 ...
## $ Assists : num 5 8 2 1 1 2 3 3 1 3 ...
## $ Passes_Attempted: num 1881 826 1504 1739 1737 ...
## $ xG : num 0.21 0.41 0.04 0.31 0.05 0.28 0.37 0.15 0.56 0.12 ...
## $ xA : num 0.24 0.21 0.05 0.09 0.09 0.14 0.09 0.28 0.07 0.26 ...
## $ Yellow_Cards : num 2 2 7 2 4 2 2 3 0 0 ...
## $ PlayerStatus : num 0 1 1 1 1 1 1 1 0 0 ...
## $ TeamStatus : num 1 1 1 1 1 1 1 1 1 1 ...Cek ulang multikolinearitas
ols2 <- lm(MarketValue ~ ., data = df)
sumols <- summary(ols2)
vif(ols2)## Age Mins Goals Assists
## 1.095599 7.764766 3.558352 2.611615
## Passes_Attempted xG xA Yellow_Cards
## 5.694075 1.961644 1.362415 2.257229
## PlayerStatus TeamStatus
## 1.096598 1.273583Terlihat bahwa nilai VIF yang besar terdapat pada peubah Starting Eleven dan Minutes Played sehingga perlu dilakukan reduksi peubah pada peubah Starting Eleven.
sumols##
## Call:
## lm(formula = MarketValue ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -573.0 -133.4 -13.8 120.5 759.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 457.58820 79.11884 5.784 1.94e-08 ***
## Age -19.37127 3.15180 -6.146 2.70e-09 ***
## Mins -0.04357 0.03604 -1.209 0.2276
## Goals 28.55271 5.83026 4.897 1.64e-06 ***
## Assists 42.53381 8.60857 4.941 1.34e-06 ***
## Passes_Attempted 0.23933 0.05288 4.526 8.88e-06 ***
## xG 62.68093 104.66047 0.599 0.5497
## xA 38.84167 146.47613 0.265 0.7911
## Yellow_Cards 3.80402 8.08896 0.470 0.6385
## PlayerStatus 50.04902 28.23576 1.773 0.0774 .
## TeamStatus 266.00743 36.12102 7.364 1.97e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 218.3 on 282 degrees of freedom
## Multiple R-squared: 0.6611, Adjusted R-squared: 0.6491
## F-statistic: 55.02 on 10 and 282 DF, p-value: < 2.2e-16Output di atas menunjukkan bahwa tingkat keyakinan 95% terlihat bahwa intersep dan peubah Age, Goals, Assists, Passes Attempted, dan Klub Asal berpengaruh signifikan terhadap peubah respon (Market Value) dikarenakan p-value < alpha (0.05). Nilai residual SE sebesar 218,3 dengan nilai Adjusted R2 dan R2 yang cukup besar sehingga menunjukkan model semakin baik. Nilai Adjusted R2 diperoleh sebesar 0,649 berarti bahwa peubah penjelas dapat menjelaskan 64.9% keragaman pada peubah respon, sisanya dijelaskan oleh faktor lainnya.
Deteksi Pencilan
outt <- ols_plot_cooksd_bar(ols2)
# Observasi ke-
cat("Observasi ke-\n", which(!is.na(outt$data$txt)), "\nmerupakan pencilan dari model, yaitu sebanyak", sum(!is.na(outt$data$txt)), "Outlier")## Observasi ke-
## 16 18 19 22 24 26 27 34 41 43 44 47 48 60 61 76 89 90 102 117 118 127 128 149 189 202 220 234 237 261 293
## merupakan pencilan dari model, yaitu sebanyak 31 OutlierKebaikan model
sse_ols <- sum(ols2$residuals^2)
mse_ols <- mean(ols2$residuals^2)
r2_ols <- sumols$r.squared
aic_ols <- AIC(ols2)
bic_ols <- BIC(ols2)
data.frame(sse_ols,mse_ols,r2_ols,aic_ols,bic_ols)Diagnostik model
par(mfrow=c(2,2))
plot(ols)
# Pengujian Asumsi ####
# Normalitas
plot(ols2,2)
qqPlot(ols2$residuals, distribution = "norm", mean = mean(ols2$residuals),
sd = sd(ols2$residuals), main = "Uji Normalitas QQ-Plot Sisaan Model")
## [1] 90 76ks.test(ols2$residuals, "pnorm",
mean = mean(ols2$residuals), sd = sd(ols2$residuals))##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: ols2$residuals
## D = 0.056537, p-value = 0.3062
## alternative hypothesis: two-sided# Heteroskedastisitas
plot(ols2, 1)
glejser(ols2)# Kebebasan Galat
plot(ols2$residuals, type = 'o'); abline(h = 0, col = 'red')
randtests::runs.test(ols2$residuals, alternative ="two.sided")##
## Runs Test
##
## data: ols2$residuals
## statistic = -2.4621, runs = 126, n1 = 146, n2 = 146, n = 292, p-value =
## 0.01381
## alternative hypothesis: nonrandomnessdwtest(ols)##
## Durbin-Watson test
##
## data: ols
## DW = 1.7426, p-value = 0.01016
## alternative hypothesis: true autocorrelation is greater than 0#exact=p-value tepatSecara keseluruhan, dapat ditunjukan bahwa semua asumsi model regresi tidak terpenuhi. Oleh karena itu, perlu dilakukan penanganan asumsi (tidak dilakukan pada kasus ini). Cara lainnya adalah dengan menggunakan regresi pendekatan LAD yang tidak memerlukan asumsi.
Regresi Robust (LAD)
Penjelasan lengkap mengenai konsep LAD dapat dilihat di sini
library(lpSolve)
n = nrow(df) # Number of observation
k = ncol(df) # Number of Columns
# Define the objective function
objective <- c(rep(0,k), 1)
# Define the constraints
constraints1 <- cbind(rep(-1,n), as.matrix(-df[,2:k]), rep(-1,n))
constraints2 <- cbind(rep(-1,n), as.matrix(-df[,2:k]), rep(1,n))
constraints <- as.matrix(rbind(constraints1, constraints2))
directions <- c(rep("<=", n),
rep(">=", n))
rhs <- rep(-df$MarketValue, 2)
# Solve the linear programming problem
solution <- lp("min", objective.in = objective,
const.mat = constraints, const.dir = directions,
const.rhs = rhs)
# Print the solution
solution$solution## [1] 0.0000000 0.0000000 0.0000000 26.9549840 33.7183501 0.2934907
## [7] 0.0000000 0.0000000 0.0000000 0.0000000 64.4369443 718.7776987Berikut adalah rangkumannya:
Coefficients <- solution$solution[-12]
data.frame(Variable=c("Intercept",colnames(df[-1])),Coefficients)Output di atas menunjukkan bahwa intersep, peubah Age, Minutes Played, xG, xA, Yellow Card, dan Negara asal tidak berpengaruh signifikan terhadap model yang ditandai oleh nilai dugaan parameternya sama dengan nol. Sehingga, model terbaik yang didapat dari pendugaan parameter LAD adalah:
\[ \hat{y}_i=26.955 Goals_i+33.718 Assits_i+0.293 Passess Attempted_i+64.637 Klub Asal_i \]
Sehingga didapat nilai R-Squared model sebesar 0.5457 yang artinya model dapat menjelaskan sebesar 54.57% dari keragaman total peubah respons.
Kebaikan Model
X <- as.matrix(df %>% dplyr::select(-MarketValue))
X <- cbind(1,X)
yhat <- X %*% Coefficientsy <- df$MarketValue
# Calculate the residuals
residuals <- y - yhat
# Calculate the Sum of Squared Errors (SSE)
sse <- sum(residuals^2)
# Calculate the Mean Squared Error (MSE)
mse <- mean(residuals^2)
# Calculate the R-squared value
rsq <- 1 - sse / sum((y - mean(y))^2)
# Calculate the Akaike Information Criterion (AIC)
aic <- n * log(sse/n) + 2 * k
# Calculate the Bayesian Information Criterion (BIC)
bic <- n * log(sse/n) + k * log(n)
# Display the results
cat("Sum Squared Error:", sse, "\n")## Sum Squared Error: 18398102cat("R-squared:", rsq, "\n")## R-squared: 0.5360914cat("Akaike Information Criterion (AIC):", aic, "\n")## Akaike Information Criterion (AIC): 3258.943cat("Bayesian Information Criterion (BIC):", bic, "\n")## Bayesian Information Criterion (BIC): 3299.424Perbandingan
Pemodelan regresi linear berganda dengan metode OLS dan regresi linear robust dengan metode LAD dibandingkan dengan menggunakan beberapa metrik evaluasi yang sering digunakan. Metrik-metrik tersebut antara lain Mean Squared Error (MSE), Mean Absolute Percentage Error (MAPE), R-Squared, Akaike Information Creiterion (AIC), dan Bayesian Information Criterion (BIC). Semakin rendah nilai MSE, MAPE, AIC, dan BIC dari suatu model maka model tersebut dapat dikatakan lebih baik dibandingkan model yang lainnya. Berbanding terbalik dengan nilai R-squared, semakin tinggi nilai R-squared maka semakin baik model yang dihasilkan karena menunjukkan perbedaan yang semakin kecil antara data observasi dengan fitted values (Jim 2020).
mape <- function(actual, predicted) {
# Calculate the absolute percentage error for each observation
ape <- abs((actual - predicted) / actual)
# Calculate the mean absolute percentage error
mean(ape)
}
compare <- data.frame(MSE=c(mse_ols,mse),MAPE=c(mape(y,ols2$fitted.values),mape(y,yhat)),R_Squared=c(r2_ols,rsq),AIC=c(aic_ols,aic),BIC=c(bic_ols,bic))
compare <- t(compare)
compare <- round(compare,4)
colnames(compare) <- c("Ordinary Least Square", "Least Absolute Deviation")
as.data.frame(compare)Model regresi linear robust dengan metode LAD dan pendekatan linear programming memiliki kinerja yang lebih baik dalam beberapa metrik tertentu daripada model regresi linear berganda dengan metode OLS. Pertama, disebutkan bahwa model regresi linear robust dengan metode LAD memiliki nilai MAPE, AIC, dan BIC yang lebih rendah daripada metode OLS. Nilai-nilai yang lebih rendah pada metrik-metrik ini menunjukkan bahwa model LAD lebih akurat dalam memprediksi data yang memiliki pencilan. Metode LAD mampu menangani pencilan dengan lebih baik karena tidak terlalu dipengaruhi oleh observasi yang jauh dari pola umum data.
Namun, apabila dibandingkan dengan metrik MSE dan R-Squared, metode LAD tidak lebih baik daripada metode OLS. Metode OLS cenderung memberikan nilai MSE yang lebih rendah dan R-Squared yang lebih tinggi dibandingkan dengan metode LAD. MSE mengukur kesalahan rata-rata dari model prediksi terhadap data aktual, sedangkan R-Squared mengindikasikan seberapa baik model cocok dengan data. Dalam hal ini, metode OLS menghasilkan model yang lebih akurat secara keseluruhan dalam hal pengurangan kesalahan kuadrat rata-rata dan kemampuan model dalam menjelaskan keragaman data.
Berdasarkan perbandingan kedua model, dapat dikatakan bahwa metode LAD lebih unggul daripada metode OLS dalam memodelkan Market Value pemain Liga Premier Inggris musim 2020-2021 yang memiliki pencilan, seperti yang ditunjukkan oleh nilai MAPE, AIC, dan BIC yang lebih rendah. Namun, dalam hal pengurangan kesalahan kuadrat rata-rata dan kemampuan model dalam menjelaskan keragaman data, metode OLS tetap menjadi pilihan yang lebih baik, sebagaimana terindikasi oleh nilai MSE yang lebih rendah dan R-Squared yang lebih tinggi.
Kesimpulan dan Saran
Model regresi linear robust dengan metode LAD dan pendekatan linear programming merupakan model yang lebih baik dalam memodelkan data yang memiliki nilai pencilan, yaitu data Market Value pemain Liga Premier Inggris musim 2020-2021 karena memiliki nilai MAPE, AIC, dan BIC yang lebih rendah daripada model regresi linear berganda dengan metode OLS. Walaupun demikian, dalam konteks pengurangan kesalahan kuadrat rata-rata dan kemampuan model dalam menjelaskan keragaman data, metode OLS lebih baik daripada metode LAD karena memiliki nilai MSE yang lebih rendah dan R-Squared yang lebih tinggi. Faktor-faktor yang berpengaruh signifikan terhadap Market Value pemain Liga Premier Inggris musim 2020-2021 antara lain Goals, Assists, Passes Attempted, dan Klub Asal. Beberapa faktor ini dapat digunakan oleh para stakeholders suatu klub sepakbola dalam memperkirakan valuasi harga pasar seorang pemain sepak bola terutama di Liga Premier Inggris. Selain itu, dengan faktor-faktor tersebut para stakeholders juga dapat menentukan strategi yang lebih efektif dalam pencarian pemain dengan valuasi dan kemampuan yang sesuai.