Definisi

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.

Memasukkan Data

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

Eksplorasi Data

# 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.

Pembentukan Model (Fungsi Im())

Ringkasan Model Regresi

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

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

Selang Kepercayaan

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

Pembentukan model (Perhitungan matriks manual)

Parameter Regresi

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"

Koefisien determinasi dan penyesuaiannya

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 dan Galat Baku Parameter Regresi

# 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"

Signifikansi Parameter (nilai t)

# 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"

Selang Kepercayaan (1-alpha)*100%

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

Perbandingan Hasil Perhitungan Fungsi Im() vs Manual

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.

Perbandingan Dua Model

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.