Pembentukan model (Perhitungan matriks manual)

Untuk pembentukan model perhitungan matriks, saya menggunakan data bangkitan yang sengaja saya buat untuk memahami model regresi tiga peubah.

set.seed (1090)
n <- 100
p <- 4

X1 <- runif(n,5,15)
X2 <- rnorm(n,50,10)
X4 <- rpois(n, lambda=3)

x0 <- rep(1,n)

beta_0 <- 5
beta_1 <- 2.5
beta_2 <- -1.2
beta_4 <- 0.8

epsilon <- rnorm(n,0,5)

Y <- beta_0*x0 + beta_1*X1 + beta_2*X2 + beta_4*X4 + epsilon

data_simulasi <- data.frame(x0,X1,X2,X4,Y)

head(data_simulasi)
##   x0        X1       X2 X4         Y
## 1  1 13.153392 44.12235  8 -12.23923
## 2  1 12.013208 72.58093  3 -46.65112
## 3  1 11.509809 42.22955  0 -13.25668
## 4  1  8.161317 52.04329  2 -48.96924
## 5  1 11.975161 53.41138  4 -23.23427
## 6  1  7.715053 48.78234  5 -26.27616
str(data_simulasi)
## 'data.frame':    100 obs. of  5 variables:
##  $ x0: num  1 1 1 1 1 1 1 1 1 1 ...
##  $ X1: num  13.15 12.01 11.51 8.16 11.98 ...
##  $ X2: num  44.1 72.6 42.2 52 53.4 ...
##  $ X4: int  8 3 0 2 4 5 2 3 6 0 ...
##  $ Y : num  -12.2 -46.7 -13.3 -49 -23.2 ...
summary(data_simulasi)
##        x0          X1               X2              X4              Y          
##  Min.   :1   Min.   : 5.080   Min.   :25.53   Min.   : 0.00   Min.   :-66.031  
##  1st Qu.:1   1st Qu.: 7.791   1st Qu.:43.79   1st Qu.: 2.00   1st Qu.:-38.971  
##  Median :1   Median : 9.678   Median :49.86   Median : 4.00   Median :-24.875  
##  Mean   :1   Mean   :10.201   Mean   :50.59   Mean   : 3.57   Mean   :-26.753  
##  3rd Qu.:1   3rd Qu.:12.928   3rd Qu.:56.55   3rd Qu.: 5.00   3rd Qu.:-15.283  
##  Max.   :1   Max.   :14.969   Max.   :76.56   Max.   :12.00   Max.   :  9.465
par(mfrow=c(1,3))
plot(data_simulasi$X1, data_simulasi$Y, main="Scatter Plot X1 Terhadap Y",
     xlab="X1", ylab="Y", pch=19, col="blue")
plot(data_simulasi$X2, data_simulasi$Y, main="Scatter Plot X2 vs Y",
     xlab="X2", ylab="Y", pch=19, col="red")
plot(data_simulasi$X4, data_simulasi$Y, main="Scatter Plot X4 vs Y",
     xlab="X4", ylab="Y", pch=19, col="green")

cor(data_simulasi)
## Warning in cor(data_simulasi): the standard deviation is zero
##    x0         X1          X2          X4          Y
## x0  1         NA          NA          NA         NA
## X1 NA 1.00000000  0.05901414  0.05878883  0.4305404
## X2 NA 0.05901414  1.00000000 -0.02871238 -0.8002630
## X4 NA 0.05878883 -0.02871238  1.00000000  0.1374310
## Y  NA 0.43054039 -0.80026302  0.13743098  1.0000000

Dari observasi terhadap scatter plot dapat terlihat hubungan X1 dan Y mungkin lemah atau tidak signifikan dalam model regresi, X2 memiliki hubungan dengan Y yang berarti jika X2 naik maka Y cenderung turun, dan kemungkinan besar X4 tidak memiliki hubungan signifikan dengan Y.

