library(lmtest)
## Warning: package 'lmtest' was built under R version 4.6.1
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(MASS)
## Warning: package 'MASS' was built under R version 4.6.1
library(car)
## Loading required package: carData
library(pastecs)
## Warning: package 'pastecs' was built under R version 4.6.1
library(pracma)
## Warning: package 'pracma' was built under R version 4.6.1
##
## Attaching package: 'pracma'
## The following object is masked from 'package:car':
##
## logit
library(Matrix)
## Warning: package 'Matrix' was built under R version 4.6.1
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:pracma':
##
## expm, lu, tril, triu
data=read.table(file.choose(),header=TRUE)
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.77
summary(data)
## Y X1 X2 X3
## Min. : 2.240 Min. :60.05 Min. :11.15 Min. : 3.30
## 1st Qu.: 3.558 1st Qu.:65.42 1st Qu.:12.10 1st Qu.:10.25
## Median : 4.510 Median :69.14 Median :12.51 Median :18.64
## Mean : 4.957 Mean :68.99 Mean :12.65 Mean :26.36
## 3rd Qu.: 5.692 3rd Qu.:72.16 3rd Qu.:12.90 3rd Qu.:41.68
## Max. :12.360 Max. :77.73 Max. :14.94 Max. :75.04
## X4
## Min. :21.39
## 1st Qu.:50.24
## Median :67.97
## Mean :61.12
## 3rd Qu.:71.33
## Max. :80.77
stat.desc(data)
## Y X1 X2 X3 X4
## nbr.val 56.0000000 5.600000e+01 56.00000000 56.0000000 56.0000000
## nbr.null 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000
## nbr.na 0.0000000 0.000000e+00 0.00000000 0.0000000 0.0000000
## min 2.2400000 6.005000e+01 11.15000000 3.3000000 21.3900000
## max 12.3600000 7.773000e+01 14.94000000 75.0400000 80.7700000
## range 10.1200000 1.768000e+01 3.79000000 71.7400000 59.3800000
## sum 277.6000000 3.863250e+03 708.36000000 1476.2800000 3422.8700000
## median 4.5100000 6.914000e+01 12.50500000 18.6400000 67.9650000
## mean 4.9571429 6.898661e+01 12.64928571 26.3621429 61.1226786
## SE.mean 0.2715868 6.008087e-01 0.10790865 2.6750203 2.2021767
## CI.mean.0.95 0.5442721 1.204048e+00 0.21625376 5.3608606 4.4132607
## var 4.1305262 2.021438e+01 0.65207948 400.7210935 271.5766018
## std.dev 2.0323696 4.496041e+00 0.80751438 20.0180192 16.4795814
## coef.var 0.4099881 6.517266e-02 0.06383873 0.7593472 0.2696148
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
plot(data$X1,data$Y,main="Scatterplot Tingkat Pengangguran Terbuka dan
Tingkat Partisipasi Angkatan Kerja"
, xlab="Tingkat Partisipasi Angkatan Kerja (X1)"
, ylab="Tingkat Pengangguran Terbuka (Y)"
,col="black")
abline(lm(data$Y~data$X1), col="darkred", lwd=3)
plot(data$X2,data$Y,main="Scatterplot Tingkat Pengangguran Terbuka dan
Harapan Lama Sekolah"
, xlab="Harapan Lama Sekolah (X2)"
, ylab="Tingkat Pengangguran Terbuka (Y)"
,col="black")
abline(lm(data$Y~data$X2), col="darkred", lwd=3)
plot(data$X3,data$Y,main="Scatterplot Tingkat Pengangguran Terbuka dan
PDRB Atas Dasar Harga Berlaku"
, xlab="PDRB Atas Dasar Harga Berlaku (X3)"
, ylab="Tingkat Pengangguran Terbuka (Y)"
,col="black")
abline(lm(data$Y~data$X3), col="darkred", lwd=3)
plot(data$X4,data$Y,main="Scatterplot Tingkat Pengangguran Terbuka dan
Indeks Pembangunan Manusia"
, xlab="TIndeks Pembangunan Manusia (X4)"
, ylab="Tingkat Penggangguran Terbuka (Y)"
,col="black")
abline(lm(data$Y~data$X4), col="darkred", lwd=3)
# ============================================================
# BAGIAN 2: PEMILIHAN 1 TITIK KNOT (GCV)
# ============================================================
GCV1 = function(data, para = 0, nk = 50)
{
library(Matrix)
library(pracma)
data = as.matrix(data)
N = nrow(data)
M = ncol(data)
m = M - para - 1 # jumlah variabel prediktor
dataA = as.matrix(data[, (para + 2):M])
y = data[, 1]
X = data[, 2:M, drop = FALSE]
# kandidat knot: nk titik per variabel, buang titik min & max
knot1 = matrix(ncol = m, nrow = nk)
for (i in 1:m)
knot1[, i] = seq(min(dataA[, i]), max(dataA[, i]), length.out = nk)
knot1 = knot1[2:(nk - 1), , drop = FALSE]
colnames(knot1) = paste0("knot_x", 1:m)
nk1 = nrow(knot1)
aa = rep(1, N)
GCV = rep(NA, nk1)
MSE = rep(NA, nk1)
Rsq = rep(NA, nk1)
I_N = diag(N)
for (i in 1:nk1)
{
data1 = matrix(0, N, m)
for (j in 1:m) data1[, j] = pmax(dataA[, j] - knot1[i, j], 0)
mx = cbind(aa, X, data1)
C = pinv(t(mx) %*% mx)
B = C %*% (t(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
A = mx %*% C %*% t(mx)
A2 = (sum(diag(I_N - A)) / N)^2
GCV[i] = MSE[i] / A2
}
knotke = 1:nk1
dataAll = cbind(GCV = GCV, Rsq = Rsq, knot_ke = knotke, knot1)
write.csv(dataAll, file = "D:/DOKUMEN KULIAH/REGRESI NONPAR- SEM 5/TUGAS PERTEMUAN 5/Data knot 1.csv")
# urutkan berdasarkan GCV
dataG = dataAll[order(GCV), -2] # kolom: GCV, knot_ke, knot_x1..x4
cat("==============================================", "\n")
cat("HASIL GCV terkecil dengan 1 knot", "\n")
cat("==============================================", "\n")
print(dataG[1, ])
cat("\nNilai GCV 10 terkecil pertama", "\n")
print(dataG[1:10, ])
# ---- estimasi parameter di knot OPTIMAL ----
best = which.min(GCV)
knotopt = knot1[best, ]
datagcv1 = matrix(0, N, m)
for (j in 1:m) datagcv1[, j] = pmax(dataA[, j] - knotopt[j], 0)
mxgcv = cbind(aa, X, datagcv1)
Bopt = pinv(t(mxgcv) %*% mxgcv) %*% t(mxgcv) %*% y
rownames(Bopt) = c("b0",
paste0("b", 1:m, "_x", 1:m),
paste0("b", (m + 1):(2 * m), "_(x", 1:m, "-k", 1:m, ")+
"))
cat("\nKnot optimal ke-", best, " GCV =", min(GCV), "\n")
print(knotopt)
cat("\n==============================================", "\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 1 (OPTIMAL)", "\n")
cat("==============================================", "\n")
print(Bopt)
cat("\n")
invisible(list(knotgcv = knotopt, mingcv = min(GCV), B = Bopt, dataAll = dataAll))
}
hasil1 = GCV1(data)
## ==============================================
## HASIL GCV terkecil dengan 1 knot
## ==============================================
## GCV knot_ke knot_x1 knot_x2 knot_x3 knot_x4
## 1.508187 45.000000 76.286735 14.630612 69.183673 75.922653
##
## Nilai GCV 10 terkecil pertama
## GCV knot_ke knot_x1 knot_x2 knot_x3 knot_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
## knot_x1 knot_x2 knot_x3 knot_x4
## 76.28673 14.63061 69.18367 75.92265
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 1 (OPTIMAL)
## ==============================================
## [,1]
## b0 10.122393144
## b1_x1 -0.200697922
## b2_x2 0.762874047
## b3_x3 -0.005244856
## b4_x4 -0.017306261
## b5_(x1-k1)+\n 0.072850412
## b6_(x2-k2)+\n -21.843749933
## b7_(x3-k3)+\n 0.415138942
## b8_(x4-k4)+\n 1.053106627
# BAGIAN 3: UJI SIGNIFIKANSI UNTUK 1 KNOT (KNOT OPTIMAL)
library(pracma)
uji = function(alpha = 0.1, para = 0)
{
data = read.table("D:/DOKUMEN KULIAH/REGRESI NONPAR- SEM 5/DATA SKRIPSIw.txt",header = TRUE)
data = as.matrix(data)
# --- Baca file knot (nama baru) & ambil knot dengan GCV minimum ---
# File tersimpan di working directory R; cek dengan getwd()
knot_all = read.csv("D:/DOKUMEN KULIAH/REGRESI NONPAR- SEM 5/TUGAS PERTEMUAN 5/Data knot 1.csv", header = TRUE, row.names = 1)
best = which.min(knot_all$GCV)
# kolom 1 = GCV, 2 = Rsq, 3 = knot_ke, 4+ = titik knot
knot = as.matrix(knot_all[best, 4:ncol(knot_all)]) # 1 baris saja
cat("Knot optimal ke-", best, " GCV =", knot_all$GCV[best], "\n")
print(knot)
cat("\n")
y = data[, 1]
n = nrow(data)
m = para + 2
dataA = data[, m:(m + 3)]
satu = rep(1, n)
n1k = ncol(knot)
data.knot = matrix(0, n, n1k)
for (i in 1:n1k) data.knot[, i] = pmax(dataA[, i] - knot[1, i], 0)
# susunan: 1, x1, (x1-k1)+, x2, (x2-k2)+, x3, (x3-k3)+, x4, (x4-k4)+
mx = cbind(satu, data[, 2], data.knot[, 1], data[, 3], data.knot[, 2],
data[, 4], data.knot[, 3], data[, 5], data.knot[, 4])
B = pinv(t(mx) %*% mx) %*% t(mx) %*% y
p = nrow(B) # jumlah parameter (9)
yhat = mx %*% B
ybar = mean(y)
res = y - yhat
SSE = sum(res^2)
SSR = sum((yhat - ybar)^2)
SST = sum((y - ybar)^2)
MSE = SSE / (n - p) # PERBAIKAN: dibagi derajat bebas galat
MSR = SSR / (p - 1)
Rsq = SSR / (SSR + SSE) * 100
# --- Uji simultan (F) ---
Fhit = MSR / MSE
pvalue = pf(Fhit, p - 1, n - p, lower.tail = FALSE)
cat('---------------------------------------', '\n')
cat('Kesimpulan hasil uji simultan', '\n')
cat('---------------------------------------', '\n')
if (pvalue <= alpha) {
cat('Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan', '
\n\n')
} else {
cat('Gagal Tolak Ho yakni semua variabel bebas tidak berpengaruh signifikan', '\n\n')
}
# --- Uji parsial (t) ---
SE = sqrt(diag(MSE * pinv(t(mx) %*% mx)))
thit = rep(NA, p)
pval = rep(NA, p)
cat('---------------------------------------------', '\n')
cat('Kesimpulan hasil uji parsial', '\n')
cat('---------------------------------------------', '\n')
for (i in 1:p)
{
thit[i] = B[i, 1] / SE[i]
pval[i] = 2 * pt(abs(thit[i]), n - p, lower.tail = FALSE)
if (pval[i] <= alpha)
cat('Parameter', i - 1, ': Tolak Ho, signifikan, pvalue =', pval[i], '\
n') else
cat('Parameter', i - 1, ': Gagal tolak Ho, tidak signifikan, pvalue =
', pval[i], '\n')
}
cat('=============================================', '\n')
cat('Estimasi parameter, t hitung, p-value', '\n')
cat('=============================================', '\n')
print(cbind(Beta = B[, 1], SE = SE, t_hitung = thit, p_value = pval))
cat('\nAnalysis of Variance', '\n')
cat('=============================================', '\n')
cat('Sumber df SS MS Fhit', '\n')
cat('Regresi ', p - 1, ' ', SSR, ' ', MSR, ' ', Fhit, '\n')
cat('Error ', n - p, ' ', SSE, ' ', MSE, '\n')
cat('Total ', n - 1, ' ', SST, '\n')
cat('=============================================', '\n')
cat('s =', sqrt(MSE), ' Rsq =', Rsq, '\n')
cat('pvalue(F) =', pvalue, '\n')
write.csv(res, file = 'output_uji_residual_knot1.csv')
write.csv(mx, file = 'output_uji_mx_knot1.csv')
write.csv(yhat, file = 'output_uji_yhat_knot1.csv')
}
uji(0.1, 0)
## Knot optimal ke- 45 GCV = 1.508187
## knot_x1 knot_x2 knot_x3 knot_x4
## 45 76.28673 14.63061 69.18367 75.92265
##
## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan
##
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Parameter 0 : Tolak Ho, signifikan, pvalue = 0.06500353
## nParameter 1 : Tolak Ho, signifikan, pvalue = 2.835411e-05
## nParameter 2 : Gagal tolak Ho, tidak signifikan, pvalue =
## 0.9156743
## Parameter 3 : Tolak Ho, signifikan, pvalue = 0.01250812
## nParameter 4 : Tolak Ho, signifikan, pvalue = 1.020449e-05
## nParameter 5 : Gagal tolak Ho, tidak signifikan, pvalue =
## 0.541134
## Parameter 6 : Tolak Ho, signifikan, pvalue = 0.09644074
## nParameter 7 : Gagal tolak Ho, tidak signifikan, pvalue =
## 0.1203336
## Parameter 8 : Tolak Ho, signifikan, pvalue = 2.335269e-06
## n=============================================
## Estimasi parameter, t hitung, p-value
## =============================================
## Beta SE t_hitung p_value
## [1,] 10.122393145 5.35720770 1.8894905 6.500353e-02
## [2,] -0.200697922 0.04328026 -4.6371698 2.835411e-05
## [3,] 0.072850412 0.68432960 0.1064552 9.156743e-01
## [4,] 0.762874047 0.29372705 2.5972210 1.250812e-02
## [5,] -21.843749933 4.41930754 -4.9427992 1.020449e-05
## [6,] -0.005244856 0.00852003 -0.6155913 5.411340e-01
## [7,] 0.415138942 0.24472964 1.6963165 9.644074e-02
## [8,] -0.017306261 0.01093879 -1.5820999 1.203336e-01
## [9,] 1.053106627 0.19589216 5.3759508 2.335269e-06
##
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 8 167.6864 20.9608 16.55933
## Error 47 59.49258 1.2658
## Total 55 227.1789
## =============================================
## s = 1.125078 Rsq = 73.81246
## pvalue(F) = 2.451288e-11
r = read.csv("D:/DOKUMEN KULIAH/REGRESI NONPAR- SEM 5/TUGAS PERTEMUAN 5/Data knot 1.csv", row.names = 1)
k48 = as.numeric(r[48, 4:7]) # knot kandidat terakhir
X = as.matrix(data[, 2:5]); y = data[, 1]
mx = cbind(1, X, sapply(1:4, function(j) pmax(X[, j] - k48[j], 0)))
print(pinv(t(mx) %*% mx) %*% t(mx) %*% y)
## [,1]
## [1,] 1.100168e+01
## [2,] -2.352118e-01
## [3,] 8.683805e-01
## [4,] -5.401921e-03
## [5,] -1.233868e-02
## [6,] 5.699848e-01
## [7,] -1.336717e+02
## [8,] 4.073604e-01
## [9,] 6.923471e+00
library(pracma)
data = read.table("D:/DOKUMEN KULIAH/REGRESI NONPAR- SEM 5/DATA SKRIPSIw.txt",header = TRUE)
nrow(data)
## [1] 56
buat_mx = function(X, knot, urut = "gcv")
{
N = nrow(X); m = ncol(X); k = nrow(knot)
if (urut == "gcv") {
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 {
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
}
GCVk = function(data, k = 2, nk = 50)
{
data = as.matrix(data)
N = nrow(data)
y = data[, 1]
X = data[, 2:ncol(data), drop = FALSE]
m = ncol(X)
# grid kandidat knot per variabel (buang titik min & max)
grid = sapply(1:m, function(j) seq(min(X[, j]), max(X[, j]), length.out = nk))
grid = grid[2:(nk - 1), , drop = FALSE]
ng = nrow(grid)
komb = t(combn(ng, k)) # tiap baris = indeks kandidat (naik)
nkomb = nrow(komb)
cat("Jumlah kombinasi knot:", nkomb, "\n")
GCV = rep(NA, nkomb)
Rsq = rep(NA, nkomb)
knotmat = matrix(NA, nkomb, k * m)
nama.knot = as.vector(t(outer(1:k, 1:m, function(s, j) paste0("k", s, "_x",j))))
for (i in 1:nkomb)
{
knot = matrix(0, k, m)
for (s in 1:k) knot[s, ] = grid[komb[i, s], ]
knotmat[i, ] = as.vector(t(knot))
mx = buat_mx(X, knot, "gcv")
C = pinv(t(mx) %*% mx)
B = C %*% (t(mx) %*% y)
yhat = mx %*% B
SSE = sum((y - yhat)^2)
SSR = sum((yhat - mean(y))^2)
trA = sum((mx %*% C) * mx) # trace matriks hat
GCV[i] = (SSE / N) / (((N - trA) / N)^2)
Rsq[i] = SSR / (SSR + SSE) * 100
if (i %% 2000 == 0) cat(" selesai", i, "dari", nkomb, "\n")
}
dataAll = cbind(GCV = GCV, Rsq = Rsq, komb_ke = 1:nkomb, knotmat)
colnames(dataAll)[4:ncol(dataAll)] = nama.knot
file.out = paste0("dataknot",k,"_v2.csv")
write.csv(dataAll, file = file.out)
dataG = dataAll[order(GCV), -2] # GCV, komb_ke, knot...
cat("\n==============================================\n")
cat("HASIL GCV terkecil dengan", k, "knot\n")
cat("==============================================\n")
print(dataG[1, ])
cat("\nNilai GCV 10 terkecil pertama\n")
print(dataG[1:10, ])
# ---- estimasi parameter pada knot OPTIMAL ----
best = which.min(GCV)
knotopt = matrix(knotmat[best, ], nrow = k, byrow = TRUE,
dimnames = list(paste0("knot_", 1:k), paste0("x", 1:m)))
mxopt = buat_mx(X, knotopt, "gcv")
Bopt = pinv(t(mxopt) %*% mxopt) %*% t(mxopt) %*% y
rownames(Bopt) = colnames(mxopt)
cat("\nKombinasi knot optimal ke-", best, " GCV =", min(GCV), "\n")
print(knotopt)
cat("\n==============================================\n")
cat("HASIL ESTIMASI PARAMETER", k, "TITIK KNOT (OPTIMAL)\n")
cat("==============================================\n")
print(Bopt)
cat("\nFile knot tersimpan di:", file.out, "\n")
invisible(list(knotopt = knotopt, mingcv = min(GCV), B = Bopt, dataAll = dataAll))
}
# ============================================================
# BAGIAN C: UJI SIGNIFIKANSI SIMULTAN & PARSIAL (k = 2 atau 3)
# ============================================================
uji_k = function(data, k = 2, alpha = 0.1)
{
data = as.matrix(data)
y = data[, 1]
X = data[, 2:ncol(data), drop = FALSE]
n = nrow(data); m = ncol(X)
# ambil knot dengan GCV minimum dari file hasil GCVk
file.knot = paste0("dataknot",k,"_v2.csv")
knot_all = read.csv(file.knot, header = TRUE, row.names = 1)
best = which.min(knot_all$GCV)
knot = matrix(as.numeric(knot_all[best, 4:ncol(knot_all)]), nrow = k, byrow= TRUE,
dimnames = list(paste0("knot_", 1:k), paste0("x", 1:m)))
cat("Kombinasi knot optimal ke-", best, " GCV =", knot_all$GCV[best], "\n")
print(knot)
cat("\n")
mx = buat_mx(X, knot, "uji")
B = pinv(t(mx) %*% mx) %*% t(mx) %*% y
p = nrow(B)
yhat = mx %*% B
ybar = mean(y)
res = y - yhat
SSE = sum(res^2)
SSR = sum((yhat - ybar)^2)
SST = sum((y - ybar)^2)
MSE = SSE / (n - p)
MSR = SSR / (p - 1)
Rsq = SSR / (SSR + SSE) * 100
# --- Uji simultan ---
Fhit = MSR / MSE
pvalue = pf(Fhit, p - 1, n - p, lower.tail = FALSE)
cat('---------------------------------------\n')
cat('Kesimpulan hasil uji simultan\n')
cat('---------------------------------------\n')
if (pvalue <= alpha) {
cat('Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan\n\n')
} else {
cat('Gagal Tolak Ho yakni semua variabel bebas tidak berpengaruh signifikan\n\n')
}
# --- Uji parsial ---
SE = sqrt(diag(MSE * pinv(t(mx) %*% mx)))
thit = B[, 1] / SE
pval = 2 * pt(abs(thit), n - p, lower.tail = FALSE)
cat('---------------------------------------------\n')
cat('Kesimpulan hasil uji parsial\n')
cat('---------------------------------------------\n')
for (i in 1:p)
{
if (pval[i] <= alpha)
cat('Parameter', i - 1, colnames(mx)[i], ': Tolak Ho, signifikan, pvalue=', pval[i], '\n') else
cat('Parameter', i - 1, colnames(mx)[i], ': Gagal tolak Ho, tidak signifikan, pvalue =', pval[i], '\n')
}
cat('=============================================\n')
cat('Estimasi parameter, t hitung, p-value\n')
cat('=============================================\n')
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('=============================================\n')
cat('Sumber df SS MS Fhit\n')
cat('Regresi ', p - 1, ' ', SSR, ' ', MSR, ' ', Fhit, '\n')
cat('Error ', n - p, ' ', SSE, ' ', MSE, '\n')
cat('Total ', n - 1, ' ', SST, '\n')
cat('=============================================\n')
cat('s =', sqrt(MSE), ' Rsq =', Rsq, '\n')
cat('pvalue(F) =', pvalue, '\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))
}
hasil2 = GCVk(data, k = 2)
## Jumlah kombinasi knot: 1128
##
## ==============================================
## HASIL GCV terkecil dengan 2 knot
## ==============================================
## GCV komb_ke k1_x1 k1_x2 k1_x3 k1_x4 k2_x1
## 1.316068 135.000000 61.132449 11.382041 7.692245 25.025510 76.286735
## k2_x2 k2_x3 k2_x4
## 14.630612 69.183673 75.922653
##
## Nilai GCV 10 terkecil pertama
## GCV komb_ke k1_x1 k1_x2 k1_x3 k1_x4 k2_x1 k2_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
## k2_x3 k2_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
##
## Kombinasi 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
##
## ==============================================
## HASIL ESTIMASI PARAMETER 2 TITIK KNOT (OPTIMAL)
## ==============================================
## [,1]
## b0 95.3587174
## x1 -2.9213119
## x2 6.1469387
## x3 0.2059982
## x4 0.7283636
## (x1-k1)+\n 2.7449916
## (x2-k1)+\n -5.3197766
## (x3-k1)+\n -0.2144538
## (x4-k1)+\n -0.7619877
## (x1-k2)+\n 0.7954129
## (x2-k2)+\n -18.0705567
## (x3-k2)+\n 0.2853691
## (x4-k2)+\n 0.9034632
##
## File knot tersimpan di: dataknot2_v2.csv
uji_k(data,k=2,alpha=0.1)
## Kombinasi 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
##
## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Parameter 0 b0 : Gagal tolak Ho, tidak signifikan, pvalue = 0.2671925
## Parameter 1 x1 : Tolak Ho, signifikan, pvalue= 0.007544198
## Parameter 2 (x1-k1)+ : Tolak Ho, signifikan, pvalue= 0.01250391
## Parameter 3 (x1-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3045796
## Parameter 4 x2 : Gagal tolak Ho, tidak signifikan, pvalue = 0.2127214
## Parameter 5 (x2-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2845971
## Parameter 6 (x2-k2)+ : Tolak Ho, signifikan, pvalue= 0.0001058266
## Parameter 7 x3 : Gagal tolak Ho, tidak signifikan, pvalue = 0.2122837
## Parameter 8 (x3-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2057346
## Parameter 9 (x3-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2224808
## Parameter 10 x4 : Tolak Ho, signifikan, pvalue= 0.05283963
## Parameter 11 (x4-k1)+ : Tolak Ho, signifikan, pvalue= 0.0465475
## Parameter 12 (x4-k2)+ : Tolak Ho, signifikan, pvalue= 2.97441e-05
## =============================================
## Estimasi parameter, t hitung, p-value
## =============================================
## Beta SE t_hitung p_value
## b0 95.3587187 84.8285640 1.124135 0.2671924526
## x1 -2.9213119 1.0418249 -2.804034 0.0075441983
## (x1-k1)+ 2.7449916 1.0529268 2.607011 0.0125039142
## (x1-k2)+ 0.7954129 0.7655043 1.039070 0.3045796353
## x2 6.1469386 4.8596443 1.264895 0.2127214069
## (x2-k1)+ -5.3197765 4.9095217 -1.083563 0.2845971114
## (x2-k2)+ -18.0705567 4.2319010 -4.270080 0.0001058266
## x3 0.2059982 0.1626993 1.266128 0.2122837025
## (x3-k1)+ -0.2144538 0.1669140 -1.284816 0.2057346473
## (x3-k2)+ 0.2853691 0.2305324 1.237870 0.2224808174
## x4 0.7283636 0.3657962 1.991173 0.0528396324
## (x4-k1)+ -0.7619877 0.3717972 -2.049471 0.0465474968
## (x4-k2)+ 0.9034632 0.1935247 4.668464 0.0000297441
##
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 12 183.7252 15.31043 15.15057
## Error 43 43.45373 1.010552
## Total 55 227.1789
## =============================================
## s = 1.005262 Rsq = 80.87246
## pvalue(F) = 9.568368e-12
#PEMILIHAN TIGA TITIK KNOT#
hasil3=GCVk(data, k=3)
## Jumlah kombinasi knot: 17296
## selesai 2000 dari 17296
## selesai 4000 dari 17296
## selesai 6000 dari 17296
## selesai 8000 dari 17296
## selesai 10000 dari 17296
## selesai 12000 dari 17296
## selesai 14000 dari 17296
## selesai 16000 dari 17296
##
## ==============================================
## HASIL GCV terkecil dengan 3 knot
## ==============================================
## GCV komb_ke k1_x1 k1_x2 k1_x3 k1_x4
## 1.444246 2200.000000 61.132449 11.382041 7.692245 25.025510
## k2_x1 k2_x2 k2_x3 k2_x4 k3_x1 k3_x2
## 61.854082 11.536735 10.620408 27.449184 76.286735 14.630612
## k3_x3 k3_x4
## 69.183673 75.922653
##
## Nilai GCV 10 terkecil pertama
## GCV komb_ke k1_x1 k1_x2 k1_x3 k1_x4 k2_x1 k2_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
## k2_x3 k2_x4 k3_x1 k3_x2 k3_x3 k3_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
##
## Kombinasi 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
##
## ==============================================
## HASIL ESTIMASI PARAMETER 3 TITIK KNOT (OPTIMAL)
## ==============================================
## [,1]
## b0 -4.87954771
## x1 -0.57070096
## x2 1.29861828
## x3 0.20162854
## x4 1.23897093
## (x1-k1)+\n -2.60712938
## (x2-k1)+\n 13.80910091
## (x3-k1)+\n -0.22960609
## (x4-k1)+\n -1.94377923
## (x1-k2)+\n 3.03208076
## (x2-k2)+\n -14.12188191
## (x3-k2)+\n 0.02361872
## (x4-k2)+\n 0.67438515
## (x1-k3)+\n 0.76491110
## (x2-k3)+\n -19.59971920
## (x3-k3)+\n 0.25980650
## (x4-k3)+\n 0.92445503
##
## File knot tersimpan di: dataknot3_v2.csv
uji_k(data,k=3,alpha=0.1)
## Kombinasi 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
##
## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Parameter 0 b0 : Gagal tolak Ho, tidak signifikan, pvalue = 0.6691134
## Parameter 1 x1 : Gagal tolak Ho, tidak signifikan, pvalue = 0.8101606
## Parameter 2 (x1-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.5935914
## Parameter 3 (x1-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2569091
## Parameter 4 (x1-k3)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3875066
## Parameter 5 x2 : Gagal tolak Ho, tidak signifikan, pvalue = 0.9145372
## Parameter 6 (x2-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.6780688
## Parameter 7 (x2-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.5260547
## Parameter 8 (x2-k3)+ : Tolak Ho, signifikan, pvalue= 0.0001404269
## Parameter 9 x3 : Gagal tolak Ho, tidak signifikan, pvalue = 0.5340848
## Parameter 10 (x3-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.7141141
## Parameter 11 (x3-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.942884
## Parameter 12 (x3-k3)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2870221
## Parameter 13 x4 : Tolak Ho, signifikan, pvalue= 0.09910266
## Parameter 14 (x4-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.1848397
## Parameter 15 (x4-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3857466
## Parameter 16 (x4-k3)+ : Tolak Ho, signifikan, pvalue= 4.564589e-05
## =============================================
## Estimasi parameter, t hitung, p-value
## =============================================
## Beta SE t_hitung p_value
## b0 -4.87954766 11.3313797 -0.43062255 6.691134e-01
## x1 -0.57070097 2.3596928 -0.24185393 8.101606e-01
## (x1-k1)+ -2.60712937 4.8454015 -0.53806261 5.935914e-01
## (x1-k2)+ 3.03208075 2.6352515 1.15058497 2.569091e-01
## (x1-k3)+ 0.76491110 0.8752551 0.87392932 3.875066e-01
## x2 1.29861831 12.0225195 0.10801549 9.145372e-01
## (x2-k1)+ 13.80910082 33.0174967 0.41823585 6.780688e-01
## (x2-k2)+ -14.12188185 22.0729359 -0.63978267 5.260547e-01
## (x2-k3)+ -19.59971920 4.6428330 -4.22149996 1.404269e-04
## x3 0.20162854 0.3213964 0.62735156 5.340848e-01
## (x3-k1)+ -0.22960609 0.6222171 -0.36901280 7.141141e-01
## (x3-k2)+ 0.02361872 0.3275433 0.07210868 9.428840e-01
## (x3-k3)+ 0.25980650 0.2406862 1.07944073 2.870221e-01
## x4 1.23897093 0.7333330 1.68950652 9.910266e-02
## (x4-k1)+ -1.94377922 1.4399711 -1.34987380 1.848397e-01
## (x4-k2)+ 0.67438515 0.7687876 0.87720604 3.857466e-01
## (x4-k3)+ 0.92445502 0.2015338 4.58709688 4.564589e-05
##
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 16 185.9148 11.61967 10.9821
## Error 39 41.26416 1.058056
## Total 55 227.1789
## =============================================
## s = 1.028618 Rsq = 81.83627
## pvalue(F) = 7.291165e-10
r = read.csv("dataknot2_v2.csv", row.names = 1)
k.last = matrix(as.numeric(r[nrow(r), 4:ncol(r)]), nrow = 2, byrow = TRUE)
X = as.matrix(data[, 2:5]); y = data[, 1]
mx = buat_mx(X, k.last, "gcv")
print(pinv(t(mx) %*% mx) %*% t(mx) %*% y)
## [,1]
## [1,] 6.044955e+00
## [2,] -1.859007e-01
## [3,] 1.031486e+00
## [4,] -4.170952e-03
## [5,] -2.228411e-02
## [6,] 9.154406e+00
## [7,] -3.547727e+02
## [8,] 4.083286e-02
## [9,] 4.845712e+00
## [10,] -2.164657e+01
## [11,] 6.943980e+02
## [12,] 2.041643e-02
## [13,] -1.020433e+01
getwd()
## [1] "D:/DOKUMEN KULIAH/REGRESI NONPAR- SEM 5"