Pertemuan 4
2026-02-24
<style>body { background-color: #fff0f5; font-family: “Poppins”, “Comic Sans MS”, sans-serif; color: #4a4a4a; line-height: 1.6; }
h1, h2, h3 { color: #d63384; font-weight: bold; text-align: left; }
code { background-color: #ffe6f0; padding: 3px 6px; border-radius: 6px; }
pre { background-color: #fff5fa; border-left: 3px solid #ff8fcf; padding: 10px; border-radius: 8px; }
strong { color: #c2185b; }Putri Dwi Suci Sri Mulyaningsih Sitanggang (M0401241033) # Membangkitkan data simulasi
set.seed(2006) #Pengetauran set data
n <- 100 #Jumlah objek amatan
p <- 3 #Jumlah parameter
x1 <- runif(n,3,8) #Bangkitan data distribusi seragam rentang dari 3 hingga 8
x2 <- rnorm(n,2,2) #Bangkitan data ditribusi normal dengan rerata 2, sd 2
x0 <- rep(1,n) #Vektor 1 matriks X
X <- data.frame(x0,x1,x2)
head(X)## x0 x1 x2
## 1 1 7.438921 1.814271
## 2 1 7.408826 2.214963
## 3 1 4.057758 4.614045
## 4 1 5.526996 4.568976
## 5 1 3.996362 2.808091
## 6 1 6.923440 1.632044Pembentukan error dan persamaan y (Komponen Acak)
e <- rnorm(n,0,2) #Bilangan acak normal rerata 0, sd 2
y <- 5 + 4*x1 + 3*x2 + e #Persamaan model y
rand.comp <- data.frame(e,y)
head(rand.comp)## e y
## 1 -0.0894841 40.10901
## 2 -1.0156221 40.26457
## 3 2.5617496 37.63492
## 4 -2.1073260 38.70759
## 5 -2.3220441 27.08768
## 6 -0.1646264 37.42526dt <- data.frame(y,x1,x2)
head(dt)## y x1 x2
## 1 40.10901 7.438921 1.814271
## 2 40.26457 7.408826 2.214963
## 3 37.63492 4.057758 4.614045
## 4 38.70759 5.526996 4.568976
## 5 27.08768 3.996362 2.808091
## 6 37.42526 6.923440 1.632044Eksplorasi data
par(mfrow=c(1,2))
plot(x1, y, pch=19)
plot(x2, y, pch=7)
Kedua plot eksplorasi memberikan gambaran pada kecenderungan korelasi antar peubah y dengan masing-masing peubah x berupa tren positif diikuti dengan keragaman yang cukup besar # Pembentukan model (Perhitungan matriks manual) ## Parameter Regresi
y <- as.matrix(y)
X <- as.matrix(cbind(x0,x1,x2))
b <- solve(t(X)%*%X)%*%t(X)%*%y;round(b,4)## [,1]
## x0 6.1623
## x1 3.7527
## x2 3.0819b0 <- b[1]
b1 <- b[2]
b2 <- b[3]
print(paste("b0:",b[1]))## [1] "b0: 6.16226224713655"print(paste("b1:",b[2]))## [1] "b1: 3.7526940698332"print(paste("b2:",b[3]))## [1] "b2: 3.08187031828733"Koefisien determinasi dan penyesuaiannya
sigma_kuadrat <- (t(y)%*%y-t(b)%*%t(X)%*%y)/(n-p)
Res_se <- sqrt(sigma_kuadrat)
round(Res_se,3)## [,1]
## [1,] 2.339nilai koefisien regresi menunjukkan bahwa, secara rata-rata, prediksi model menyimpang sekitar 2.339 unit dari nilai aktual \(y\)
#derajat bebas
df <- n-p
df## [1] 97y_duga <- b0+b1*x1+b2*x2
Y <- data.frame(y,y_duga);head(Y)## y y_duga
## 1 40.10901 39.66961
## 2 40.26457 40.79155
## 3 37.63492 35.60968
## 4 38.70759 40.98438
## 5 27.08768 29.81356
## 6 37.42526 37.17356R_squared <- (cor(y,y_duga))^2;round(R_squared,4)## [,1]
## [1,] 0.9364\(R^2\) mengukur proporsi variasi dalam variabel dependen yang dapat dijelaskan oleh variabel independen dalam model. Nilai \(R^2\) sebesar 0.9364 menunjukkan bahwa sekitar 93.64% dari variasi dalam \(y\) dapat dijelaskan oleh model regresi. Ini menunjukkan bahwa model memiliki kecocokan yang sangat baik dengan data.
R_squared_adj <- 1-((1-R_squared)*(n-1)/(n-p));round(R_squared_adj,4)## [,1]
## [1,] 0.9351Adjusted \(R^2\) memperhitungkan jumlah variabel independen dalam model dan memberikan penalti untuk model yang terlalu kompleks. Nilai adjusted \(R^2\) sebesar 0.9351 menunjukkan bahwa meskipun ada penalti untuk jumlah variabel, model masih mampu menjelaskan 93.51% dari variasi dalam \(y\) . # Uji F dan Std.Error parameter regresi
galat <- y-(b0+b1*x1+b2*x2)
KTReg <- sum((y_duga-mean(y))^2)/(p-1);print(paste("KTR:",KTReg))## [1] "KTR: 3908.26723296631"KTG <- sum(galat^2)/(n-p);print(paste("KTG",KTG))## [1] "KTG 5.46965945502427"Fhit <- KTReg/KTG;print(paste("Fhit:",round(Fhit,0)))## [1] "Fhit: 715"dbreg <- p-1;print(paste("dbr:",dbreg)) #Derajat bebas regresi## [1] "dbr: 2"dbg <- n-p;print(paste("dbg:",dbg)) #Derajat bebas galat## [1] "dbg: 97"pv<- pf(Fhit, dbreg, dbg, lower.tail = F);print(paste("pvalue:",pv))## [1] "pvalue: 9.01669301493964e-59"se_b <- sqrt(sigma_kuadrat[1]*solve(t(X)%*%X))## Warning in sqrt(sigma_kuadrat[1] * solve(t(X) %*% X)): NaNs producedse_b## x0 x1 x2
## x0 0.9269223 NaN NaN
## x1 NaN 0.1621392 NaN
## x2 NaN NaN 0.1091336se_b0 <- se_b[1,1];print(paste("se_b0:",round(se_b0,4)))## [1] "se_b0: 0.9269"se_b1 <- se_b[2,2];print(paste("se_b1:",round(se_b1,4)))## [1] "se_b1: 0.1621"se_b2 <- se_b[3,3];print(paste("se_b2:",round(se_b2,4)))## [1] "se_b2: 0.1091"#Signifikansi Parameter (nilai-t) ## t-value
t_b0 <- b0/se_b0;print(paste("t_b0:",round(t_b0,2)))## [1] "t_b0: 6.65"t_b1 <- b1/se_b1;print(paste("t_b1:",round(t_b1,2)))## [1] "t_b1: 23.14"t_b2 <- b2/se_b2;print(paste("t_b2:",round(t_b2,2)))## [1] "t_b2: 28.24"p-value
print(paste("p_b0:",2*pt(-abs(t_b0),df <- n-p)))## [1] "p_b0: 1.74733534165223e-09"print(paste("p_b1:",2*pt(-abs(t_b1),df <- n-p)))## [1] "p_b1: 2.77708386441887e-41"print(paste("p_b2:",2*pt(-abs(t_b2),df <- n-p)))## [1] "p_b2: 1.37860874619101e-48"Selang Kepercayaan (1-alfa)*100%
t <- qt(.975, df <- 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
Batas.Bawah <- as.matrix(c(round(BB_b0,6),round(BB_b1,6),round(BB_b2,6)))
Batas.Atas <- as.matrix(c(round(BA_b0,6),round(BA_b1,6),round(BA_b2,6)))
Fit <- (c(round(b0,6),round(b1,6),round(b2,6)))
Selang.Kepercayaan <- cbind(Batas.Bawah, Fit, Batas.Atas)
colnames(Selang.Kepercayaan ) <- c("Batas bawah Selang (2.5%)", "Fit", "Batas atas Selang (97.5%)")
rownames(Selang.Kepercayaan ) <- c("Intersep", "b1", "b2")
Selang.Kepercayaan## Batas bawah Selang (2.5%) Fit Batas atas Selang (97.5%)
## Intersep 4.322578 6.162262 8.001946
## b1 3.430893 3.752694 4.074496
## b2 2.865270 3.081870 3.298470Pembentukan model (Fungsi lm())
reg <- lm(y~x1+x2, data=dt)
summary(reg)##
## Call:
## lm(formula = y ~ x1 + x2, data = dt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1949 -1.5580 -0.1707 1.5380 6.4589
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1623 0.9269 6.648 1.75e-09 ***
## x1 3.7527 0.1621 23.145 < 2e-16 ***
## x2 3.0819 0.1091 28.239 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.339 on 97 degrees of freedom
## Multiple R-squared: 0.9364, Adjusted R-squared: 0.9351
## F-statistic: 714.5 on 2 and 97 DF, p-value: < 2.2e-16anova(reg)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## x1 1 3454.7 3454.7 631.61 < 2.2e-16 ***
## x2 1 4361.9 4361.9 797.46 < 2.2e-16 ***
## Residuals 97 530.6 5.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1confint(reg)## 2.5 % 97.5 %
## (Intercept) 4.322578 8.001946
## x1 3.430893 4.074496
## x2 2.865270 3.298470Perbandingan manual (matriks) dan fungsi lm()
koef <- as.matrix(reg$coefficients)
penduga <- cbind(b, koef)
colnames(penduga) <- c('matriks', 'fungsi lm')
rownames(penduga) <- c("intersep", "b1", "b2")
penduga## matriks fungsi lm
## intersep 6.162262 6.162262
## b1 3.752694 3.752694
## b2 3.081870 3.081870Berdasar perbandingan ini diperoleh hasil yang sama yang artinya penggunaan matriks secara manual sudah tepat dalam memodelkan regresi linear berganda bergitu pula sebaliknya.
Perbandingan 2 model
# Model 1
a <- summary(reg)
r.sq1 <- a$r.squared
adj_r.sq1 <- a$adj.r.squared
se1 <- a$sigma
# Model 2
reg2 <- lm(y~.-1, data=dt)
b <- summary(reg2)
r.sq2 <- b$r.squared
adj_r.sq2 <- b$adj.r.squared
se2 <- b$sigma
##===Membandingkan kebaikan model===
model1 <- as.matrix(c(r.sq1,adj_r.sq1,se1))
model2 <- as.matrix(c(r.sq2,adj_r.sq2,se2))
tabel <- cbind(model1, model2)
colnames(tabel) <- c("Model 1", "Model 2")
rownames(tabel) <- c("R-Square", "Adj R-Square", "Standar Error Sisaan")
tabel## Model 1 Model 2
## R-Square 0.9364381 0.9933879
## Adj R-Square 0.9351275 0.9932529
## Standar Error Sisaan 2.3387303 2.8072429