y ~ x1 + x2 + x4
Pada RUntuk 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.
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.
Standar Error Koefisien Regresi
se_b0 = 3.4371
se_b1 = 0.1945
se_b2 = 0.054
se_b3 = 0.2536
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.
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)
5.Uji Pengaruh Masing-Masing Variabel (ANOVA Table)
X1 (p = 0.001731) → Signifikan pada tingkat 1%, berarti X1 berkontribusi secara signifikan terhadap model.
X2 (p = 0.064461) → Tidak signifikan pada tingkat 5%, berarti X2 kurang berpengaruh dalam model ini.
X4 (p = 0.036619) → Signifikan pada tingkat 5%, berarti X4 memiliki pengaruh terhadap Y.
Kesimpulan
Model secara keseluruhan signifikan (F-test p-value < 0.05).
Variabel X4 dan X1 berpengaruh signifikan terhadap Y, sedangkan X2 tidak signifikan.
Model hanya mampu menjelaskan 38.6% variasi Y, sehingga mungkin masih ada variabel lain yang berpengaruh.