Memanggil Data

File data_spline.csv harus berada satu folder dengan file Rmd ini. Jika ingin memakai file asli (.txt), ganti baris read.csv(...) dengan read.table(file.choose(), header = TRUE).

data <- read.table(file.choose(),header=TRUE)  # kolom: X1, X2, X3, X4, X5, X6, Y

# PENTING: fungsi GCV dan uji di bawah menganggap KOLOM 1 = respons (Y)
# dan kolom 2 dst = prediktor. Urutan kolom file asli (Y di kolom terakhir)
# harus disusun ulang agar Y menjadi kolom pertama.
data <- data[, c("Y", "X1", "X2", "X3", "X4", "X5", "X6")]
stopifnot(colnames(data)[1] == "Y", nrow(data) == 514)
head(data)

Analisis Statistika Deskriptif

round(stat.desc(data), 4)

Tabel 1. Analisis Statistika Deskriptif

tabel1 <- data.frame(
  Variabel = colnames(data),
  Data     = nrow(data),
  Rata2    = sapply(data, mean),
  Minimum  = sapply(data, min),
  Maksimum = sapply(data, max),
  SD       = sapply(data, sd),
  row.names = NULL
)
knitr::kable(tabel1, digits = 4)

Variabel	Data	Rata2	Minimum	Maksimum	SD
Y	514	75.1104	14.14	96.37	14.4301
X1	514	11.4999	2.27	40.01	7.1544
X2	514	28.6938	0.00	99.62	20.0524
X3	514	8.8510	1.28	12.74	1.5454
X4	514	4.0321	0.00	86.82	8.7419
X5	514	70.2002	55.72	77.93	3.3852
X6	514	22.3482	0.00	50.20	9.1727

Scatterplot

for (v in paste0("X", 1:6)) {
  plot(data[[v]], data$Y,
       main = paste("Scatterplot Y dan", v),
       xlab = paste0("(", v, ")"), ylab = "(Y)", col = "black")
  abline(lm(data$Y ~ data[[v]]), col = "yellow", lwd = 3)
}

Interpretasi:

for (v in paste0("X", 1:6)) {
  r    <- cor(data[[v]], data$Y)
  arah <- ifelse(r < 0, "negatif", "positif")
  kuat <- ifelse(abs(r) >= 0.7, "kuat", ifelse(abs(r) >= 0.3, "sedang", "lemah"))
  efek <- ifelse(r < 0, "semakin kecil", "semakin besar")
  cat(sprintf(
    "- Hubungan antara variabel %s dan Y bersifat %s dan %s (r = %s). Artinya, semakin besar nilai %s maka cenderung %s nilai Y.\n",
    v, arah, kuat, fmt(r, 4), v, efek))
}

Hubungan antara variabel X1 dan Y bersifat negatif dan sedang (r = -0,6811). Artinya, semakin besar nilai X1 maka cenderung semakin kecil nilai Y.
Hubungan antara variabel X2 dan Y bersifat negatif dan sedang (r = -0,6669). Artinya, semakin besar nilai X2 maka cenderung semakin kecil nilai Y.
Hubungan antara variabel X3 dan Y bersifat positif dan sedang (r = 0,4573). Artinya, semakin besar nilai X3 maka cenderung semakin besar nilai Y.
Hubungan antara variabel X4 dan Y bersifat negatif dan sedang (r = -0,5922). Artinya, semakin besar nilai X4 maka cenderung semakin kecil nilai Y.
Hubungan antara variabel X5 dan Y bersifat positif dan sedang (r = 0,6291). Artinya, semakin besar nilai X5 maka cenderung semakin besar nilai Y.
Hubungan antara variabel X6 dan Y bersifat negatif dan lemah (r = -0,1121). Artinya, semakin besar nilai X6 maka cenderung semakin kecil nilai Y.

Model Regresi Linear Awal

\[Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{3i} + \beta_4 X_{4i} + \beta_5 X_{5i} + \beta_6 X_{6i} + \varepsilon_i,\quad i = 1,2,\ldots,514\]