X <- as.matrix(cbind(1,X1,X2,X4))
Y <- as.matrix(Y)
b <- solve(t(X)%*%X) %*% t(X) %*% Y

round(b,4)
##       [,1]
##     7.8389
## X1  2.5651
## X2 -1.2438
## X4  0.6064
b0 <- b[1]
b1 <- b[2]
b2 <- b[3]
b3 <- b[4]

print(paste("b0:", b0))
## [1] "b0: 7.8388915271183"
print(paste("b1:", b1))
## [1] "b1: 2.56513577544418"
print(paste("b2:", b2))
## [1] "b2: -1.24384647407756"
print(paste("b3:", b3))
## [1] "b3: 0.606356555637488"
sigma_kuadrat <- (t(Y) %*% Y - t(b) %*% t(X) %*% Y) / (n - p)

sigma_kuadrat <- max(sigma_kuadrat, 0)

Res_se <- sqrt(sigma_kuadrat)
print(round(Res_se, 3))
## [1] 5.418
df <- n - p
print(df)
## [1] 96
if (exists("x0") & exists("X1") & exists("X2") & exists("X4")) {
  Y_duga <- x0 + X1 * b1 + X2 * b2 + X4 * b3
} else {
  stop("Salah satu variabel X0, X1, X2, atau X4 tidak ditemukan dalam dataset")
}

Y_df <- data.frame(Y, Y_duga)
print(head(Y_df))
##           Y    Y_duga
## 1 -12.23923 -15.29034
## 2 -46.65112 -56.64495
## 3 -13.25668 -22.00285
## 4 -48.96924 -41.58627
## 5 -23.23427 -32.29222
## 6 -26.27616 -36.85580
if (!any(is.na(Y_df))) {
  R_squared <- (cor(Y_df$Y, Y_df$Y_duga))^2
  print(round(R_squared, 4))

  R_squared_adj <- 1 - ((1 - R_squared) * (n - 1) / (n - p))
  print(round(R_squared_adj, 4))
} else {
  print("Terdapat NA dalam data, cek kembali perhitungan!")
}
## [1] 0.8768
## [1] 0.873

Interpretasi Hasil Residual Standard Error (RSE) = 5.418 → Rata-rata penyimpangan antara nilai aktual Y dan nilai prediksi Ŷ (Y_duga) adalah 5.418 unit. Semakin kecil nilai ini, semakin baik model dalam memprediksi data.

Degree of Freedom (DF) = 96 → Derajat bebas model, yaitu jumlah data dikurangi jumlah parameter yang diestimasi.

R-Squared = 0.8768 (87.68%) → Model regresi mampu menjelaskan 87.68% variasi dalam Y berdasarkan variabel independen X0, X1, X2, X4.

Adjusted R-Squared = 0.873 (87.3%) → Setelah memperhitungkan penalti jumlah variabel dalam model, model tetap menjelaskan 87.3% variasi dalam Y, menunjukkan model cukup baik dan tidak overfitting.

galat <- Y - (b0 + b1 * X1 + b2 * X2 + b3 * X4)

KTReg <- sum((Y_duga - mean(Y))^2) / (p - 1)
print(paste("KTR:", KTReg))
## [1] "KTR: 8245.49147835991"
KTG <- sum(galat^2) / (n - p)
print(paste("KTG", KTG))
## [1] "KTG 29.3558072281818"
Fhit <- KTReg / KTG
print(paste("Fhit:", round(Fhit, 0)))
## [1] "Fhit: 281"
dbreg <- p - 1
print(paste("dbr:", dbreg))
## [1] "dbr: 3"
dbg <- n - p
print(paste("dbg:", dbg))
## [1] "dbg: 96"
p_value_F <- pf(Fhit, dbreg, dbg, lower.tail = F)
print(paste("p-value F:", p_value_F))
## [1] "p-value F: 2.20045169307955e-47"
se_b <- sqrt(sigma_kuadrat[1] * solve(t(X) %*% X))
## Warning in sqrt(sigma_kuadrat[1] * solve(t(X) %*% X)): NaNs produced
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: 3.4371"
print(paste("se_b1:", round(se_b1, 4)))
## [1] "se_b1: 0.1945"
print(paste("se_b2:", round(se_b2, 4)))
## [1] "se_b2: 0.0541"
print(paste("se_b3:", round(se_b3, 4)))
## [1] "se_b3: 0.2536"
t_b0 <- b0 / se_b0
t_b1 <- b1 / se_b1
t_b2 <- b2 / se_b2
t_b3 <- b3 / se_b3

