library(readxl)
data1=read_excel("D:/Artikel spline TPT/spline tpt jabar.xlsx")
data=as.data.frame(data1)
data1## # A tibble: 27 × 6
## Y x1 x2 x3 x4 x5
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 8.47 64.2 73.0 89.5 43.9 5.19
## 2 7.32 67.8 69.7 87.9 47.1 5.17
## 3 7.71 72.3 68.2 106. 49.1 5.16
## 4 6.52 67.1 74.0 85.3 45.0 4.97
## 5 7.33 70.1 69.2 107. 49.0 4.94
## 6 3.89 68.4 69.4 89.7 49.0 4.69
## 7 3.52 66.3 73.1 83.9 47.1 4.99
## 8 9.49 62.0 71.0 107. 47.3 5.25
## 9 7.65 66.2 71.8 90.0 43.3 4.85
## 10 4.12 68.5 70.8 97.4 44.8 6.15
## # ℹ 17 more rowsData yang digunakan dalam penelitian ini merupakan data sekunder yang didapatkan dari Website Badan Pusat Statistik (BPS) Jawa Barat tahun 2023.
Y :Tingkat Pengangguran Terbuka
X1: Tingkat Partisipasi Angakatan Kerja
X2: Indeks Pembangunan Manusia
X3: APK Sekolah Menengah
X4: Dependency Ratio
X5: Laju Pertumbuhan PDRB
fungsi membuat matriks x truncated untuk x=1 variabel.k dibuat menjadi vektor yang elemennya adalah knot knot yang akan digunakan
truncated<-function(x,orde,k)
{
library(pracma)
orde=1
x<-as.vector(x)
k<-as.vector(k)
d<-length(k)
trun<-matrix(0,nrow=length(x),ncol=orde+d+1)
xtrun<-matrix(0,nrow=length(x),ncol=orde+d+1)
for (j in 1:(orde+1)){
xtrun[,j]<-x^(j-1)}
for (r in 1:d){
for (i in 1:length(x)){
trun[i,r+orde+1]<-ifelse(x[i]>=k[r],(x[i]-k[r])^(orde),0)
}
}
for (j in (orde+2):(orde+1+d)){
xtrun[,j]=trun[,j]
}
H=xtrun%*%pinv(t(xtrun)%*%xtrun)%*%t(xtrun)
hasil=list(xtrun=xtrun,H=H)
return(hasil)
}
mtruncated<-function(x,knot){
# fungsi membuat matriks x truncated untuk x=p-variabel
#knot=list dari k untuk masing masing variabel
library(pracma)
x=as.matrix(x)
p=ncol(x)
xp=vector("list",p)
for (r in 1:p){
xp[[r]]=truncated(x[,r],orde,knot[[r]])$xtrun[,-1]
}
xt=do.call("cbind",xp)
xtruncated=cbind(rep(1,nrow(x)),xt)
H=xtruncated%*%pinv(t(xtruncated)%*%xtruncated)%*%t(xtruncated)
hasil=list(xtruncated=xtruncated,H=H)
return(hasil)
}
gcv1knot <- function(x, y, b) {
x <- as.matrix(x)
y <- as.matrix(y)
gcv <- as.vector(0)
knot <- matrix(0, nrow = b, ncol = ncol(x))
knott <- matrix(0, nrow = b - 2, ncol = ncol(x))
knots <- vector("list", b)
for (j in 1:ncol(x)) {
knot[, j] <- seq(min(x[, j]), max(x[, j]), length.out = b)
knott[, j] <- knot[c(-1, -b), j]
}
kn <- vector("list", ncol(x))
# Perhitungan tanpa output ke konsol
for (i in 1:nrow(knott)) {
knots[[i]] <- knott[i, ]
kn <- as.list(knots[[i]])
H <- mtruncated(x, kn)$H
identity <- diag(1, nrow = nrow(x))
numerator <- (1 / nrow(x)) * t(y) %*% (identity - H) %*% y
denominator <- ((1 / nrow(x)) * sum(diag(identity - H)))^2
gcv[i] <- numerator / denominator
}
optimum <- min(gcv)
indeks.knotoptimum <- which.min(gcv)
knot.optimum <- knott[indeks.knotoptimum, ]
hasil <- list(knot = knott, gcv = gcv, gcv.minimum = optimum, knot.optimum = knot.optimum)
return(hasil)
}
gcv2knot <- function(x, y, b) {
x <- as.matrix(x)
y <- as.matrix(y)
gcv <- as.vector(0)
knot <- matrix(0, nrow = b, ncol = ncol(x))
knott <- matrix(0, nrow = b - 2, ncol = ncol(x))
kn <- vector("list", ncol(x))
a <- vector("list", ncol(x))
for (j in 1:ncol(x)) {
knot[, j] <- seq(min(x[, j]), max(x[, j]), length.out = b)
knott[, j] <- knot[c(-1, -b), j]
a[[j]] <- t(combn(knott[, j], 2))
}
total_iter <- nrow(a[[1]])
# Perhitungan tanpa output ke konsol
for (i in 1:total_iter) {
kn <- lapply(a, "[", i, 1:2)
H <- mtruncated(x, kn)$H
identity <- diag(1, nrow(x))
numerator <- (1 / nrow(x)) * t(y) %*% (identity - H) %*% y
denominator <- ((1 / nrow(x)) * sum(diag(identity - H)))^2
gcv[i] <- numerator / denominator
}
optimum <- min(gcv)
indeks.knotoptimum <- which.min(gcv)
knot.optimum <- lapply(a, "[", indeks.knotoptimum, 1:2)
hasil <- list(gcv = gcv, gcv.minimum = optimum, knot.optimum = knot.optimum)
return(hasil)
}
gcv3knot <- function(x, y, b) {
x <- as.matrix(x)
y <- as.matrix(y)
gcv <- as.vector(0)
knot <- matrix(0, nrow = b, ncol = ncol(x))
knott <- matrix(0, nrow = b - 2, ncol = ncol(x))
kn <- vector("list", ncol(x))
a <- vector("list", ncol(x))
for (j in 1:ncol(x)) {
knot[, j] <- seq(min(x[, j]), max(x[, j]), length.out = b)
knott[, j] <- knot[c(-1, -b), j]
a[[j]] <- t(combn(knott[, j], 3))
}
total_iter <- nrow(a[[1]])
# Perhitungan tanpa output ke konsol
for (i in 1:total_iter) {
kn <- lapply(a, "[", i, 1:3)
H <- mtruncated(x, kn)$H
identity <- diag(1, nrow(x))
numerator <- (1 / nrow(x)) * t(y) %*% (identity - H) %*% y
denominator <- ((1 / nrow(x)) * sum(diag(identity - H)))^2
gcv[i] <- numerator / denominator
}
optimum <- min(gcv)
indeks.knotoptimum <- which.min(gcv)
knot.optimum <- lapply(a, "[", indeks.knotoptimum, 1:3)
hasil <- list(gcv = gcv, gcv.minimum = optimum, knot.optimum = knot.optimum)
return(hasil)
}
kombinasi.knot <- function(x, y, b) {
x <- as.matrix(x)
y <- as.matrix(y)
gcv <- as.vector(0)
# Dapatkan hasil GCV untuk masing-masing jenis (tanpa output)
k1 <- gcv1knot(x, y, b)
k2 <- gcv2knot(x, y, b)
k3 <- gcv3knot(x, y, b)
# Simpan semua knot ke dalam list
knott <- vector("list", ncol(x))
for (i in 1:ncol(x)) {
knott[[i]] <- list(k1$knot.optimum[i], k2$knot.optimum[[i]], k3$knot.optimum[[i]])
}
# Buat semua kombinasi knot (1-2-3)
kom <- expand.grid(knott)
# Perhitungan tanpa output ke konsol
for (i in 1:nrow(kom)) {
kno <- lapply(kom, "[[", i)
H <- mtruncated(x, kno)$H
identity <- diag(1, nrow(x))
numerator <- (1 / nrow(x)) * t(y) %*% (identity - H) %*% y
denominator <- ((1 / nrow(x)) * sum(diag(identity - H)))^2
gcv[i] <- numerator / denominator
}
# Pilih kombinasi dengan GCV minimum
optimum <- min(gcv)
indeks.knotoptimum <- which.min(gcv)
knot.optimum <- lapply(kom, "[[", indeks.knotoptimum)
combo.optimum <- sapply(1:length(knott), function(j) {
if (identical(knot.optimum[[j]], knott[[j]][[1]])) {
return(1)
} else if (identical(knot.optimum[[j]], knott[[j]][[2]])) {
return(2)
} else {
return(3)
}
})
hasil <- list(
gcv = gcv,
gcv.1knot = k1$gcv.minimum,
gcv.2knot = k2$gcv.minimum,
gcv.3knot = k3$gcv.minimum,
satuknot.optimum = k1$knot.optimum,
duaknot.optimum = k2$knot.optimum,
tigaknot.optimum = k3$knot.optimum,
gcv.minimum = optimum,
knot.optimum = knot.optimum,
kombinasi.knot.optimum = combo.optimum
)
return(hasil)
}fungsi untuk penaksiran beta dan pengujian hipotesis pada suatu titik knot yang ditentukan dan suatu taraf alpha.sedangkan knot adalah suatu list dari vektor knot maksimal 3 knot untuk masing masing variabel
hitung<-function(x,y,knot,taraf.alpha){
library(pracma)
library(car)
x=as.matrix(x)
y=as.matrix(y)
I=diag(1,nrow=nrow(x),ncol=nrow(x))
J=matrix(1,nrow=nrow(x),ncol=nrow(x))
matriksX=mtruncated(x,knot)$xtruncated
Ak=mtruncated(x,knot)$H
beta=pinv(t(matriksX)%*%matriksX)%*%t(matriksX)%*%y
SSE=t(y)%*%(I-Ak)%*%y
SST=t(y)%*%(I-(1/nrow(x))*J)%*%y
SSR=SST-SSE
R.square=as.numeric((SSR/SST)*100)
dbr=ncol(matriksX)-1
dbe=nrow(x)-ncol(matriksX)
dbt=nrow(x)-1
MSE=SSE/dbe
MSR=SSR/dbr
Se=sqrt(MSE*diag(pinv(t(matriksX)%*%matriksX)))
Fhitung=MSR/MSE
p.valueF=pf(Fhitung,dbr,dbe,lower.tail=FALSE )
keputusan=NULL
keputusan.1=NULL
if (p.valueF >= taraf.alpha)
{keputusan.1=c("gagal tolak H0")}
else {keputusan.1=c("tolak H0")}
for (i in 1:length(Se)){
thitung=beta/Se
p.value=2*pt(abs(thitung),df=nrow(x)-1,lower=FALSE)
if (p.value[i]>=taraf.alpha) {
keputusan[i]=c("gagal tolak H0")
} else{
keputusan[i]=c("tolak H0")
}
}
ytopi=Ak%*%y
b<-data.frame(Sumber=c("regresi","error","total"),Sum.of.Square=c(round(SSR,4),round(SSE,4),round(SST,4)),db=c(dbr,dbe,dbt),mean.of.square=c(round(MSR,4),round(MSE,4),"-"),Fhitung=c(round(Fhitung,4),"-","-"),p.value=c(p.valueF,"-","-"),keputusan=c(keputusan.1,"-","-"))
inferensi=data.frame(beta=beta,standar.error=Se,t.hitung=thitung,p.value=p.value,keputusan=keputusan)
hasil<-list(ytopi=ytopi,beta=beta,tabel.anova=b,inferensi=inferensi,R.square=R.square)
return(hasil)
}Fungsi untuk spline linier multivariabel. Fungsi ini memanggil fungsi fungsi yang lain dalam membentuk matriks truncated, penaksiran parameter, dan mencari gcv minimum. b adalah banyaknya titik yang dicoba sebagai knot untuk masing masing variabel dalam trial error pencarian knot optimum
spline.truncated.linier.multivariabel=function(x,y,b,taraf.alpha){
gcv=kombinasi.knot(x,y,b)
gcvsatuknot=gcv$gcv.1knot
satuknot=gcv$satuknot.optimum
gcvduaknot=gcv$gcv.2knot
duaknot=gcv$duaknot.optimum
gcvtigaknot=gcv$gcv.3knot
tigaknot=gcv$tigaknot.optimum
knot.optimum=gcv$knot.optimum
gcv.minimum=gcv$gcv.minimum
beta=hitung(x,y,knot.optimum,taraf.alpha)$beta
tabel.anova=hitung(x,y,knot.optimum,taraf.alpha)$tabel.anova
rsquare=hitung(x,y,knot.optimum,taraf.alpha)$R.square
inferensi=hitung(x,y,knot.optimum,taraf.alpha)$inferensi
ytopi=hitung(x,y,knot.optimum,taraf.alpha)$ytopi
residual=y-ytopi
kenormalan=ks.test(residual,"pnorm",mean=mean(residual),sd=sd(residual))
if(kenormalan$p.value >= taraf.alpha){
keputusan.kenormalan=c("residual berdistribusi normal")}
else { keputusan.kenormalan=c("residual tidak berdistribusi normal")}
gleyjser=hitung(x,abs(residual),knot.optimum,taraf.alpha)$tabel.anova
uji.kenormalan.KS=data.frame(D=as.numeric(kenormalan$statistic),p.value=kenormalan$p.value,keputusan=keputusan.kenormalan)
durbin=durbinWatsonTest(as.vector(residual))
matrikstruncated=mtruncated(x,knot.optimum)$xtruncated
p.value.dw=durbinWatsonTest(lm(y~matrikstruncated[,-1]))$p
if(p.value.dw >= taraf.alpha){
keputusan.dw=c("gagal tolak H0")}
else { keputusan.dw=c("tolak H0")}
cat("=================================================","\n")
cat("satu knot optimum untuk masing-masing variabel","\n")
cat("=================================================","\n")
cat("","\n")
cat("GCV minimum = ",gcvsatuknot,"\n")
cat("","\n")
for (i in 1:length(satuknot)){
cat("knot optimum variabel ke-",i,"\n")
cat(satuknot[[i]],"\n")
cat("","\n")
}
cat("=================================================","\n")
cat("dua knot optimum untuk masing-masing variabel","\n")
cat("=================================================","\n")
cat("","\n")
cat("GCV minimum = ",gcvduaknot,"\n")
cat("","\n")
for (i in 1:length(duaknot)){
cat("knot optimum variabel ke-",i,"\n")
cat(duaknot[[i]],"\n")
cat("","\n")
}
cat("=================================================","\n")
cat("tiga knot optimum untuk masing-masing variabel","\n")
cat("=================================================","\n")
cat("","\n")
cat("GCV minimum = ",gcvtigaknot,"\n")
cat("","\n")
for (i in 1:length(tigaknot)){
cat("knot optimum variabel ke-",i,"\n")
cat(tigaknot[[i]],"\n")
cat("","\n")
}
cat("==============================================","\n")
cat("hasil kombinasi knot optimum dan GCV minimum","\n")
cat("==============================================","\n")
cat("GCV minimum = ",gcv.minimum,"\n")
cat("","\n")
cat("knot optimum","\n")
for (i in 1:length(knot.optimum)){
cat("knot optimum variabel ke-",i,"\n")
cat(knot.optimum[[i]],"\n")
cat("","\n")
}
cat("","\n")
cat("=====================================================================================================================","\n")
cat("tabel an0va","\n")
cat("=====================================================================================================================","\n")
cat("Sumber db SS MS Fhit P-value keputusan","\n")
cat("=====================================================================================================================","\n")
cat("Regresi ",tabel.anova[1,3]," ",round(tabel.anova[1,2],3)," ",as.character(tabel.anova[1,4])," ",as.character(tabel.anova[1,5])," ",as.character(tabel.anova[1,6])," ",as.character(tabel.anova[1,7]),"\n")
cat("Error ",tabel.anova[2,3]," ",round(tabel.anova[2,2],3)," ",as.character(tabel.anova[2,4]),"\n")
cat("Total ",tabel.anova[3,3]," ",round(tabel.anova[3,2],3),"\n")
cat("=====================================================================================================================","\n")
cat("","\n")
cat("","\n")
cat("=====================================================================================================================","\n")
cat("parameter beta dan uji individu","\n")
cat("=====================================================================================================================","\n")
cat(" beta standar error t hitung p-value keputusan","\n")
cat("=====================================================================================================================","\n")
for (i in 1:nrow(beta)){
cat("beta",i-1 ," ",beta[i]," ",inferensi[i,2]," ",inferensi[i,3]," ",inferensi[i,4]," ",as.character(inferensi[i,5]),"\n")
}
cat("=====================================================================================================================","\n")
cat("","\n")
cat("R square =",rsquare,"\n")
cat("","\n")
cat("","\n")
cat("=====================================================================================================================","\n")
cat("Uji kolmogorov-smirnov untuk kenormalan residual ","\n")
cat("=====================================================================================================================","\n")
cat("D P-value keputusan","\n")
cat("=====================================================================================================================","\n")
cat(as.numeric(kenormalan$statistic)," ",kenormalan$p.value," ",keputusan.kenormalan,"\n")
cat("=====================================================================================================================","\n")
cat("","\n")
cat("","\n")
cat("","\n")
cat("=====================================================================================================================","\n")
cat("uji Gleyjser","\n")
cat("=====================================================================================================================","\n")
cat("Sumber db SS MS Fhit P-value keputusan","\n")
cat("=====================================================================================================================","\n")
cat("Regresi ",gleyjser[1,3]," ",round(gleyjser[1,2],3)," ",as.character(gleyjser[1,4])," ",as.character(gleyjser[1,5])," ",as.character(gleyjser[1,6])," ",as.character(gleyjser[1,7]),"\n")
cat("Error ",gleyjser[2,3]," ",round(gleyjser[2,2],3)," ",as.character(gleyjser[2,4]),"\n")
cat("Total ",gleyjser[3,3]," ",round(gleyjser[3,2],3),"\n")
cat("=====================================================================================================================","\n")
cat("","\n")
cat("","\n")
cat("=====================================================================================================================","\n")
cat("uji Durbin Watson ","\n")
cat("=====================================================================================================================","\n")
cat("DW P-value keputusan","\n")
cat("=====================================================================================================================","\n")
cat(as.numeric(durbin)," ",p.value.dw," ",keputusan.dw,"\n")
cat("=====================================================================================================================","\n")
hasil=list(beta=beta,ytopi=ytopi)
}x=data[,c(2,3,4,5,6)] #menggunakan variabel X2 dan X3
y=data[,1]
model=spline.truncated.linier.multivariabel(x,y,b=50,taraf.alpha=0.05)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:pracma':
##
## logit## Warning in MSE * diag(pinv(t(matriksX) %*% matriksX)): Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.
## Warning in MSE * diag(pinv(t(matriksX) %*% matriksX)): Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.
## Warning in MSE * diag(pinv(t(matriksX) %*% matriksX)): Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.
## Warning in MSE * diag(pinv(t(matriksX) %*% matriksX)): Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.
## Warning in MSE * diag(pinv(t(matriksX) %*% matriksX)): Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.
## Warning in MSE * diag(pinv(t(matriksX) %*% matriksX)): Recycling array of length 1 in array-vector arithmetic is deprecated.
## Use c() or as.vector() instead.## =================================================
## satu knot optimum untuk masing-masing variabel
## =================================================
##
## GCV minimum = 3.103661
##
## knot optimum variabel ke- 1
## 62.69286
##
## knot optimum variabel ke- 2
## 68.79673
##
## knot optimum variabel ke- 3
## 80.05143
##
## knot optimum variabel ke- 4
## 40.39224
##
## knot optimum variabel ke- 5
## 4.839388
##
## =================================================
## dua knot optimum untuk masing-masing variabel
## =================================================
##
## GCV minimum = 1.173644
##
## knot optimum variabel ke- 1
## 65.66429 79.40714
##
## knot optimum variabel ke- 2
## 71.26367 82.67327
##
## knot optimum variabel ke- 3
## 84.97714 107.7586
##
## knot optimum variabel ke- 4
## 41.88122 48.76776
##
## knot optimum variabel ke- 5
## 5.676939 9.550612
##
## =================================================
## tiga knot optimum untuk masing-masing variabel
## =================================================
##
## GCV minimum = 0.1821971
##
## knot optimum variabel ke- 1
## 62.69286 65.29286 69.75
##
## knot optimum variabel ke- 2
## 68.79673 70.95531 74.65571
##
## knot optimum variabel ke- 3
## 80.05143 84.36143 91.75
##
## knot optimum variabel ke- 4
## 40.39224 41.6951 43.92857
##
## knot optimum variabel ke- 5
## 4.839388 5.572245 6.828571
##
## ==============================================
## hasil kombinasi knot optimum dan GCV minimum
## ==============================================
## GCV minimum = 0.1821971
##
## knot optimum
## knot optimum variabel ke- 1
## 62.69286 65.29286 69.75
##
## knot optimum variabel ke- 2
## 68.79673 70.95531 74.65571
##
## knot optimum variabel ke- 3
## 80.05143 84.36143 91.75
##
## knot optimum variabel ke- 4
## 40.39224 41.6951 43.92857
##
## knot optimum variabel ke- 5
## 4.839388 5.572245 6.828571
##
##
## =====================================================================================================================
## tabel an0va
## =====================================================================================================================
## Sumber db SS MS Fhit P-value keputusan
## =====================================================================================================================
## Regresi 20 106.046 5.3023 96.2151 6.46969095438813e-06 tolak H0
## Error 6 0.331 0.0551
## Total 26 106.377
## =====================================================================================================================
##
##
## =====================================================================================================================
## parameter beta dan uji individu
## =====================================================================================================================
## beta standar error t hitung p-value keputusan
## =====================================================================================================================
## beta 0 -0.1316722 0.02487688 -5.292956 1.553596e-05 tolak H0
## beta 1 2.734932 0.7324064 3.734173 0.0009314847 tolak H0
## beta 2 -2.497948 0.7550908 -3.308142 0.002751953 tolak H0
## beta 3 -0.6079231 0.1434221 -4.238699 0.0002502843 tolak H0
## beta 4 3.083237 0.3364033 9.165302 1.26054e-09 tolak H0
## beta 5 9.206173 1.026016 8.972737 1.926392e-09 tolak H0
## beta 6 -7.48363 0.9762012 -7.666074 3.897513e-08 tolak H0
## beta 7 -2.371142 0.3270661 -7.249731 1.064894e-07 tolak H0
## beta 8 0.806834 0.1006068 8.019679 1.689227e-08 tolak H0
## beta 9 -4.586374 1.092282 -4.19889 0.0002778525 tolak H0
## beta 10 13.95131 1.960695 7.115493 1.479353e-07 tolak H0
## beta 11 -9.335357 1.000364 -9.331962 8.766981e-10 tolak H0
## beta 12 0.07329268 0.05882518 1.24594 0.2238969 gagal tolak H0
## beta 13 -12.85148 1.123454 -11.43925 1.20149e-11 tolak H0
## beta 14 13.40361 1.305824 10.26448 1.22763e-10 tolak H0
## beta 15 0.3745942 0.5909591 0.6338751 0.5316999 gagal tolak H0
## beta 16 -0.8829472 0.2750863 -3.20971 0.003517033 tolak H0
## beta 17 8.847561 1.10683 7.993604 1.79564e-08 tolak H0
## beta 18 -5.136792 1.335426 -3.846556 0.0006966888 tolak H0
## beta 19 -14.99066 1.100068 -13.62703 2.382174e-13 tolak H0
## beta 20 15.50645 0.9948482 15.58675 1.046411e-14 tolak H0
## =====================================================================================================================
##
## R square = 99.68917
##
##
## =====================================================================================================================
## Uji kolmogorov-smirnov untuk kenormalan residual
## =====================================================================================================================
## D P-value keputusan
## =====================================================================================================================
## 0.1451742 0.5703618 residual berdistribusi normal
## =====================================================================================================================
##
##
##
## =====================================================================================================================
## uji Gleyjser
## =====================================================================================================================
## Sumber db SS MS Fhit P-value keputusan
## =====================================================================================================================
## Regresi 20 0.116 0.0058 1.3249 0.386784976407208 gagal tolak H0
## Error 6 0.026 0.0044
## Total 26 0.142
## =====================================================================================================================
##
##
## =====================================================================================================================
## uji Durbin Watson
## =====================================================================================================================
## DW P-value keputusan
## =====================================================================================================================
## 2.19786 0.3 gagal tolak H0
## =====================================================================================================================library(openxlsx)
# === Buat workbook ===
wb <- createWorkbook()
# === Data GCV 1 Knot ===
k1 <- gcv1knot(x, y, b = 50)
df_k1 <- as.data.frame(k1$knot)
colnames(df_k1) <- paste0("Knot_Var", 1:ncol(x))
df_k1$GCV <- round(k1$gcv, 6)
addWorksheet(wb, "GCV - 1 Knot")
writeData(wb, "GCV - 1 Knot", df_k1)
# === Data GCV 2 Knot ===
k2 <- gcv2knot(x, y, b = 50)
a2 <- vector("list", ncol(x))
for (j in 1:ncol(x)) {
a2[[j]] <- t(combn(seq(min(x[, j]), max(x[, j]), length.out = 50)[-c(1, 50)], 2))
}
df_k2 <- data.frame(GCV = round(k2$gcv, 6))
for (j in 1:ncol(x)) {
df_k2[[paste0("Knot_Var", j)]] <- apply(a2[[j]], 1, function(row) paste(round(row, 5), collapse = " | "))
}
addWorksheet(wb, "GCV - 2 Knot")
writeData(wb, "GCV - 2 Knot", df_k2)
# === Data GCV 3 Knot ===
k3 <- gcv3knot(x, y, b = 50)
a3 <- vector("list", ncol(x))
for (j in 1:ncol(x)) {
a3[[j]] <- t(combn(seq(min(x[, j]), max(x[, j]), length.out = 50)[-c(1, 50)], 3))
}
df_k3 <- data.frame(GCV = round(k3$gcv, 6))
for (j in 1:ncol(x)) {
df_k3[[paste0("Knot_Var", j)]] <- apply(a3[[j]], 1, function(row) paste(round(row, 5), collapse = " | "))
}
addWorksheet(wb, "GCV - 3 Knot")
writeData(wb, "GCV - 3 Knot", df_k3)
# === Data Kombinasi Knot ===
kombi <- kombinasi.knot(x, y, b = 50)
knott <- vector("list", ncol(x))
for (i in 1:ncol(x)) {
knott[[i]] <- list(k1$knot.optimum[i], k2$knot.optimum[[i]], k3$knot.optimum[[i]])
}
kom <- expand.grid(knott)
df_komb <- data.frame(GCV = round(kombi$gcv, 6))
for (j in 1:ncol(x)) {
df_komb[[paste0("Knot_Var", j)]] <- sapply(1:nrow(kom), function(i) paste(round(kom[[j]][[i]], 5), collapse = " | "))
}
addWorksheet(wb, "GCV - Kombinasi")
writeData(wb, "GCV - Kombinasi", df_komb)
# === Simpan Excel ===
saveWorkbook(wb, "Hasil_Iterasi_GCV_Spline_tpt.xlsx", overwrite = TRUE)