Metode statistik guna menganalisis hubungan antara satu variabel dependen dengan dua/lebih variabel dependen adalah Regresi Linear Berganda. model ini merupakan pengembangan dari regresi linear sederhana, yang memungkinkan untuk memahami bagaimana bebeapa faktor secara simultan memengaruhi variabel target. Dalam analisis kali ini, kita meggunakan dataset “Data Liver” yang terdiri dari 7 variabel. Namun, untuk analisis kali ini hanya akan digunakan 4 variabel saja dari Data Liver, yaitu X1, X2, X4, Y.
n <- 36 #Jumlah objek amatan
p <- 4 #Jumlah parameter
# Memasukkan Data Liver
data_liver <- read.csv("D:/Documents/Downloads/data liver (1).csv", header = TRUE)
# mengambil 4 variabel yang dibuthkan dan menambahkan satu variabel x0 sebagai data untuk intersep.
df <- data_liver[, c("Y", "X1", "X2", "X4"), drop = FALSE]
colnames(df) <- c("y", "x1", "x2", "x3")
df <- as.data.frame(df)
print(df)
## y x1 x2 x3
## 1 158.76 16.36 8.90 6.02
## 2 197.19 26.68 21.22 12.07
## 3 144.73 12.49 16.62 8.88
## 4 140.06 8.45 22.86 7.46
## 5 129.71 10.19 14.23 2.06
## 6 162.59 19.53 17.35 7.54
## 7 178.48 20.65 10.48 4.88
## 8 120.90 22.96 14.23 3.69
## 9 191.24 21.22 21.64 11.94
## 10 150.03 8.11 3.16 8.82
## 11 173.44 24.74 7.84 3.68
## 12 211.98 11.38 15.71 7.20
## 13 193.49 15.82 15.04 9.89
## 14 164.04 8.36 9.01 3.40
## 15 156.97 12.04 9.72 6.03
## 16 208.36 10.97 4.58 5.55
## 17 154.62 7.97 9.33 4.17
## 18 137.38 7.46 6.11 2.99
## 19 180.15 29.09 15.71 9.35
## 20 228.47 10.30 8.54 10.78
## 21 153.62 7.82 4.41 4.19
## 22 121.31 14.71 6.29 6.16
## 23 157.37 8.54 6.73 5.52
## 24 211.27 23.05 11.34 3.00
## 25 178.16 13.12 5.86 10.92
## 26 174.89 7.41 9.11 5.50
## 27 142.98 14.59 5.59 3.75
## 28 165.59 8.52 6.52 6.92
## 29 141.54 18.97 6.35 5.61
## 30 238.22 35.41 36.36 15.00
## 31 138.42 4.55 1.27 2.83
## 32 247.45 22.59 28.70 10.35
## 33 140.27 9.21 4.55 7.92
## 34 216.06 18.32 11.61 8.07
## 35 144.18 5.69 6.88 2.78
## 36 156.22 11.21 11.92 10.29
# Membuat scatterplot untuk melihat persebaran data
par(mfrow=c(1,3))
plot(df$x1, df$y, pch = 19, col = "red",
xlab = "X1", ylab = "Y", main = "Scatterplot X1 vs Y")
plot(df$x2, df$y, pch = 19, col = "Green",
xlab = "X2", ylab = "Y", main = "Scatterplot X2 vs Y")
plot(df$x3, df$y, pch = 19, col = "Blue",
xlab = "X3", ylab = "Y", main = "Scatterplot X3 vs Y")
Ketiga plot eksplorasi menggambarkan hubungan antara peubah y dengan tiap variabel x1, x2, x3. Tampak ada kecenderungan tren positif, terutama untuk x3 vs y, meskipun penyebaran data memperlihatkan keragaman yang cukup besar. Ini mengindikasikan bahwa meskipun x1, x2, x3 memiliki hubungan positif terhadap y, hubungan tersebut tidak begitu kuat atau mungkin ada faktor lain yang lebih memengaruhi y.
model <- lm(y ~ x1 + x2 + x3, data = df)
summary(model)
##
## Call:
## lm(formula = y ~ x1 + x2 + x3, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -45.799 -16.805 -1.352 9.696 51.355
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 122.1006 11.9983 10.177 1.48e-11 ***
## x1 0.9242 0.8127 1.137 0.2639
## x2 0.6180 0.9055 0.682 0.4999
## x3 3.9523 1.8119 2.181 0.0366 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.97 on 32 degrees of freedom
## Multiple R-squared: 0.386, Adjusted R-squared: 0.3284
## F-statistic: 6.706 on 3 and 32 DF, p-value: 0.001223
anova(model)
## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 8501.9 8501.9 11.6922 0.001731 **
## x2 1 2666.8 2666.8 3.6674 0.064461 .
## x3 1 3459.9 3459.9 4.7582 0.036619 *
## Residuals 32 23268.7 727.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(model)
## 2.5 % 97.5 %
## (Intercept) 97.6609117 146.540358
## x1 -0.7310951 2.579564
## x2 -1.2265068 2.462407
## x3 0.2616205 7.643053
Dalam regresi linear berganda, estimasi parameter dilakukan memakai metode Ordinary Least Square yang dinotasikan dalam bentuk matriks berikut: \[ b_{(k+1)\times1}=(X'X)_{(k+1)\times(k+1)}^{-1}X'_{(k+1)\times{n}}y_{(n\times1)} \]
# Menentukan estimasi parameter regressi dengan rumus OLS
x0 <- rep(1, nrow(df))
Y <- as.matrix(df$y)
X <- as.matrix(cbind(x0, df$x1, df$x2, df$x3))
b <- solve(t(X)%*%X)%*%t(X)%*%Y;round(b,4)
## [,1]
## x0 122.1006
## 0.9242
## 0.6180
## 3.9523
# Menampilkan nilai koefisien regresi
b0 <- b[1]
b1 <- b[2]
b2 <- b[3]
b3 <- b[4]
print(paste("b0:",b[1]))
## [1] "b0: 122.1006349907"
print(paste("b1:",b[2]))
## [1] "b1: 0.924234526233408"
print(paste("b2:",b[3]))
## [1] "b2: 0.617950334052231"
print(paste("b3:",b[4]))
## [1] "b3: 3.95233664011293"
Dalam RLB, varians residual diperlukan untuk memgnukur seberapa jauh oenyebaran data observasi dari model yang siestimasi. Varians residual bisa dihitung dengan rumus di bawah: \[ \hat{\sigma}^2=\frac{SSR_{esidual}}{n-p}=\frac{y'y-\hat{\beta}x'y}{n-p} \]
# Koefisien regresi dengan metode OLS
sigma_kuadrat <- (t(Y)%*%Y-t(b)%*%t(X)%*%Y)/(n-p)
Res_se <- sqrt(sigma_kuadrat)
round(Res_se,3)
## [,1]
## [1,] 26.966
Nilai galat baku residual menunjukkan bahwa secara rata-rata penyimpangan prediksi model dari nilai aktual y yaitusekitar 26.966
# Memprediksi nilai y_duga
y_duga <- b0 + b1 * df$x1 + b2 * df$x2 + b3 * df$x3
Y <- data.frame(y = df$y, y_duga = y_duga)
head(Y)
## y y_duga
## 1 158.76 166.5139
## 2 197.19 207.5768
## 3 144.73 179.0114
## 4 140.06 173.5212
## 5 129.71 148.4538
## 6 162.59 180.6730
R_squared <- (cor(df$y, y_duga))^2 # Gunakan df$y, bukan y
round(R_squared, 4)
## [1] 0.386
Koefisien determinasi (\(R^2\)) menunjukkan besaran proporsi variasi dalam variabel dependen (\(y\)) yang bisa dijelaskan oleh variabel independen dalam model regresi. Dengan nilai \(R^2\) sebesar 0.386, sekitar 38.6% dari variasi dalam \(y\) bisa dijelaskna oleh model ini. Hal ini mengisyaratkan bahwa model ini kurang baik dalam menjelaskan y secara keseluruhan,
R_squared_adj <- 1-((1-R_squared)*(n-1)/(n-p));
print(round(R_squared_adj,4))
## [1] 0.3284
Adjusted \(R^2\) menunjukkan jumlah variabel independen dalam model dan memberikan penalti jika model terlalu kompleks. Dengan nilai adjusted \(R^2\) sebesar 0.3824, ini berarti setelah memperhitungkan jumlah variabel, model masih mampu menjelaskan sekitar 32.84% dari variasi dalam y. Artinya model ini masih memiliki keterbatasan dalam menjelaskan variabel dependen secara menyeluruh.
# Uji F
dbr <- p-1
dbg <- n-p
galat <- df$y -(b0+b1*df$x1+b2*df$x2+b3*df$x3)
KTR <- sum((y_duga-mean(df$y))^2)/(p-1)
KTG <- sum(galat^2)/(n-p)
Fhit <- KTR/KTG
pval <- pf(Fhit, dbr, dbg, lower.tail = F)
print(paste("db regresi:",dbr)) #Derajat bebas regresi
## [1] "db regresi: 3"
print(paste("db galat:",dbg)) #Derajat bebas regresi
## [1] "db galat: 32"
print(paste("KT regresi:",KTR)) #Kuadrat tengah regresi
## [1] "KT regresi: 4876.19884277776"
print(paste("KT galat:",KTG)) #Kuadrat tengah galat
## [1] "KT galat: 727.147978975695"
print(paste("Fhit:",round(Fhit, 4))) #Nilai F hitung
## [1] "Fhit: 6.7059"
print(paste("P-value:",round(pval, 4))) #Nilai p-value
## [1] "P-value: 0.0012"
#p-value untuk uji F
pf(Fhit, dbr, dbg, lower.tail = F)
## [1] 0.001223205
#Standar error koefisien regresi
se_b <- sqrt(sigma_kuadrat[1]*solve(t(X)%*%X))
se_b
## x0
## x0 11.998293 NaN 1.1829188 NaN
## NaN 0.8126577 NaN NaN
## 1.182919 NaN 0.9055069 NaN
## NaN NaN NaN 1.811898
se_b0 <- se_b[1,1]
se_b1 <- se_b[2,2]
se_b2 <- se_b[3,3]
se_b3 <- se_b[4,4]
print(paste("se_b0:",round(se_b0,4)))
## [1] "se_b0: 11.9983"
print(paste("se_b1:",round(se_b1,4)))
## [1] "se_b1: 0.8127"
print(paste("se_b2:",round(se_b2,4)))
## [1] "se_b2: 0.9055"
print(paste("se_b3:",round(se_b3,4)))
## [1] "se_b3: 1.8119"
# t-value
t_b0 <- b0/se_b0
t_b1 <- b1/se_b1
t_b2 <- b2/se_b2
t_b3 <- b2/se_b3
print(paste("t_b0:",round(t_b0,2)))
## [1] "t_b0: 10.18"
print(paste("t_b1:",round(t_b1,2)))
## [1] "t_b1: 1.14"
print(paste("t_b2:",round(t_b2,2)))
## [1] "t_b2: 0.68"
print(paste("t_b3:",round(t_b3,2)))
## [1] "t_b3: 0.34"
# p-value
print(paste("p_b0:",2*pt(-abs(t_b0),db <- n-p)))
## [1] "p_b0: 1.47506344832897e-11"
print(paste("p_b1:",2*pt(-abs(t_b1),db <- n-p)))
## [1] "p_b1: 0.263858037514289"
print(paste("p_b2:",2*pt(-abs(t_b2),db <- n-p)))
## [1] "p_b2: 0.499874045483506"
print(paste("p_b3:",2*pt(-abs(t_b3),db <- n-p)))
## [1] "p_b3: 0.735294176287717"
t <- qt(.975, db <- n-p)
BB_b0 <- b0-t*se_b0
BA_b0 <- b0+t*se_b0
BB_b1 <- b1-t*se_b1
BA_b1 <- b1+t*se_b1
BB_b2 <- b2-t*se_b2
BA_b2 <- b2+t*se_b2
BB_b3 <- b3-t*se_b3
BA_b3 <- b3+t*se_b3
BB <- as.matrix(c(round(BB_b0,5),round(BB_b1,5),round(BB_b2,5),round(BB_b3,5)))
BA <- as.matrix(c(round(BA_b0,5),round(BA_b1,5),round(BA_b2,5),round(BA_b3,5)))
Fit <- (c(round(b0,5),round(b1,5),round(b2,5),round(b3,5)))
Selang.Kepercayaan <- cbind(BB, Fit, BA)
colnames(Selang.Kepercayaan ) <- c("Batas bawah Selang (2.5%)", "Fit", "Batas atas Selang (97.5%)")
rownames(Selang.Kepercayaan ) <- c("Intersep", "b1", "b2", "b3")
Selang.Kepercayaan
## Batas bawah Selang (2.5%) Fit Batas atas Selang (97.5%)
## Intersep 97.66091 122.10063 146.54036
## b1 -0.73110 0.92423 2.57956
## b2 -1.22651 0.61795 2.46241
## b3 0.26162 3.95234 7.64305
koef <- as.matrix(model$coefficients)
penduga <- cbind(b, koef)
colnames(penduga) <- c('matriks', 'fungsi lm')
rownames(penduga) <- c("intersep", "b1", "b2", "b3")
penduga
## matriks fungsi lm
## intersep 122.1006350 122.1006350
## b1 0.9242345 0.9242345
## b2 0.6179503 0.6179503
## b3 3.9523366 3.9523366
Berdasarkan perbandingan di atas, didapatkan hasil yang sama, ini berarti penggunaan matriks secara manual sudah tepat dalam memodelkan RLB, begitupun sebaliknya.
p <- summary(model)
r_sq1 <- p$r.sq
adj_r_sq1 <- p$adj.r.squared
res_se1 <- p$sigma
df <- as.data.frame(df)
model2 <- lm(y ~ . -1, data = df)
q <- summary(model2)
r_sq2 <- q$r.squared
adj_r_sq2 <- q$adj.r.squared
res_se2 <- q$sigma
# Membandingkan performa kedua model
model1 <- as.matrix(c(r_sq1, adj_r_sq1, res_se1))
model2 <- as.matrix(c(r_sq2, adj_r_sq2, res_se2))
Tabel <- cbind(model1, model2)
# Menetapkan nama kolom dan baris agar lebih informatif
colnames(Tabel) <- c("Model 1", "Model 2")
rownames(Tabel) <- c("R-Squared", "Adjusted R-Squared", "Residual Std. Error")
# Menampilkan hasil perbandingan kedua model
Tabel
## Model 1 Model 2
## R-Squared 0.386006 0.9082997
## Adjusted R-Squared 0.328444 0.8999633
## Residual Std. Error 26.965682 54.6540148
Perbandingan ini menunjukkan cara mengevaluasi dua model regresi berganda. Hasil analisis mengindikasikan bahwa Adjusted \(R^2\) lebih efektif dibandingkan \(R^2\) dalam menentukan model yang paling sesuai.