model <- lm(Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
model

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
## 
## Coefficients:
## (Intercept)           X1           X2           X3           X4           X5  
##    26.47106     -0.61512     -0.21398     -0.21664     -0.42017      0.92267  
##          X6  
##     0.03101

summary(model)

## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.041  -3.610   2.047   5.431  21.434 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.47106   11.16110   2.372   0.0181 *  
## X1          -0.61512    0.07038  -8.740  < 2e-16 ***
## X2          -0.21398    0.02453  -8.721  < 2e-16 ***
## X3          -0.21664    0.29665  -0.730   0.4656    
## X4          -0.42017    0.04933  -8.518  < 2e-16 ***
## X5           0.92267    0.14182   6.506 1.85e-10 ***
## X6           0.03101    0.04248   0.730   0.4657    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.209 on 507 degrees of freedom
## Multiple R-squared:  0.6801, Adjusted R-squared:  0.6763 
## F-statistic: 179.7 on 6 and 507 DF,  p-value: < 2.2e-16

Estimasi Parameter

b <- coef(model)
cat(sprintf("Konstanta sebesar %s menyatakan nilai Y ketika seluruh X bernilai nol.\n\n", fmt(b[1])))

Konstanta sebesar 26,47106 menyatakan nilai Y ketika seluruh X bernilai nol.

for (v in paste0("X", 1:6)) {
  arah <- ifelse(b[v] < 0, "menurunkan", "meningkatkan")
  cat(sprintf("- Koefisien %s sebesar %s: setiap kenaikan satu satuan %s %s Y sebesar %s satuan.\n",
              v, fmt(b[v]), v, arah, fmt(abs(b[v]))))
}

Koefisien X1 sebesar -0,61512: setiap kenaikan satu satuan X1 menurunkan Y sebesar 0,61512 satuan.
Koefisien X2 sebesar -0,21398: setiap kenaikan satu satuan X2 menurunkan Y sebesar 0,21398 satuan.
Koefisien X3 sebesar -0,21664: setiap kenaikan satu satuan X3 menurunkan Y sebesar 0,21664 satuan.
Koefisien X4 sebesar -0,42017: setiap kenaikan satu satuan X4 menurunkan Y sebesar 0,42017 satuan.
Koefisien X5 sebesar 0,92267: setiap kenaikan satu satuan X5 meningkatkan Y sebesar 0,92267 satuan.
Koefisien X6 sebesar 0,03101: setiap kenaikan satu satuan X6 meningkatkan Y sebesar 0,03101 satuan.

Koefisien Determinasi

r2 <- summary(model)$r.squared
cat(sprintf("Nilai R² = %s, artinya %s%% variasi Y dijelaskan oleh X1–X6, sedangkan %s%% dijelaskan oleh faktor lain.\n",
            fmt(r2, 4), fmt(100 * r2, 2), fmt(100 * (1 - r2), 2)))

Nilai R² = 0,6801, artinya 68,01% variasi Y dijelaskan oleh X1–X6, sedangkan 31,99% dijelaskan oleh faktor lain.

Identifikasi Multikolinearitas

vif(model)

##       X1       X2       X3       X4       X5       X6 
## 1.930060 1.842417 1.599757 1.415399 1.754443 1.155668

Seluruh nilai VIF berada di bawah 10, sehingga tidak terdapat multikolinearitas dan analisis dapat dilanjutkan.

Model Regresi Nonparametrik Spline Truncated

Model dengan 6 prediktor dan \(r\) titik knot pada tiap prediktor:

\[y_i = \beta_0 + \sum_{j=1}^{6}\beta_{j1}x_{ji} + \sum_{j=1}^{6}\sum_{v=1}^{r}\beta_{j(1+v)}(x_{ji}-k_{vj})_+ + \varepsilon_i\]

dengan \((x_{ji}-k_{vj})_+ = x_{ji}-k_{vj}\) jika \(x_{ji} \ge k_{vj}\) dan \(0\) jika \(x_{ji} < k_{vj}\).

Fungsi Pembantu

# Grid kandidat knot: 50 titik sama jarak, titik pertama dan terakhir dibuang
buat_grid <- function(x, nk = 50) {
  g <- seq(min(x), max(x), length.out = nk)
  g[2:(nk - 1)]
}

# Matriks desain: 1, lalu per variabel [x_j, (x_j - knot_1)+, ..., (x_j - knot_r)+]
# X: matriks n x p ; knot: matriks p x r
bangun_mx <- function(X, knot) {
  n <- nrow(X); p <- ncol(X); r <- ncol(knot)
  mx <- matrix(1, nrow = n, ncol = 1)
  nama <- "satu"
  for (j in 1:p) {
    mx <- cbind(mx, X[, j])
    nama <- c(nama, colnames(X)[j])
    for (v in 1:r) {
      mx <- cbind(mx, pmax(X[, j] - knot[j, v], 0))
      nama <- c(nama, paste0(colnames(X)[j], "_knot", v))
    }
  }
  colnames(mx) <- nama
  mx
}

# Pemilihan knot optimal berdasarkan GCV minimum untuk r titik knot
# (kombinasi indeks grid yang sama dipakai untuk semua variabel)
gcv_knot <- function(data, r, nk = 50) {
  y <- data[, 1]
  X <- as.matrix(data[, -1])
  n <- nrow(X); p <- ncol(X)
  grid <- sapply(1:p, function(j) buat_grid(X[, j], nk))     # (nk-2) x p
  kombinasi <- matrix(combn(nk - 2, r), nrow = r)             # r x banyak kombinasi
  nc <- ncol(kombinasi)

  GCV <- rep(NA_real_, nc)
  Rsq <- rep(NA_real_, nc)
  for (i in 1:nc) {
    knot <- matrix(NA_real_, nrow = p, ncol = r)
    idx <- kombinasi[, i]                      # indeks grid untuk kombinasi ke-i
    for (j in 1:p) knot[j, ] <- grid[idx, j]
    mx  <- bangun_mx(X, knot)
    fit <- lm.fit(mx, y)
    sse <- sum(fit$residuals^2)
    ssr <- sum((fit$fitted.values - mean(y))^2)
    GCV[i] <- (sse / n) / (((n - fit$rank) / n)^2)
    Rsq[i] <- 100 * ssr / (ssr + sse)
  }

  # tabel semua kombinasi knot: GCV, Rsq, knot_ke, lalu knot per variabel
  kolom_knot <- NULL
  nama_knot  <- character(0)
  for (j in 1:p) {
    for (v in 1:r) {
      kolom_knot <- cbind(kolom_knot, grid[kombinasi[v, ], j])
      nama_knot  <- c(nama_knot, paste0(colnames(X)[j], "_k", v))
    }
  }
  colnames(kolom_knot) <- nama_knot
  dataAll <- data.frame(GCV = GCV, Rsq = Rsq, knot_ke = 1:nc, kolom_knot)
  write.csv(dataAll, file = file.path(dir_out, paste0("Data Knot ", r, ".csv")))

  best <- which.min(GCV)
  knot_terbaik <- matrix(kolom_knot[best, ], nrow = p, byrow = TRUE)
  rownames(knot_terbaik) <- colnames(X)
  colnames(knot_terbaik) <- paste0("knot", 1:r)

  list(tabel = dataAll, best = best, mingcv = GCV[best],
       knot = knot_terbaik, top10 = dataAll[order(dataAll$GCV)[1:10], ])
}

# Estimasi parameter + uji simultan + uji parsial untuk knot terbaik
uji_spline <- function(data, knot, alpha = 0.05) {
  y <- data[, 1]
  X <- as.matrix(data[, -1])
  n <- nrow(X); r <- ncol(knot)
  mx <- bangun_mx(X, knot)
  k  <- ncol(mx)

  XtX <- crossprod(mx)
  B   <- drop(solve(XtX, crossprod(mx, y)))
  names(B) <- colnames(mx)
  yhat <- drop(mx %*% B)
  res  <- y - yhat

  SSE <- sum(res^2)
  SSR <- sum((yhat - mean(y))^2)
  SST <- sum((y - mean(y))^2)
  MSE <- SSE / (n - k)        # kuadrat tengah galat (db = n - k)
  MSR <- SSR / (k - 1)
  Fhit   <- MSR / MSE
  pvalue <- pf(Fhit, k - 1, n - k, lower.tail = FALSE)
  Rsq    <- 100 * SSR / SST

  SE   <- sqrt(diag(MSE * solve(XtX)))
  thit <- B / SE
  pval <- 2 * pt(abs(thit), n - k, lower.tail = FALSE)

  parsial <- data.frame(
    Parameter = colnames(mx), Estimasi = B, SE = SE, t_hitung = thit,
    p_value = pval,
    Keputusan = ifelse(pval < alpha, "Tolak H0", "Gagal Tolak H0"),
    row.names = NULL
  )
  anova <- data.frame(
    Sumber = c("Regresi", "Galat", "Total"),
    db = c(k - 1, n - k, n - 1),
    JK = c(SSR, SSE, SST),
    KT = c(MSR, MSE, NA),
    F_hitung = c(Fhit, NA, NA),
    p_value = c(pvalue, NA, NA)
  )

  # keluaran CSV (sesuai nama berkas pada laporan)
  write.csv(mx,                  file.path(dir_out, paste0("Output Uji mx Knot ", r, ".csv")))
  write.csv(data.frame(x = res), file.path(dir_out, paste0("Output Uji Residual Knot ", r, ".csv")))
  write.csv(data.frame(x = yhat), file.path(dir_out, paste0("Output Uji yhat Knot ", r, ".csv")))

  list(B = B, mx = mx, yhat = yhat, res = res, parsial = parsial, anova = anova,
       Rsq = Rsq, pvalue = pvalue, k = k)
}

Pemilihan Titik Knot Optimal: Satu Titik Knot

g1 <- gcv_knot(data, r = 1)
cat("GCV minimum 1 knot =", fmt(g1$mingcv), "pada knot ke-", g1$best, "\n")

## GCV minimum 1 knot = 60,89140 pada knot ke- 15

g1$knot

##        knot1
## X1 13.823061
## X2 30.495918
## X3  4.788163
## X4 26.577551
## X5 62.518980
## X6 15.367347

Tabel 5. Nilai GCV dengan Satu Titik Knot (10 terkecil)

knitr::kable(g1$top10, digits = 5, row.names = FALSE)

GCV	Rsq	knot_ke	X1_k1	X2_k1	X3_k1	X4_k1	X5_k1	X6_k1
60.89140	72.16346	15	13.82306	30.49592	4.78816	26.57755	62.51898	15.36735
60.90658	72.15653	16	14.59327	32.52898	5.02204	28.34939	62.97224	16.39184
60.93708	72.14258	17	15.36347	34.56204	5.25592	30.12122	63.42551	17.41633
61.00914	72.10964	18	16.13367	36.59510	5.48980	31.89306	63.87878	18.44082
61.03068	72.09979	19	16.90388	38.62816	5.72367	33.66490	64.33204	19.46531
61.03475	72.09793	14	13.05286	28.46286	4.55429	24.80571	62.06571	14.34286
61.08612	72.07445	20	17.67408	40.66122	5.95755	35.43673	64.78531	20.48980
61.22630	72.01037	21	18.44429	42.69429	6.19143	37.20857	65.23857	21.51429
61.33258	71.96178	22	19.21449	44.72735	6.42531	38.98041	65.69184	22.53878
61.42555	71.91928	13	12.28265	26.42980	4.32041	23.03388	61.61245	13.31837

u1 <- uji_spline(data, g1$knot, alpha)
u1$B

##         satu           X1     X1_knot1           X2     X2_knot1           X3 
## 55.235597561 -0.060206537 -0.866308131 -0.139780329 -0.081699055  4.602967603 
##     X3_knot1           X4     X4_knot1           X5     X5_knot1           X6 
## -5.057603076 -0.613837239  0.720183134  0.008921615  1.020919185  0.059789517 
##     X6_knot1 
## -0.129832245

Pemilihan Titik Knot Optimal: Dua Titik Knot

g2 <- gcv_knot(data, r = 2)
cat("GCV minimum 2 knot =", fmt(g2$mingcv), "pada knot ke-", g2$best, "\n")

## GCV minimum 2 knot = 59,65374 pada knot ke- 245

g2$knot

##        knot1     knot2
## X1  6.891224 22.295306
## X2 12.198367 52.859592
## X3  2.683265  7.360816
## X4 10.631020 46.067755
## X5 58.439592 67.504898
## X6  6.146939 26.636735

Tabel 6. Nilai GCV dengan Dua Titik Knot (10 terkecil)

knitr::kable(g2$top10, digits = 5, row.names = FALSE)

GCV	Rsq	knot_ke	X1_k1	X1_k2	X2_k1	X2_k2	X3_k1	X3_k2	X4_k1	X4_k2	X5_k1	X5_k2	X6_k1	X6_k2
59.65374	73.37854	245	6.89122	22.29531	12.19837	52.85959	2.68327	7.36082	10.63102	46.06776	58.43959	67.50490	6.14694	26.63673
59.71889	73.34947	239	6.89122	17.67408	12.19837	40.66122	2.68327	5.95755	10.63102	35.43673	58.43959	64.78531	6.14694	20.48980
59.71968	73.34911	246	6.89122	23.06551	12.19837	54.89265	2.68327	7.59469	10.63102	47.83959	58.43959	67.95816	6.14694	27.66122
59.72621	73.34620	244	6.89122	21.52510	12.19837	50.82653	2.68327	7.12694	10.63102	44.29592	58.43959	67.05163	6.14694	25.61224
59.72924	73.34485	238	6.89122	16.90388	12.19837	38.62816	2.68327	5.72367	10.63102	33.66490	58.43959	64.33204	6.14694	19.46531
59.73756	73.34114	905	23.06551	29.22714	54.89265	71.15714	7.59469	9.46571	47.83959	62.01429	67.95816	71.58429	27.66122	35.85714
59.76477	73.32900	237	6.89122	16.13367	12.19837	36.59510	2.68327	5.48980	10.63102	31.89306	58.43959	63.87878	6.14694	18.44082
59.76563	73.32861	241	6.89122	19.21449	12.19837	44.72735	2.68327	6.42531	10.63102	38.98041	58.43959	65.69184	6.14694	22.53878
59.76975	73.32677	240	6.89122	18.44429	12.19837	42.69429	2.68327	6.19143	10.63102	37.20857	58.43959	65.23857	6.14694	21.51429
59.77079	73.32631	242	6.89122	19.98469	12.19837	46.76041	2.68327	6.65918	10.63102	40.75224	58.43959	66.14510	6.14694	23.56327

u2 <- uji_spline(data, g2$knot, alpha)
u2$B

##         satu           X1     X1_knot1     X1_knot2           X2     X2_knot1 
## 232.61514822  -0.21284620  -0.04625130  -1.12800600  -0.47339109   0.33226449 
##     X2_knot2           X3     X3_knot1     X3_knot2           X4     X4_knot1 
##  -0.04946566  -3.37097474   6.15885749  -4.00900376  -0.93861870   0.87768543 
##     X4_knot2           X5     X5_knot1     X5_knot2           X6     X6_knot1 
##   0.02702941  -2.77940374   4.11508131  -0.43791153   0.78647431  -0.91535404 
##     X6_knot2 
##   0.05822887

Pemilihan Titik Knot Optimal: Tiga Titik Knot

Proses ini mengevaluasi 17.296 kombinasi knot sehingga memerlukan waktu sekitar satu menit.

g3 <- gcv_knot(data, r = 3)
cat("GCV minimum 3 knot =", fmt(g3$mingcv), "pada knot ke-", g3$best, "\n")

## GCV minimum 3 knot = 57,87155 pada knot ke- 11127

g3$knot

##        knot1     knot2     knot3
## X1 13.052857 24.605918 28.456939
## X2 28.462857 58.958776 69.124082
## X3  4.554286  8.062449  9.231837
## X4 24.805714 51.383265 60.242449
## X5 62.065714 68.864694 71.131020
## X6 14.342857 29.710204 34.832653

Tabel 7. Nilai GCV dengan Tiga Titik Knot (10 terkecil)

knitr::kable(g3$top10, digits = 5, row.names = FALSE)

GCV	Rsq	knot_ke	X1_k1	X1_k2	X1_k3	X2_k1	X2_k2	X2_k3	X3_k1	X3_k2	X3_k3	X4_k1	X4_k2	X4_k3	X5_k1	X5_k2	X5_k3	X6_k1	X6_k2	X6_k3
57.87155	74.79617	11127	13.05286	24.60592	28.45694	28.46286	58.95878	69.12408	4.55429	8.06245	9.23184	24.80571	51.38327	60.24245	62.06571	68.86469	71.13102	14.34286	29.71020	34.83265
57.90291	74.78251	11144	13.05286	25.37612	27.68673	28.46286	60.99184	67.09102	4.55429	8.29633	8.99796	24.80571	53.15510	58.47061	62.06571	69.31796	70.67776	14.34286	30.73469	33.80816
57.91615	74.77674	11126	13.05286	24.60592	27.68673	28.46286	58.95878	67.09102	4.55429	8.06245	8.99796	24.80571	51.38327	58.47061	62.06571	68.86469	70.67776	14.34286	29.71020	33.80816
57.97474	74.75123	11655	13.82306	24.60592	28.45694	30.49592	58.95878	69.12408	4.78816	8.06245	9.23184	26.57755	51.38327	60.24245	62.51898	68.86469	71.13102	15.36735	29.71020	34.83265
58.00390	74.73853	11672	13.82306	25.37612	27.68673	30.49592	60.99184	67.09102	4.78816	8.29633	8.99796	26.57755	53.15510	58.47061	62.51898	69.31796	70.67776	15.36735	30.73469	33.80816
58.01200	74.73500	11654	13.82306	24.60592	27.68673	30.49592	58.95878	67.09102	4.78816	8.06245	8.99796	26.57755	51.38327	58.47061	62.51898	68.86469	70.67776	15.36735	29.71020	33.80816
58.01834	74.73224	10566	12.28265	24.60592	28.45694	26.42980	58.95878	69.12408	4.32041	8.06245	9.23184	23.03388	51.38327	60.24245	61.61245	68.86469	71.13102	13.31837	29.71020	34.83265
58.05155	74.71777	11145	13.05286	25.37612	28.45694	28.46286	60.99184	69.12408	4.55429	8.29633	9.23184	24.80571	53.15510	60.24245	62.06571	69.31796	71.13102	14.34286	30.73469	34.83265
58.07181	74.70895	10583	12.28265	25.37612	27.68673	26.42980	60.99184	67.09102	4.32041	8.29633	8.99796	23.03388	53.15510	58.47061	61.61245	69.31796	70.67776	13.31837	30.73469	33.80816
58.08304	74.70406	10565	12.28265	24.60592	27.68673	26.42980	58.95878	67.09102	4.32041	8.06245	8.99796	23.03388	51.38327	58.47061	61.61245	68.86469	70.67776	13.31837	29.71020	33.80816

u3 <- uji_spline(data, g3$knot, alpha)
u3$B

##          satu            X1      X1_knot1      X1_knot2      X1_knot3 
## 167.379165414   0.006721307  -0.884363855   0.860349570  -2.205508265 
##            X2      X2_knot1      X2_knot2      X2_knot3            X3 
##  -0.162048384  -0.027020113  -0.313063519   0.513304370   3.088515428 
##      X3_knot1      X3_knot2      X3_knot3            X4      X4_knot1 
##  -1.193279338  -6.518365931   6.397219439  -0.572073896   0.756445243 
##      X4_knot2      X4_knot3            X5      X5_knot1      X5_knot2 
##   0.935416562  -2.197128367  -1.828642799   3.616134322  -2.564442077 
##      X5_knot3            X6      X6_knot1      X6_knot2      X6_knot3 
##   2.261039366   0.121548312  -0.235273051  -0.185090874   0.665478503

Pemilihan Model Terbaik Berdasarkan GCV Minimum

Tabel 8. Perbandingan Nilai GCV Minimum dari 1, 2, 3 Titik Knot

tabel8 <- data.frame(
  Titik_Knot = c("1 Titik Knot", "2 Titik Knot", "3 Titik Knot"),
  GCV = c(g1$mingcv, g2$mingcv, g3$mingcv),
  R2_persen = c(u1$Rsq, u2$Rsq, u3$Rsq)
)
knitr::kable(tabel8, digits = 5)

Titik_Knot	GCV	R2_persen
1 Titik Knot	60.89140	72.16346
2 Titik Knot	59.65374	73.37854
3 Titik Knot	57.87155	74.79617

idx <- which.min(tabel8$GCV)
cat(sprintf("Nilai GCV terkecil adalah %s, yaitu pada model dengan **%s**. Model inilah yang digunakan sebagai model terbaik.\n",
            fmt(tabel8$GCV[idx]), tolower(tabel8$Titik_Knot[idx])))

Nilai GCV terkecil adalah 57,87155, yaitu pada model dengan 3 titik knot. Model inilah yang digunakan sebagai model terbaik.

# objek model terbaik (otomatis mengikuti GCV terkecil)
hasil <- list(g1, g2, g3)[[idx]]
uji   <- list(u1, u2, u3)[[idx]]
r_opt <- idx

Penaksiran Parameter Model Terbaik

knitr::kable(round(hasil$knot, 5))

	knot1	knot2	knot3
X1	13.05286	24.60592	28.45694
X2	28.46286	58.95878	69.12408
X3	4.55429	8.06245	9.23184
X4	24.80571	51.38327	60.24245
X5	62.06571	68.86469	71.13102
X6	14.34286	29.71020	34.83265

B <- uji$B
vars <- paste0("X", 1:6)
suku <- character(0)
tanda <- function(x) ifelse(x < 0, "-", "+")
for (v in vars) {
  suku <- c(suku, sprintf("%s %s·%s", tanda(B[v]), fmt(abs(B[v])), tolower(v)))
  for (m in 1:r_opt) {
    suku <- c(suku, sprintf("%s %s·(%s − %s)+", tanda(B[paste0(v, "_knot", m)]),
                            fmt(abs(B[paste0(v, "_knot", m)])), tolower(v),
                            fmt(hasil$knot[v, m])))
  }
}
cat("ŷ = ", fmt(B["satu"]), " ", paste(suku, collapse = " "), "\n\n", sep = "")

ŷ = 167,37917 + 0,00672·x1 - 0,88436·(x1 − 13,05286)+ + 0,86035·(x1 − 24,60592)+ - 2,20551·(x1 − 28,45694)+ - 0,16205·x2 - 0,02702·(x2 − 28,46286)+ - 0,31306·(x2 − 58,95878)+ + 0,51330·(x2 − 69,12408)+ + 3,08852·x3 - 1,19328·(x3 − 4,55429)+ - 6,51837·(x3 − 8,06245)+ + 6,39722·(x3 − 9,23184)+ - 0,57207·x4 + 0,75645·(x4 − 24,80571)+ + 0,93542·(x4 − 51,38327)+ - 2,19713·(x4 − 60,24245)+ - 1,82864·x5 + 3,61613·(x5 − 62,06571)+ - 2,56444·(x5 − 68,86469)+ + 2,26104·(x5 − 71,13102)+ + 0,12155·x6 - 0,23527·(x6 − 14,34286)+ - 0,18509·(x6 − 29,71020)+ + 0,66548·(x6 − 34,83265)+

cat(sprintf("Model memiliki %d parameter dengan R² = %s%%.\n", uji$k, fmt(uji$Rsq, 5)))

Model memiliki 25 parameter dengan R² = 74,79617%.

Pengujian Signifikansi Parameter Secara Simultan

Hipotesis

\(H_0: \beta_{1,1} = \beta_{1,2} = \cdots = \beta_{6,r} = 0\) (secara simultan X1–X6 tidak berpengaruh signifikan terhadap Y)

\(H_1:\) minimal terdapat satu \(\beta_{j(1+h)} \neq 0\)

Taraf signifikansi \(\alpha = 0{,}05\) ; Statistik uji \(F = KT_{regresi}/KT_{galat}\) ; Daerah kritis \(H_0\) ditolak jika \(p\text{-value} < \alpha\).

Tabel 9. Analisis Varians Model Regresi Secara Simultan

knitr::kable(uji$anova, digits = 5)

Sumber	db	JK	KT	F_hitung	p_value
Regresi	24	79897.37	3329.05694	60.46587	0
Galat	489	26922.77	55.05679	NA	NA
Total	513	106820.14	NA	NA	NA

if (uji$pvalue < alpha) {
  cat(sprintf("**Keputusan:** H0 ditolak karena p-value = %s < α = %s.\n\n**Kesimpulan:** secara simultan X1–X6 berpengaruh signifikan terhadap Y.\n",
              format(uji$pvalue, digits = 4, scientific = TRUE), alpha))
} else {
  cat("**Keputusan:** H0 gagal ditolak; secara simultan X1–X6 tidak berpengaruh signifikan terhadap Y.\n")
}

Keputusan: H0 ditolak karena p-value = 1.157e-129 < α = 0.05.

Kesimpulan: secara simultan X1–X6 berpengaruh signifikan terhadap Y.

Pengujian Signifikansi Parameter Secara Parsial

\(H_0: \beta = 0\) (parameter tidak berpengaruh signifikan) lawan \(H_1: \beta \neq 0\). Statistik uji \(t = \hat\beta / SE(\hat\beta)\); \(H_0\) ditolak jika \(p\text{-value} < \alpha = 0{,}05\).

knitr::kable(uji$parsial, digits = 6)

Parameter	Estimasi	SE	t_hitung	p_value	Keputusan
satu	167.379165	77.340267	2.164192	0.030933	Tolak H0
X1	0.006721	0.130590	0.051469	0.958973	Gagal Tolak H0
X1_knot1	-0.884364	0.262858	-3.364418	0.000827	Tolak H0
X1_knot2	0.860350	0.897723	0.958368	0.338351	Gagal Tolak H0
X1_knot3	-2.205508	1.154933	-1.909642	0.056764	Gagal Tolak H0
X2	-0.162048	0.054276	-2.985645	0.002972	Tolak H0
X2_knot1	-0.027020	0.089657	-0.301371	0.763260	Gagal Tolak H0
X2_knot2	-0.313064	0.314288	-0.996104	0.319692	Gagal Tolak H0
X2_knot3	0.513304	0.416129	1.233522	0.217974	Gagal Tolak H0
X3	3.088515	2.182112	1.415379	0.157594	Gagal Tolak H0
X3_knot1	-1.193279	2.630559	-0.453622	0.650302	Gagal Tolak H0
X3_knot2	-6.518366	1.699148	-3.836257	0.000141	Tolak H0
X3_knot3	6.397219	1.432649	4.465308	0.000010	Tolak H0
X4	-0.572074	0.092514	-6.183673	0.000000	Tolak H0
X4_knot1	0.756445	0.243692	3.104102	0.002019	Tolak H0
X4_knot2	0.935417	1.166628	0.801812	0.423051	Gagal Tolak H0
X4_knot3	-2.197128	1.406362	-1.562278	0.118869	Gagal Tolak H0
X5	-1.828643	1.205513	-1.516900	0.129938	Gagal Tolak H0
X5_knot1	3.616134	1.352379	2.673907	0.007748	Tolak H0
X5_knot2	-2.564442	0.735499	-3.486668	0.000533	Tolak H0
X5_knot3	2.261039	0.709611	3.186310	0.001533	Tolak H0
X6	0.121548	0.153186	0.793471	0.427888	Gagal Tolak H0
X6_knot1	-0.235273	0.206243	-1.140756	0.254530	Gagal Tolak H0
X6_knot2	-0.185091	0.353273	-0.523931	0.600564	Gagal Tolak H0
X6_knot3	0.665479	0.484792	1.372710	0.170472	Gagal Tolak H0

sig <- uji$parsial[uji$parsial$p_value < alpha, ]
cat(sprintf("Dari %d parameter, terdapat %d parameter yang berpengaruh signifikan:\n\n",
            nrow(uji$parsial), nrow(sig)))

Dari 25 parameter, terdapat 10 parameter yang berpengaruh signifikan:

for (i in seq_len(nrow(sig))) cat("- ", sig$Parameter[i], " (p = ", format(sig$p_value[i], digits = 3), ")\n", sep = "")

satu (p = 0.0309)
X1_knot1 (p = 0.000827)
X2 (p = 0.00297)
X3_knot2 (p = 0.000141)
X3_knot3 (p = 9.94e-06)
X4 (p = 1.32e-09)
X4_knot1 (p = 0.00202)
X5_knot1 (p = 0.00775)
X5_knot2 (p = 0.000533)
X5_knot3 (p = 0.00153)

Interpretasi Model Regresi Nonparametrik Spline Truncated

Fungsi di bawah adalah kontribusi masing-masing variabel terhadap ŷ (variabel lain dan konstanta dianggap tetap), dipecah menurut interval knot.

fungsi_segmen <- function(B, knot_j, nama_var) {
  r <- length(knot_j)
  b_lin <- unname(B[nama_var])
  b_tr  <- unname(B[paste0(nama_var, "_knot", 1:r)])
  tab <- data.frame(Interval = character(r + 1), Kemiringan = numeric(r + 1),
                    Konstanta = numeric(r + 1), stringsAsFactors = FALSE)
  for (i in 0:r) {
    if (i == 0) { sl <- b_lin; ic <- 0 } else {
      sl <- b_lin + sum(b_tr[1:i]); ic <- -sum(b_tr[1:i] * knot_j[1:i])
    }
    tab$Interval[i + 1] <- if (i == 0) {
      paste0(nama_var, " < ", fmt(knot_j[1]))
    } else if (i == r) {
      paste0(nama_var, " >= ", fmt(knot_j[r]))
    } else {
      paste0(fmt(knot_j[i]), " <= ", nama_var, " < ", fmt(knot_j[i + 1]))
    }
    tab$Kemiringan[i + 1] <- sl
    tab$Konstanta[i + 1]  <- ic
  }
  tab
}

for (v in paste0("X", 1:6)) {
  cat("\n\n**Pengaruh", v, "terhadap Y**\n\n")
  print(knitr::kable(fungsi_segmen(uji$B, hasil$knot[v, ], v), digits = 5))
}

Pengaruh X1 terhadap Y

Interval	Kemiringan	Konstanta
X1 < 13,05286	0.00672	0.00000
13,05286 <= X1 < 24,60592	-0.87764	11.54348
24,60592 <= X1 < 28,45694	-0.01729	-9.62622
X1 >= 28,45694	-2.22280	53.13580

Pengaruh X2 terhadap Y

Interval	Kemiringan	Konstanta
X2 < 28,46286	-0.16205	0.00000
28,46286 <= X2 < 58,95878	-0.18907	0.76907
58,95878 <= X2 < 69,12408	-0.50213	19.22691
X2 >= 69,12408	0.01117	-16.25478

Pengaruh X3 terhadap Y

Interval	Kemiringan	Konstanta
X3 < 4,55429	3.08852	0.00000
4,55429 <= X3 < 8,06245	1.89524	5.43454
8,06245 <= X3 < 9,23184	-4.62313	57.98853
X3 >= 9,23184	1.77409	-1.06956

Pengaruh X4 terhadap Y

Interval	Kemiringan	Konstanta
X4 < 24,80571	-0.57207	0.00000
24,80571 <= X4 < 51,38327	0.18437	-18.76416
51,38327 <= X4 < 60,24245	1.11979	-66.82892
X4 >= 60,24245	-1.07734	65.53147

Pengaruh X5 terhadap Y

Interval	Kemiringan	Konstanta
X5 < 62,06571	-1.82864	0.00000
62,06571 <= X5 < 68,86469	1.78749	-224.43796
68,86469 <= X5 < 71,13102	-0.77695	-47.83844
X5 >= 71,13102	1.48409	-208.66848

Pengaruh X6 terhadap Y

Interval	Kemiringan	Konstanta
X6 < 14,34286	0.12155	0.00000
14,34286 <= X6 < 29,71020	-0.11372	3.37449
29,71020 <= X6 < 34,83265	-0.29882	8.87358
X6 >= 34,83265	0.36666	-14.30681

sig_var <- unique(sub("_knot[0-9]+$", "", sig$Parameter[sig$Parameter != "satu"]))
cat(sprintf("\nBerdasarkan GCV minimum (%s), model terbaik memiliki %d parameter, dengan %d parameter signifikan secara parsial. Variabel yang berpengaruh signifikan adalah %s.\n",
            fmt(hasil$mingcv), uji$k, nrow(sig[sig$Parameter != "satu", ]), paste(sig_var, collapse = ", ")))

Berdasarkan GCV minimum (57,87155), model terbaik memiliki 25 parameter, dengan 9 parameter signifikan secara parsial. Variabel yang berpengaruh signifikan adalah X1, X2, X3, X4, X5.

TUGAS 6_SPLINE TRUNCATED

yumnaa

2026-10-11