TUGAS 5 REGRESI SPLINE TRUNCATED
Muhammad Farhan Nur Ismail
2026-10-11
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(MASS)
library(car)## Loading required package: carDatalibrary(pastecs)## Warning: package 'pastecs' was built under R version 4.6.1library(pracma)## Warning: package 'pracma' was built under R version 4.6.1##
## Attaching package: 'pracma'## The following object is masked from 'package:car':
##
## logitlibrary(Matrix)##
## Attaching package: 'Matrix'## The following objects are masked from 'package:pracma':
##
## expm, lu, tril, triulibrary(ggplot2)
library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.6.1data=read.table(file.choose(), header=TRUE) # ganti dengan nama file saat knit
data## Y X1 X2 X3 X4
## 1 4.52 67.88 13.00 14.86 72.04
## 2 4.97 71.02 12.89 14.38 51.19
## 3 5.70 61.98 13.58 18.56 73.59
## 4 5.45 68.96 12.78 29.61 43.00
## 5 5.08 67.40 13.31 14.89 74.71
## 6 6.22 69.04 12.55 55.70 71.41
## 7 3.49 76.22 12.50 10.46 67.09
## 8 9.00 62.90 14.13 15.59 80.01
## 9 8.26 65.16 14.70 75.04 30.11
## 10 9.46 69.24 12.90 31.21 80.02
## 11 3.71 74.28 12.60 20.67 67.03
## 12 3.91 75.81 12.08 8.67 47.87
## 13 3.38 71.78 12.39 10.74 65.98
## 14 7.55 64.14 12.33 8.54 25.74
## 15 3.52 70.38 11.56 19.95 65.77
## 16 7.30 60.75 11.79 28.13 67.17
## 17 4.50 75.57 12.02 14.71 36.88
## 18 4.02 74.09 12.04 10.28 65.69
## 19 3.39 77.53 11.57 6.56 64.76
## 20 2.70 73.93 11.15 25.23 25.55
## 21 3.71 65.53 11.81 4.21 22.68
## 22 7.14 67.71 13.64 28.93 67.95
## 23 12.36 60.05 13.64 37.68 79.44
## 24 8.78 63.84 12.89 10.14 41.94
## 25 3.57 72.03 11.96 10.16 69.38
## 26 4.96 64.68 11.92 17.31 68.86
## 27 3.87 72.55 12.28 11.73 59.18
## 28 2.93 74.61 12.38 5.76 66.22
## 29 3.73 70.17 11.86 6.35 70.11
## 30 2.24 73.15 12.10 4.65 48.85
## 31 3.90 71.15 12.19 4.83 68.84
## 32 4.49 70.08 12.88 3.30 65.59
## 33 3.07 69.27 12.59 14.48 72.19
## 34 6.95 70.16 12.36 15.39 50.71
## 35 2.46 76.50 12.37 9.17 68.82
## 36 8.32 62.07 13.92 21.93 77.10
## 37 5.54 66.82 14.80 6.11 79.10
## 38 5.08 66.44 13.29 11.21 71.94
## 39 4.45 67.38 12.99 18.72 41.10
## 40 4.83 67.81 12.20 5.75 66.97
## 41 4.14 66.91 12.63 26.30 45.79
## 42 5.86 65.65 13.73 38.11 75.83
## 43 4.76 73.01 12.71 33.79 72.87
## 44 5.25 67.41 12.69 56.63 71.31
## 45 4.98 70.04 12.90 45.81 69.48
## 46 4.21 64.68 12.54 45.53 70.22
## 47 5.29 71.54 12.48 51.50 30.59
## 48 4.70 65.50 12.11 66.27 68.03
## 49 2.83 70.50 12.47 67.99 70.51
## 50 4.30 65.04 11.98 40.96 37.58
## 51 5.69 64.55 12.51 48.22 68.68
## 52 2.63 72.77 12.40 43.84 68.45
## 53 2.49 71.22 11.77 51.25 70.81
## 54 2.91 77.73 12.82 54.56 21.39
## 55 3.10 63.93 11.74 62.92 67.98
## 56 5.95 62.71 14.94 61.01 80.77stat.desc(data)## Y X1 X2 X3
## nbr.val 56.0000000000 5.600000000e+01 56.00000000000 56.0000000000
## nbr.null 0.0000000000 0.000000000e+00 0.00000000000 0.0000000000
## nbr.na 0.0000000000 0.000000000e+00 0.00000000000 0.0000000000
## min 2.2400000000 6.005000000e+01 11.15000000000 3.3000000000
## max 12.3600000000 7.773000000e+01 14.94000000000 75.0400000000
## range 10.1200000000 1.768000000e+01 3.79000000000 71.7400000000
## sum 277.6000000000 3.863250000e+03 708.36000000000 1476.2800000000
## median 4.5100000000 6.914000000e+01 12.50500000000 18.6400000000
## mean 4.9571428571 6.898660714e+01 12.64928571429 26.3621428571
## SE.mean 0.2715868131 6.008087232e-01 0.10790864858 2.6750203387
## CI.mean.0.95 0.5442721359 1.204047587e+00 0.21625376425 5.3608605550
## var 4.1305262338 2.021438282e+01 0.65207948052 400.7210935065
## std.dev 2.0323696105 4.496040794e+00 0.80751438409 20.0180192204
## coef.var 0.4099881059 6.517266149e-02 0.06383873385 0.7593471945
## X4
## nbr.val 56.0000000000
## nbr.null 0.0000000000
## nbr.na 0.0000000000
## min 21.3900000000
## max 80.7700000000
## range 59.3800000000
## sum 3422.8700000000
## median 67.9650000000
## mean 61.1226785714
## SE.mean 2.2021766902
## CI.mean.0.95 4.4132607078
## var 271.5766017857
## std.dev 16.4795813595
## coef.var 0.2696148426model=(lm(formula=Y~X1+X2+X3+X4,data=data))
model##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = data)
##
## Coefficients:
## (Intercept) X1 X2 X3 X4
## 10.498612415 -0.237421748 0.919838930 -0.008344086 -0.009454529summary(model)##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1295732 -0.9921587 -0.1119246 0.7638639 4.6374335
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.498612415 5.796395481 1.81123 0.0759963 .
## X1 -0.237421748 0.049781592 -4.76927 1.5862e-05 ***
## X2 0.919838930 0.278410032 3.30390 0.0017485 **
## X3 -0.008344086 0.010186319 -0.81915 0.4165145
## X4 -0.009454529 0.012466490 -0.75840 0.4517049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.453091 on 51 degrees of freedom
## Multiple R-squared: 0.5259901, Adjusted R-squared: 0.4888128
## F-statistic: 14.14817 on 4 and 51 DF, p-value: 7.786214e-08vif(model)## X1 X2 X3 X4
## 1.304894529 1.316581210 1.083063970 1.099406028par(mfrow = c(2, 2))
for (j in 2:5) {
plot(
data[, j], data[, 1],
pch = 19,
xlab = names(data)[j],
ylab = names(data)[1],
main = paste(names(data)[1], "vs", names(data)[j])
)
}
par(mfrow = c(1, 1))buat_mx <- function(X, knot, urut = "gcv")
{
X <- as.matrix(X)
knot <- as.matrix(knot)
N <- nrow(X)
m <- ncol(X)
k <- nrow(knot)
if (urut == "gcv")
{
# urutan kolom: b0, x1..xm, lalu per knot: (x1-ks)+ ... (xm-ks)+
trunc <- do.call(
cbind,
lapply(
1:k,
function(s)
sapply(
1:m,
function(j)
pmax(X[, j] - knot[s, j], 0)
)
)
)
mx <- cbind(1, X, trunc)
nama <- c(
"b0",
paste0("x", 1:m),
unlist(
lapply(
1:k,
function(s) paste0("(x", 1:m, "-k", s, ")+")
)
)
)
}
else
{
# urutan kolom: b0, lalu per variabel: xj, (xj-k1)+, (xj-k2)+, ...
kol <- list(rep(1, N))
nama <- "b0"
for (j in 1:m)
{
kol[[length(kol) + 1]] <- X[, j]
nama <- c(nama, paste0("x", j))
for (s in 1:k)
{
kol[[length(kol) + 1]] <- pmax(X[, j] - knot[s, j], 0)
nama <- c(nama, paste0("(x", j, "-k", s, ")+"))
}
}
mx <- do.call(cbind, kol)
}
colnames(mx) <- nama
mx
}GCV untuk k titik knot (k = 1, 2, 3, …)
GCVk <- function(data, k = 1, nk = 50, simpan = TRUE)
{
data <- as.matrix(data)
N <- nrow(data)
M <- ncol(data)
y <- data[, 1]
X <- data[, 2:M, drop = FALSE]
m <- ncol(X)
# kandidat knot tiap variabel
grid <- sapply(
1:m,
function(j) seq(min(X[, j]), max(X[, j]), length.out = nk)
)
grid <- grid[2:(nk - 1), , drop = FALSE]
nk1 <- nrow(grid)
# semua kombinasi indeks knot (k baris x nkomb kolom)
kombinasi <- combn(nk1, k)
nkomb <- ncol(kombinasi)
GCV <- rep(NA_real_, nkomb)
MSE <- rep(NA_real_, nkomb)
Rsq <- rep(NA_real_, nkomb)
for (i in 1:nkomb)
{
knot <- grid[kombinasi[, i], , drop = FALSE]
mx <- buat_mx(X, knot, "gcv")
C <- pinv(crossprod(mx))
B <- C %*% crossprod(mx, y)
yhat <- mx %*% B
SSE <- sum((y - yhat)^2)
SSR <- sum((yhat - mean(y))^2)
MSE[i] <- SSE / N
Rsq[i] <- SSR / (SSR + SSE) * 100
# trace matriks hat = trace(mx C mx')
trA <- sum((mx %*% C) * mx)
GCV[i] <- MSE[i] / ((N - trA) / N)^2
}
# tabel semua kombinasi knot
knotmat <- t(
apply(
kombinasi, 2,
function(idx) as.vector(t(grid[idx, , drop = FALSE]))
)
)
if (nrow(knotmat) != nkomb) knotmat <- t(knotmat)
colnames(knotmat) <- paste0(
"knot", rep(1:k, each = m),
"_x", rep(1:m, times = k)
)
dataAll <- cbind(
GCV = GCV,
Rsq = Rsq,
knot_ke = 1:nkomb,
knotmat
)
if (simpan)
{
write.csv(
dataAll,
file = paste0("dataAll_knot", k, "_v2.csv")
)
}
cat("\n=== ", k, " titik knot: 10 GCV terkecil ===\n", sep = "")
print(head(dataAll[order(GCV), -2], 10))
# knot optimal
best <- which.min(GCV)
knotopt <- grid[kombinasi[, best], , drop = FALSE]
dimnames(knotopt) <- list(paste0("knot_", 1:k), paste0("x", 1:m))
mxopt <- buat_mx(X, knotopt, "gcv")
Bopt <- pinv(crossprod(mxopt)) %*% crossprod(mxopt, y)
rownames(Bopt) <- colnames(mxopt)
colnames(Bopt) <- "Estimasi"
cat("\nKnot optimal ke-", best, " GCV = ", min(GCV), "\n", sep = "")
print(knotopt)
cat("\nEstimasi parameter\n")
print(Bopt)
invisible(
list(
knotopt = knotopt,
mingcv = min(GCV),
Rsq = Rsq[best],
B = Bopt,
dataAll = dataAll
)
)
}Uji signifikansi parameter
uji_k <- function(data, hasil, alpha = 0.1)
{
data <- as.matrix(data)
y <- data[, 1]
X <- data[, 2:ncol(data), drop = FALSE]
n <- nrow(data)
knot <- hasil$knotopt
k <- nrow(knot)
print(knot)
mx <- buat_mx(X, knot, "gcv")
p <- ncol(mx)
C <- pinv(crossprod(mx))
B <- C %*% crossprod(mx, y)
yhat <- mx %*% B
res <- y - yhat
SSE <- sum(res^2)
SST <- sum((y - mean(y))^2)
SSR <- SST - SSE
df_reg <- p - 1
df_error <- n - p
MSR <- SSR / df_reg
MSE <- SSE / df_error
Fhit <- MSR / MSE
pvalue <- pf(Fhit, df_reg, df_error, lower.tail = FALSE)
Rsq <- SSR / SST * 100
SE <- sqrt(diag(MSE * C))
thit <- B[, 1] / SE
pval <- 2 * pt(abs(thit), df_error, lower.tail = FALSE)
tab <- cbind(
Beta = B[, 1],
SE = SE,
t_hitung = thit,
p_value = pval
)
rownames(tab) <- colnames(mx)
print(tab)
cat("\nAnalysis of Variance\n")
cat("Regresi:", df_reg, SSR, MSR, Fhit, "\n")
cat("Error :", df_error, SSE, MSE, "\n")
cat("Total :", n - 1, SST, "\n")
cat("s =", sqrt(MSE), "\n")
cat("R-squared =", Rsq, "\n")
cat("p-value F =", pvalue, "\n")
cat("\nUji simultan:\n")
if (pvalue <= alpha) {
cat("Tolak H0: minimal terdapat satu parameter yang signifikan.\n")
} else {
cat("Gagal menolak H0.\n")
}
cat("\nUji parsial:\n")
for (i in 1:p)
{
if (pval[i] <= alpha) {
cat(colnames(mx)[i], ": signifikan, p-value =", pval[i], "\n")
} else {
cat(colnames(mx)[i], ": tidak signifikan, p-value =", pval[i], "\n")
}
}
write.csv(res, file = paste0("output_uji_residual_knot", k, "_v2.csv"))
write.csv(mx, file = paste0("output_uji_mx_knot", k, "_v2.csv"))
write.csv(yhat, file = paste0("output_uji_yhat_knot", k, "_v2.csv"))
invisible(
list(
B = B,
tab = tab,
Rsq = Rsq,
pvalue = pvalue,
residual = res,
yhat = yhat
)
)
}Pemilihan knot optimal dengan GCV
Catatan: 3 knot memerlukan waktu lebih lama karena jumlah kombinasinya jauh lebih banyak. Kurangi nk (mis. nk = 30) bila terlalu lambat.
hasil1 <- GCVk(data, k = 1)##
## === 1 titik knot: 10 GCV terkecil ===
## GCV knot_ke knot1_x1 knot1_x2 knot1_x3 knot1_x4
## [1,] 1.508186732 45 76.28673469 14.63061224 69.183673469 75.92265306
## [2,] 1.515930572 44 75.92591837 14.55326531 67.719591837 74.71081633
## [3,] 1.603413515 46 76.64755102 14.70795918 70.647755102 77.13448980
## [4,] 1.615274874 43 75.56510204 14.47591837 66.255510204 73.49897959
## [5,] 1.734460815 42 75.20428571 14.39857143 64.791428571 72.28714286
## [6,] 1.763515462 47 77.00836735 14.78530612 72.111836735 78.34632653
## [7,] 1.830094998 3 61.13244898 11.38204082 7.692244898 25.02551020
## [8,] 1.861405941 4 61.49326531 11.45938776 9.156326531 26.23734694
## [9,] 1.869770285 41 74.84346939 14.32122449 63.327346939 71.07530612
## [10,] 1.912401666 5 61.85408163 11.53673469 10.620408163 27.44918367
##
## Knot optimal ke-45 GCV = 1.508186732
## x1 x2 x3 x4
## knot_1 76.28673469 14.63061224 69.18367347 75.92265306
##
## Estimasi parameter
## Estimasi
## b0 10.122393144263
## x1 -0.200697921687
## x2 0.762874046882
## x3 -0.005244855787
## x4 -0.017306260651
## (x1-k1)+ 0.072850412161
## (x2-k1)+ -21.843749932747
## (x3-k1)+ 0.415138942192
## (x4-k1)+ 1.053106626506hasil2 <- GCVk(data, k = 2)##
## === 2 titik knot: 10 GCV terkecil ===
## GCV knot_ke knot1_x1 knot1_x2 knot1_x3 knot1_x4
## [1,] 1.316067641 135 61.13244898 11.38204082 7.692244898 25.02551020
## [2,] 1.317893528 179 61.49326531 11.45938776 9.156326531 26.23734694
## [3,] 1.321113653 178 61.49326531 11.45938776 9.156326531 26.23734694
## [4,] 1.321402209 134 61.13244898 11.38204082 7.692244898 25.02551020
## [5,] 1.358738225 221 61.85408163 11.53673469 10.620408163 27.44918367
## [6,] 1.359629986 222 61.85408163 11.53673469 10.620408163 27.44918367
## [7,] 1.359989299 138 61.13244898 11.38204082 7.692244898 25.02551020
## [8,] 1.388798597 136 61.13244898 11.38204082 7.692244898 25.02551020
## [9,] 1.393179639 180 61.49326531 11.45938776 9.156326531 26.23734694
## [10,] 1.396086664 177 61.49326531 11.45938776 9.156326531 26.23734694
## knot2_x1 knot2_x2 knot2_x3 knot2_x4
## [1,] 76.28673469 14.63061224 69.18367347 75.92265306
## [2,] 76.28673469 14.63061224 69.18367347 75.92265306
## [3,] 75.92591837 14.55326531 67.71959184 74.71081633
## [4,] 75.92591837 14.55326531 67.71959184 74.71081633
## [5,] 75.92591837 14.55326531 67.71959184 74.71081633
## [6,] 76.28673469 14.63061224 69.18367347 75.92265306
## [7,] 77.36918367 14.86265306 73.57591837 79.55816327
## [8,] 76.64755102 14.70795918 70.64775510 77.13448980
## [9,] 76.64755102 14.70795918 70.64775510 77.13448980
## [10,] 75.56510204 14.47591837 66.25551020 73.49897959
##
## Knot optimal ke-135 GCV = 1.316067641
## x1 x2 x3 x4
## knot_1 61.13244898 11.38204082 7.692244898 25.02551020
## knot_2 76.28673469 14.63061224 69.183673469 75.92265306
##
## Estimasi parameter
## Estimasi
## b0 95.3587173978
## x1 -2.9213118661
## x2 6.1469386905
## x3 0.2059982446
## x4 0.7283636092
## (x1-k1)+ 2.7449916135
## (x2-k1)+ -5.3197765896
## (x3-k1)+ -0.2144538392
## (x4-k1)+ -0.7619876831
## (x1-k2)+ 0.7954128596
## (x2-k2)+ -18.0705567159
## (x3-k2)+ 0.2853691492
## (x4-k2)+ 0.9034632270hasil3 <- GCVk(data, k = 3)##
## === 3 titik knot: 10 GCV terkecil ===
## GCV knot_ke knot1_x1 knot1_x2 knot1_x3 knot1_x4
## [1,] 1.444245765 2200 61.13244898 11.38204082 7.692244898 25.02551020
## [2,] 1.445793403 3146 61.49326531 11.45938776 9.156326531 26.23734694
## [3,] 1.447104419 2157 61.13244898 11.38204082 7.692244898 25.02551020
## [4,] 1.452450984 3106 61.13244898 11.38204082 7.692244898 25.02551020
## [5,] 1.453903061 2199 61.13244898 11.38204082 7.692244898 25.02551020
## [6,] 1.454858315 3145 61.49326531 11.45938776 9.156326531 26.23734694
## [7,] 1.458180136 2156 61.13244898 11.38204082 7.692244898 25.02551020
## [8,] 1.458930923 2242 61.13244898 11.38204082 7.692244898 25.02551020
## [9,] 1.464583798 3037 61.13244898 11.38204082 7.692244898 25.02551020
## [10,] 1.465579618 3101 61.13244898 11.38204082 7.692244898 25.02551020
## knot2_x1 knot2_x2 knot2_x3 knot2_x4 knot3_x1 knot3_x2
## [1,] 61.85408163 11.53673469 10.620408163 27.44918367 76.28673469 14.63061224
## [2,] 61.85408163 11.53673469 10.620408163 27.44918367 76.28673469 14.63061224
## [3,] 61.49326531 11.45938776 9.156326531 26.23734694 76.28673469 14.63061224
## [4,] 77.00836735 14.78530612 72.111836735 78.34632653 77.36918367 14.86265306
## [5,] 61.85408163 11.53673469 10.620408163 27.44918367 75.92591837 14.55326531
## [6,] 61.85408163 11.53673469 10.620408163 27.44918367 75.92591837 14.55326531
## [7,] 61.49326531 11.45938776 9.156326531 26.23734694 75.92591837 14.55326531
## [8,] 62.21489796 11.61408163 12.084489796 28.66102041 76.28673469 14.63061224
## [9,] 73.03938776 13.93448980 56.006938776 65.01612245 76.28673469 14.63061224
## [10,] 76.28673469 14.63061224 69.183673469 75.92265306 76.64755102 14.70795918
## knot3_x3 knot3_x4
## [1,] 69.18367347 75.92265306
## [2,] 69.18367347 75.92265306
## [3,] 69.18367347 75.92265306
## [4,] 73.57591837 79.55816327
## [5,] 67.71959184 74.71081633
## [6,] 67.71959184 74.71081633
## [7,] 67.71959184 74.71081633
## [8,] 69.18367347 75.92265306
## [9,] 69.18367347 75.92265306
## [10,] 70.64775510 77.13448980
##
## Knot optimal ke-2200 GCV = 1.444245765
## x1 x2 x3 x4
## knot_1 61.13244898 11.38204082 7.692244898 25.02551020
## knot_2 61.85408163 11.53673469 10.620408163 27.44918367
## knot_3 76.28673469 14.63061224 69.183673469 75.92265306
##
## Estimasi parameter
## Estimasi
## b0 -4.87954771438
## x1 -0.57070096202
## x2 1.29861828003
## x3 0.20162854466
## x4 1.23897093494
## (x1-k1)+ -2.60712938204
## (x2-k1)+ 13.80910091002
## (x3-k1)+ -0.22960609217
## (x4-k1)+ -1.94377922501
## (x1-k2)+ 3.03208075509
## (x2-k2)+ -14.12188191258
## (x3-k2)+ 0.02361871805
## (x4-k2)+ 0.67438515487
## (x1-k3)+ 0.76491110081
## (x2-k3)+ -19.59971920121
## (x3-k3)+ 0.25980650425
## (x4-k3)+ 0.92445502511Perbandingan GCV
gcv_perbandingan <-
data.frame(
Model = c("1 Knot", "2 Knot", "3 Knot"),
GCV = c(hasil1$mingcv, hasil2$mingcv, hasil3$mingcv)
)
gcv_perbandingan## Model GCV
## 1 1 Knot 1.508186732
## 2 2 Knot 1.316067641
## 3 3 Knot 1.444245765ggplot(
gcv_perbandingan,
aes(x = Model, y = GCV, group = 1)
) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
geom_text(aes(label = round(GCV, 5)), vjust = -0.8) +
labs(
title = "Perbandingan Nilai GCV",
x = "Jumlah Titik Knot",
y = "GCV"
) +
theme_minimal()
gcv_minimum <- min(hasil1$mingcv, hasil2$mingcv, hasil3$mingcv)
model_terbaik <-
c("1 Knot", "2 Knot", "3 Knot")[
which.min(c(hasil1$mingcv, hasil2$mingcv, hasil3$mingcv))
]
cat("Model dengan GCV minimum:", model_terbaik, "\n")## Model dengan GCV minimum: 2 Knotcat("Nilai GCV:", gcv_minimum, "\n")## Nilai GCV: 1.316067641knot_optimal <- hasil2$knotopt
knot_optimal## x1 x2 x3 x4
## knot_1 61.13244898 11.38204082 7.692244898 25.02551020
## knot_2 76.28673469 14.63061224 69.183673469 75.92265306B_optimal <- hasil2$B
B_optimal## Estimasi
## b0 95.3587173978
## x1 -2.9213118661
## x2 6.1469386905
## x3 0.2059982446
## x4 0.7283636092
## (x1-k1)+ 2.7449916135
## (x2-k1)+ -5.3197765896
## (x3-k1)+ -0.2144538392
## (x4-k1)+ -0.7619876831
## (x1-k2)+ 0.7954128596
## (x2-k2)+ -18.0705567159
## (x3-k2)+ 0.2853691492
## (x4-k2)+ 0.9034632270Uji signifikansi parameter
uji2 <- uji_k(data, hasil2, alpha = 0.1)## x1 x2 x3 x4
## knot_1 61.13244898 11.38204082 7.692244898 25.02551020
## knot_2 76.28673469 14.63061224 69.183673469 75.92265306
## Beta SE t_hitung p_value
## b0 95.3587173978 84.8285642104 1.124134521 2.671924603e-01
## x1 -2.9213118661 1.0418248755 -2.804033513 7.544198504e-03
## x2 6.1469386905 4.8596442762 1.264894783 2.127214042e-01
## x3 0.2059982446 0.1626993431 1.266128312 2.122837024e-01
## x4 0.7283636092 0.3657961901 1.991173306 5.283963215e-02
## (x1-k1)+ 2.7449916135 1.0529267776 2.607010926 1.250391458e-02
## (x2-k1)+ -5.3197765896 4.9095217464 -1.083563097 2.845971079e-01
## (x3-k1)+ -0.2144538392 0.1669139901 -1.284816444 2.057346472e-01
## (x4-k1)+ -0.7619876831 0.3717972255 -2.049471139 4.654749657e-02
## (x1-k2)+ 0.7954128596 0.7655042699 1.039070441 3.045796347e-01
## (x2-k2)+ -18.0705567159 4.2319010493 -4.270080161 1.058266457e-04
## (x3-k2)+ 0.2853691492 0.2305324459 1.237869785 2.224808174e-01
## (x4-k2)+ 0.9034632270 0.1935247451 4.668463594 2.974409879e-05
##
## Analysis of Variance
## Regresi: 12 183.7252095 15.31043413 15.15056629
## Error : 43 43.45373333 1.010551938
## Total : 55 227.1789429
## s = 1.005262124
## R-squared = 80.87246433
## p-value F = 9.568367676e-12
##
## Uji simultan:
## Tolak H0: minimal terdapat satu parameter yang signifikan.
##
## Uji parsial:
## b0 : tidak signifikan, p-value = 0.2671924603
## x1 : signifikan, p-value = 0.007544198504
## x2 : tidak signifikan, p-value = 0.2127214042
## x3 : tidak signifikan, p-value = 0.2122837024
## x4 : signifikan, p-value = 0.05283963215
## (x1-k1)+ : signifikan, p-value = 0.01250391458
## (x2-k1)+ : tidak signifikan, p-value = 0.2845971079
## (x3-k1)+ : tidak signifikan, p-value = 0.2057346472
## (x4-k1)+ : signifikan, p-value = 0.04654749657
## (x1-k2)+ : tidak signifikan, p-value = 0.3045796347
## (x2-k2)+ : signifikan, p-value = 0.0001058266457
## (x3-k2)+ : tidak signifikan, p-value = 0.2224808174
## (x4-k2)+ : signifikan, p-value = 2.974409879e-05uji3 <- uji_k(data, hasil3, alpha = 0.1)## x1 x2 x3 x4
## knot_1 61.13244898 11.38204082 7.692244898 25.02551020
## knot_2 61.85408163 11.53673469 10.620408163 27.44918367
## knot_3 76.28673469 14.63061224 69.183673469 75.92265306
## Beta SE t_hitung p_value
## b0 -4.87954771438 11.3313795043 -0.43062256564 6.691134312e-01
## x1 -0.57070096202 2.3596927737 -0.24185392623 8.101605857e-01
## x2 1.29861828003 12.0225195725 0.10801548479 9.145372133e-01
## x3 0.20162854466 0.3213964214 0.62735155478 5.340847884e-01
## x4 1.23897093494 0.7333330296 1.68950652009 9.910266384e-02
## (x1-k1)+ -2.60712938204 4.8454014797 -0.53806261317 5.935914409e-01
## (x2-k1)+ 13.80910091002 33.0174968574 0.41823585142 6.780687838e-01
## (x3-k1)+ -0.22960609217 0.6222171416 -0.36901280408 7.141140668e-01
## (x4-k1)+ -1.94377922501 1.4399710751 -1.34987379860 1.848396688e-01
## (x1-k2)+ 3.03208075509 2.6352514927 1.15058496825 2.569090837e-01
## (x2-k2)+ -14.12188191258 22.0729359865 -0.63978266966 5.260546626e-01
## (x3-k2)+ 0.02361871805 0.3275433283 0.07210868307 9.428839547e-01
## (x4-k2)+ 0.67438515487 0.7687876279 0.87720604546 3.857465713e-01
## (x1-k3)+ 0.76491110081 0.8752551052 0.87392932213 3.875065807e-01
## (x2-k3)+ -19.59971920121 4.6428329720 -4.22149995909 1.404269405e-04
## (x3-k3)+ 0.25980650425 0.2406862162 1.07944072705 2.870220588e-01
## (x4-k3)+ 0.92445502511 0.2015337913 4.58709687949 4.564589289e-05
##
## Analysis of Variance
## Regresi: 16 185.9147781 11.61967363 10.98210214
## Error : 39 41.26416471 1.058055505
## Total : 55 227.1789429
## s = 1.028618251
## R-squared = 81.83627224
## p-value F = 7.291164671e-10
##
## Uji simultan:
## Tolak H0: minimal terdapat satu parameter yang signifikan.
##
## Uji parsial:
## b0 : tidak signifikan, p-value = 0.6691134312
## x1 : tidak signifikan, p-value = 0.8101605857
## x2 : tidak signifikan, p-value = 0.9145372133
## x3 : tidak signifikan, p-value = 0.5340847884
## x4 : signifikan, p-value = 0.09910266384
## (x1-k1)+ : tidak signifikan, p-value = 0.5935914409
## (x2-k1)+ : tidak signifikan, p-value = 0.6780687838
## (x3-k1)+ : tidak signifikan, p-value = 0.7141140668
## (x4-k1)+ : tidak signifikan, p-value = 0.1848396688
## (x1-k2)+ : tidak signifikan, p-value = 0.2569090837
## (x2-k2)+ : tidak signifikan, p-value = 0.5260546626
## (x3-k2)+ : tidak signifikan, p-value = 0.9428839547
## (x4-k2)+ : tidak signifikan, p-value = 0.3857465713
## (x1-k3)+ : tidak signifikan, p-value = 0.3875065807
## (x2-k3)+ : signifikan, p-value = 0.0001404269405
## (x3-k3)+ : tidak signifikan, p-value = 0.2870220588
## (x4-k3)+ : signifikan, p-value = 4.564589289e-05Plot nilai aktual vs prediksi (model 2 knot)
X <- as.matrix(data[, 2:5])
y <- data[, 1]
knot <- hasil2$knotopt
mx <- buat_mx(X, knot, "gcv")
yhat <- as.vector(mx %*% hasil2$B)
data_plot <-
data.frame(
Observasi = 1:nrow(data),
Aktual = y,
Prediksi = yhat
)
ggplot(data_plot, aes(x = Observasi)) +
geom_point(aes(y = Aktual), size = 2) +
geom_line(aes(y = Aktual, colour = "Aktual"), linewidth = 0.7) +
geom_line(aes(y = Prediksi, colour = "Prediksi"), linewidth = 1) +
scale_colour_manual(
name = NULL,
values = c(Aktual = "black", Prediksi = "red")
) +
labs(
title = "Nilai Aktual dan Prediksi Model Spline 2 Knot",
x = "Pengamatan",
y = "Nilai"
) +
theme_minimal()