print(paste("t_b0:", round(t_b0, 2)))
## [1] "t_b0: 2.28"
print(paste("t_b1:", round(t_b1, 2)))
## [1] "t_b1: 13.19"
print(paste("t_b2:", round(t_b2, 2)))
## [1] "t_b2: -23"
print(paste("t_b3:", round(t_b3, 2)))  
## [1] "t_b3: 2.39"
p_b0 <- 2 * pt(-abs(t_b0), df = n - p)
p_b1 <- 2 * pt(-abs(t_b1), df = n - p)
p_b2 <- 2 * pt(-abs(t_b2), df = n - p)
p_b3 <- 2 * pt(-abs(t_b3), df = n - p)

print(paste("p_b0:", p_b0))
## [1] "p_b0: 0.0247801612933666"
print(paste("p_b1:", p_b1))
## [1] "p_b1: 2.78890702883156e-23"
print(paste("p_b2:", p_b2))
## [1] "p_b2: 7.77718097038952e-41"
print(paste("p_b3:", p_b3)) 
## [1] "p_b3: 0.0187534653616471"
t_value <- qt(0.975, df = n - p)

BB_b0 <- b0 - t_value * se_b0
BA_b0 <- b0 + t_value * se_b0

BB_b1 <- b1 - t_value * se_b1
BA_b1 <- b1 + t_value * se_b1

BB_b2 <- b2 - t_value * se_b2
BA_b2 <- b2 + t_value * se_b2

BB_b3 <- b3 - t_value * se_b3
BA_b3 <- b3 + t_value * se_b3  

Batas.Bawah <- c(BB_b0, BB_b1, BB_b2, BB_b3)
Batas.Atas <- c(BA_b0, BA_b1, BA_b2, BA_b3)
Fit <- c(b0, b1, b2, b3)

Selang.Kepercayaan <- data.frame(
  "Batas bawah Selang (2.5%)" = round(Batas.Bawah, 6),
  "Fit" = round(Fit, 6),
  "Batas atas Selang (97.5%)" = round(Batas.Atas, 6)
)

rownames(Selang.Kepercayaan) <- c("Intercept", "b1", "b2", "b3")
print(Selang.Kepercayaan)
##           Batas.bawah.Selang..2.5..       Fit Batas.atas.Selang..97.5..
## Intercept                  1.016264  7.838892                 14.661519
## b1                         2.179142  2.565136                  2.951129
## b2                        -1.351192 -1.243846                 -1.136501
## b3                         0.102963  0.606357                  1.109751

Berdasarkan output yang diberikan, bentuk model regresi liniernya adalah:

𝑌=7.8389 +2.5651 X1 − 1.2438 X2 + 0.6064 X4 di mana:

