library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(MASS)
library(car)
## Loading required package: carData
library(pastecs)
library(pracma)
## 
## Attaching package: 'pracma'
## The following object is masked from 'package:car':
## 
##     logit
library(Matrix)
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:pracma':
## 
##     expm, lu, tril, triu
library(ggplot2)
library(gridExtra)
data=read.table(file.choose(), header=TRUE)  # ganti dengan nama file saat knit
data
stat.desc(data)
model=(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.498612    -0.237422     0.919839    -0.008344    -0.009455
summary(model)
## 
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1296 -0.9922 -0.1119  0.7639  4.6374 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.498612   5.796395   1.811  0.07600 .  
## X1          -0.237422   0.049782  -4.769 1.59e-05 ***
## X2           0.919839   0.278410   3.304  0.00175 ** 
## X3          -0.008344   0.010186  -0.819  0.41651    
## X4          -0.009455   0.012466  -0.758  0.45170    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.453 on 51 degrees of freedom
## Multiple R-squared:  0.526,  Adjusted R-squared:  0.4888 
## F-statistic: 14.15 on 4 and 51 DF,  p-value: 7.786e-08
vif(model)
##       X1       X2       X3       X4 
## 1.304895 1.316581 1.083064 1.099406
par(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.508187      45 76.28673 14.63061 69.183673 75.92265
##  [2,] 1.515931      44 75.92592 14.55327 67.719592 74.71082
##  [3,] 1.603414      46 76.64755 14.70796 70.647755 77.13449
##  [4,] 1.615275      43 75.56510 14.47592 66.255510 73.49898
##  [5,] 1.734461      42 75.20429 14.39857 64.791429 72.28714
##  [6,] 1.763515      47 77.00837 14.78531 72.111837 78.34633
##  [7,] 1.830095       3 61.13245 11.38204  7.692245 25.02551
##  [8,] 1.861406       4 61.49327 11.45939  9.156327 26.23735
##  [9,] 1.869770      41 74.84347 14.32122 63.327347 71.07531
## [10,] 1.912402       5 61.85408 11.53673 10.620408 27.44918
## 
## Knot optimal ke-45  GCV = 1.508187
##              x1       x2       x3       x4
## knot_1 76.28673 14.63061 69.18367 75.92265
## 
## Estimasi parameter
##               Estimasi
## b0        10.122393144
## x1        -0.200697922
## x2         0.762874047
## x3        -0.005244856
## x4        -0.017306261
## (x1-k1)+   0.072850412
## (x2-k1)+ -21.843749933
## (x3-k1)+   0.415138942
## (x4-k1)+   1.053106627
hasil2 <- GCVk(data, k = 2)
## 
## === 2 titik knot: 10 GCV terkecil ===
##            GCV knot_ke knot1_x1 knot1_x2  knot1_x3 knot1_x4 knot2_x1 knot2_x2
##  [1,] 1.316068     135 61.13245 11.38204  7.692245 25.02551 76.28673 14.63061
##  [2,] 1.317894     179 61.49327 11.45939  9.156327 26.23735 76.28673 14.63061
##  [3,] 1.321114     178 61.49327 11.45939  9.156327 26.23735 75.92592 14.55327
##  [4,] 1.321402     134 61.13245 11.38204  7.692245 25.02551 75.92592 14.55327
##  [5,] 1.358738     221 61.85408 11.53673 10.620408 27.44918 75.92592 14.55327
##  [6,] 1.359630     222 61.85408 11.53673 10.620408 27.44918 76.28673 14.63061
##  [7,] 1.359989     138 61.13245 11.38204  7.692245 25.02551 77.36918 14.86265
##  [8,] 1.388799     136 61.13245 11.38204  7.692245 25.02551 76.64755 14.70796
##  [9,] 1.393180     180 61.49327 11.45939  9.156327 26.23735 76.64755 14.70796
## [10,] 1.396087     177 61.49327 11.45939  9.156327 26.23735 75.56510 14.47592
##       knot2_x3 knot2_x4
##  [1,] 69.18367 75.92265
##  [2,] 69.18367 75.92265
##  [3,] 67.71959 74.71082
##  [4,] 67.71959 74.71082
##  [5,] 67.71959 74.71082
##  [6,] 69.18367 75.92265
##  [7,] 73.57592 79.55816
##  [8,] 70.64776 77.13449
##  [9,] 70.64776 77.13449
## [10,] 66.25551 73.49898
## 
## Knot optimal ke-135  GCV = 1.316068
##              x1       x2        x3       x4
## knot_1 61.13245 11.38204  7.692245 25.02551
## knot_2 76.28673 14.63061 69.183673 75.92265
## 
## Estimasi parameter
##             Estimasi
## b0        95.3587174
## x1        -2.9213119
## x2         6.1469387
## x3         0.2059982
## x4         0.7283636
## (x1-k1)+   2.7449916
## (x2-k1)+  -5.3197766
## (x3-k1)+  -0.2144538
## (x4-k1)+  -0.7619877
## (x1-k2)+   0.7954129
## (x2-k2)+ -18.0705567
## (x3-k2)+   0.2853691
## (x4-k2)+   0.9034632
hasil3 <- GCVk(data, k = 3)
## 
## === 3 titik knot: 10 GCV terkecil ===
##            GCV knot_ke knot1_x1 knot1_x2 knot1_x3 knot1_x4 knot2_x1 knot2_x2
##  [1,] 1.444246    2200 61.13245 11.38204 7.692245 25.02551 61.85408 11.53673
##  [2,] 1.445793    3146 61.49327 11.45939 9.156327 26.23735 61.85408 11.53673
##  [3,] 1.447104    2157 61.13245 11.38204 7.692245 25.02551 61.49327 11.45939
##  [4,] 1.452451    3106 61.13245 11.38204 7.692245 25.02551 77.00837 14.78531
##  [5,] 1.453903    2199 61.13245 11.38204 7.692245 25.02551 61.85408 11.53673
##  [6,] 1.454858    3145 61.49327 11.45939 9.156327 26.23735 61.85408 11.53673
##  [7,] 1.458180    2156 61.13245 11.38204 7.692245 25.02551 61.49327 11.45939
##  [8,] 1.458931    2242 61.13245 11.38204 7.692245 25.02551 62.21490 11.61408
##  [9,] 1.464584    3037 61.13245 11.38204 7.692245 25.02551 73.03939 13.93449
## [10,] 1.465580    3101 61.13245 11.38204 7.692245 25.02551 76.28673 14.63061
##        knot2_x3 knot2_x4 knot3_x1 knot3_x2 knot3_x3 knot3_x4
##  [1,] 10.620408 27.44918 76.28673 14.63061 69.18367 75.92265
##  [2,] 10.620408 27.44918 76.28673 14.63061 69.18367 75.92265
##  [3,]  9.156327 26.23735 76.28673 14.63061 69.18367 75.92265
##  [4,] 72.111837 78.34633 77.36918 14.86265 73.57592 79.55816
##  [5,] 10.620408 27.44918 75.92592 14.55327 67.71959 74.71082
##  [6,] 10.620408 27.44918 75.92592 14.55327 67.71959 74.71082
##  [7,]  9.156327 26.23735 75.92592 14.55327 67.71959 74.71082
##  [8,] 12.084490 28.66102 76.28673 14.63061 69.18367 75.92265
##  [9,] 56.006939 65.01612 76.28673 14.63061 69.18367 75.92265
## [10,] 69.183673 75.92265 76.64755 14.70796 70.64776 77.13449
## 
## Knot optimal ke-2200  GCV = 1.444246
##              x1       x2        x3       x4
## knot_1 61.13245 11.38204  7.692245 25.02551
## knot_2 61.85408 11.53673 10.620408 27.44918
## knot_3 76.28673 14.63061 69.183673 75.92265
## 
## Estimasi parameter
##              Estimasi
## b0        -4.87954771
## x1        -0.57070096
## x2         1.29861828
## x3         0.20162854
## x4         1.23897093
## (x1-k1)+  -2.60712938
## (x2-k1)+  13.80910091
## (x3-k1)+  -0.22960609
## (x4-k1)+  -1.94377923
## (x1-k2)+   3.03208076
## (x2-k2)+ -14.12188191
## (x3-k2)+   0.02361872
## (x4-k2)+   0.67438515
## (x1-k3)+   0.76491110
## (x2-k3)+ -19.59971920
## (x3-k3)+   0.25980650
## (x4-k3)+   0.92445503

Perbandingan GCV

gcv_perbandingan <-
  data.frame(
    Model = c("1 Knot", "2 Knot", "3 Knot"),
    GCV   = c(hasil1$mingcv, hasil2$mingcv, hasil3$mingcv)
  )

gcv_perbandingan
ggplot(
  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 Knot
cat("Nilai GCV:", gcv_minimum, "\n")
## Nilai GCV: 1.316068
knot_optimal <- hasil2$knotopt
knot_optimal
##              x1       x2        x3       x4
## knot_1 61.13245 11.38204  7.692245 25.02551
## knot_2 76.28673 14.63061 69.183673 75.92265
B_optimal <- hasil2$B
B_optimal
##             Estimasi
## b0        95.3587174
## x1        -2.9213119
## x2         6.1469387
## x3         0.2059982
## x4         0.7283636
## (x1-k1)+   2.7449916
## (x2-k1)+  -5.3197766
## (x3-k1)+  -0.2144538
## (x4-k1)+  -0.7619877
## (x1-k2)+   0.7954129
## (x2-k2)+ -18.0705567
## (x3-k2)+   0.2853691
## (x4-k2)+   0.9034632

Uji signifikansi parameter

uji2 <- uji_k(data, hasil2, alpha = 0.1)
##              x1       x2        x3       x4
## knot_1 61.13245 11.38204  7.692245 25.02551
## knot_2 76.28673 14.63061 69.183673 75.92265
##                 Beta         SE  t_hitung      p_value
## b0        95.3587174 84.8285642  1.124135 0.2671924603
## x1        -2.9213119  1.0418249 -2.804034 0.0075441985
## x2         6.1469387  4.8596443  1.264895 0.2127214042
## x3         0.2059982  0.1626993  1.266128 0.2122837024
## x4         0.7283636  0.3657962  1.991173 0.0528396322
## (x1-k1)+   2.7449916  1.0529268  2.607011 0.0125039146
## (x2-k1)+  -5.3197766  4.9095217 -1.083563 0.2845971079
## (x3-k1)+  -0.2144538  0.1669140 -1.284816 0.2057346472
## (x4-k1)+  -0.7619877  0.3717972 -2.049471 0.0465474966
## (x1-k2)+   0.7954129  0.7655043  1.039070 0.3045796347
## (x2-k2)+ -18.0705567  4.2319010 -4.270080 0.0001058266
## (x3-k2)+   0.2853691  0.2305324  1.237870 0.2224808174
## (x4-k2)+   0.9034632  0.1935247  4.668464 0.0000297441
## 
## Analysis of Variance
## Regresi: 12 183.7252 15.31043 15.15057 
## Error  : 43 43.45373 1.010552 
## Total  : 55 227.1789 
## s = 1.005262 
## R-squared = 80.87246 
## p-value F = 9.568368e-12 
## 
## Uji simultan:
## Tolak H0: minimal terdapat satu parameter yang signifikan.
## 
## Uji parsial:
## b0 : tidak signifikan, p-value = 0.2671925 
## x1 : signifikan, p-value = 0.007544199 
## x2 : tidak signifikan, p-value = 0.2127214 
## x3 : tidak signifikan, p-value = 0.2122837 
## x4 : signifikan, p-value = 0.05283963 
## (x1-k1)+ : signifikan, p-value = 0.01250391 
## (x2-k1)+ : tidak signifikan, p-value = 0.2845971 
## (x3-k1)+ : tidak signifikan, p-value = 0.2057346 
## (x4-k1)+ : signifikan, p-value = 0.0465475 
## (x1-k2)+ : tidak signifikan, p-value = 0.3045796 
## (x2-k2)+ : signifikan, p-value = 0.0001058266 
## (x3-k2)+ : tidak signifikan, p-value = 0.2224808 
## (x4-k2)+ : signifikan, p-value = 2.97441e-05
uji3 <- uji_k(data, hasil3, alpha = 0.1)
##              x1       x2        x3       x4
## knot_1 61.13245 11.38204  7.692245 25.02551
## knot_2 61.85408 11.53673 10.620408 27.44918
## knot_3 76.28673 14.63061 69.183673 75.92265
##                  Beta         SE    t_hitung      p_value
## b0        -4.87954771 11.3313795 -0.43062257 6.691134e-01
## x1        -0.57070096  2.3596928 -0.24185393 8.101606e-01
## x2         1.29861828 12.0225196  0.10801548 9.145372e-01
## x3         0.20162854  0.3213964  0.62735155 5.340848e-01
## x4         1.23897093  0.7333330  1.68950652 9.910266e-02
## (x1-k1)+  -2.60712938  4.8454015 -0.53806261 5.935914e-01
## (x2-k1)+  13.80910091 33.0174969  0.41823585 6.780688e-01
## (x3-k1)+  -0.22960609  0.6222171 -0.36901280 7.141141e-01
## (x4-k1)+  -1.94377923  1.4399711 -1.34987380 1.848397e-01
## (x1-k2)+   3.03208076  2.6352515  1.15058497 2.569091e-01
## (x2-k2)+ -14.12188191 22.0729360 -0.63978267 5.260547e-01
## (x3-k2)+   0.02361872  0.3275433  0.07210868 9.428840e-01
## (x4-k2)+   0.67438515  0.7687876  0.87720605 3.857466e-01
## (x1-k3)+   0.76491110  0.8752551  0.87392932 3.875066e-01
## (x2-k3)+ -19.59971920  4.6428330 -4.22149996 1.404269e-04
## (x3-k3)+   0.25980650  0.2406862  1.07944073 2.870221e-01
## (x4-k3)+   0.92445503  0.2015338  4.58709688 4.564589e-05
## 
## Analysis of Variance
## Regresi: 16 185.9148 11.61967 10.9821 
## Error  : 39 41.26416 1.058056 
## Total  : 55 227.1789 
## s = 1.028618 
## R-squared = 81.83627 
## p-value F = 7.291165e-10 
## 
## Uji simultan:
## Tolak H0: minimal terdapat satu parameter yang signifikan.
## 
## Uji parsial:
## b0 : tidak signifikan, p-value = 0.6691134 
## x1 : tidak signifikan, p-value = 0.8101606 
## x2 : tidak signifikan, p-value = 0.9145372 
## x3 : tidak signifikan, p-value = 0.5340848 
## x4 : signifikan, p-value = 0.09910266 
## (x1-k1)+ : tidak signifikan, p-value = 0.5935914 
## (x2-k1)+ : tidak signifikan, p-value = 0.6780688 
## (x3-k1)+ : tidak signifikan, p-value = 0.7141141 
## (x4-k1)+ : tidak signifikan, p-value = 0.1848397 
## (x1-k2)+ : tidak signifikan, p-value = 0.2569091 
## (x2-k2)+ : tidak signifikan, p-value = 0.5260547 
## (x3-k2)+ : tidak signifikan, p-value = 0.942884 
## (x4-k2)+ : tidak signifikan, p-value = 0.3857466 
## (x1-k3)+ : tidak signifikan, p-value = 0.3875066 
## (x2-k3)+ : signifikan, p-value = 0.0001404269 
## (x3-k3)+ : tidak signifikan, p-value = 0.2870221 
## (x4-k3)+ : signifikan, p-value = 4.564589e-05

Plot 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()