Intercept (7.8389) adalah nilai prediksi Y ketika semua variabel independen (X1, X2, dan X3) bernilai nol. Koefisien X1 (2.5651) menunjukkan bahwa setiap kenaikan 1 unit pada X1 akan meningkatkan Y sebesar 2.5651, dengan asumsi variabel lain tetap. Koefisien X2 (-1.2438) menunjukkan bahwa setiap kenaikan 1 unit pada X2 akan menurunkan Y sebesar 1.2438, dengan asumsi variabel lain tetap. Koefisien X4 (0.6064) menunjukkan bahwa setiap kenaikan 1 unit pada X4 akan meningkatkan Y sebesar 0.6064, dengan asumsi variabel lain tetap.

  1. Uji F (Signifikansi Model)

    KTR (KT Regresi) = 8245.49

    KTG (KT Galat) = 29.36

    F hitung = 281 p-value F = 2.2 × 10^−47 Karena p-value sangat kecil (jauh lebih kecil dari 0.05), kita bisa menyimpulkan bahwa model regresi yang digunakan secara keseluruhan signifikan.

  2. Standar Error Koefisien Regresi

    se_b0 = 3.4371

    se_b1 = 0.1945

    se_b2 = 0.054

    se_b3 = 0.2536

  3. Uji t (Signifikansi Parameter)

    t_b0 = 2.28 (p = 0.0248) → Signifikan

    t_b1 = 13.19 (p = 2.79 × 10^−23) → Signifikan

    t_b2 = -23 (p = 7.78 × 10^−41) → Signifikan

    t_b3 = 2.39 (p = 0.0188) → Signifikan Semua parameter signifikan karena p-value < 0.05.

Pembentukan model (Fungsi lm())

Pada pembentukan model fungsi lm saya menggunakan data yang diberikan oleh asisten dosen.

library(readr)
## Warning: package 'readr' was built under R version 4.4.2
url <- "https://raw.githubusercontent.com/Fhar1st/STA1231_Anreg_P2/refs/heads/main/Minggu%204/data%20liver.csv"
data_liver <- read_csv(url)
## Rows: 36 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (8): No, X1, X2, X3, X4, X5, X6, Y
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
data_liver <- data_liver[, c("X1", "X2", "X4", "Y")]

head(data_liver)
## # A tibble: 6 × 4
##      X1    X2    X4     Y
##   <dbl> <dbl> <dbl> <dbl>
## 1 16.4    8.9  6.02  159.
## 2 26.7   21.2 12.1   197.
## 3 12.5   16.6  8.88  145.
## 4  8.45  22.9  7.46  140.
## 5 10.2   14.2  2.06  130.
## 6 19.5   17.4  7.54  163.
str(data_liver)
## tibble [36 × 4] (S3: tbl_df/tbl/data.frame)
##  $ X1: num [1:36] 16.36 26.68 12.49 8.45 10.19 ...
##  $ X2: num [1:36] 8.9 21.2 16.6 22.9 14.2 ...
##  $ X4: num [1:36] 6.02 12.07 8.88 7.46 2.06 ...
##  $ Y : num [1:36] 159 197 145 140 130 ...
summary(data_liver)
##        X1               X2               X4               Y        
##  Min.   : 4.550   Min.   : 1.270   Min.   : 2.060   Min.   :120.9  
##  1st Qu.: 8.502   1st Qu.: 6.335   1st Qu.: 4.065   1st Qu.:143.9  
##  Median :12.265   Median : 9.220   Median : 6.095   Median :160.7  
##  Mean   :14.680   Mean   :11.549   Mean   : 6.811   Mean   :169.7  
##  3rd Qu.:19.810   3rd Qu.:15.207   3rd Qu.: 8.998   3rd Qu.:191.8  
##  Max.   :35.410   Max.   :36.360   Max.   :15.000   Max.   :247.4
model <- lm(Y ~ X1 + X2 + X4, data = data_liver)

summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + X4, data = data_liver)
## 
## 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    
## X4            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 . 
## X4         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
## X4           0.2616205   7.643053
data_liver$Y_pred <- predict(model)

data_liver$residuals <- residuals(model)
  1. Persamaan Regresi Dari tabel Coefficients, model regresi yang diperoleh adalah: 𝑌=122.1006 + 0.9242 X1 + 0.6180 X2 + 3.9523 X4 Ini berarti:
  1. Signifikansi Koefisien (T-test dan P-value)
  1. Kualitas Model (R-squared & Adjusted R-squared)
  1. Uji ANOVA (F-Test)

5.Uji Pengaruh Masing-Masing Variabel (ANOVA Table)

Kesimpulan