#MEMANGGIL LIBRARY YANG DIGUNAKAN#

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

#MEMANGGIL DATA#

data=read.table(file.choose(),header=TRUE)
data

#ANALISIS STATISTIKA DESKRIPTIF & DETEKSI MULTIKOLINIERITAS#

stat.desc(data)

#SCATTERPLOT#

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)

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

#PEMILIHAN 1 TITIK KNOT#

GCV1=function(data) 
{ 
  para=0 
  data=as.matrix(data) 
  N=nrow(data) 
  M=ncol(data) 
  m=ncol(data)-para-1 
  dataA=data[,(para+2):M] 
  dataA=as.matrix(dataA) 
  F=diag(N) 
  nk=50 #nk=banyaknya alternatif titik knot yang akan dicoba  
  knot1=matrix(ncol=m,nrow=nk) 
  for (i in (1:m)) #membuat knot 
  { 
    a=seq(min(dataA[,i]),max(dataA[,i]),length.out=nk) 
    knot1[,i]=t(as.matrix(a)) 
  } 
  a1=length(knot1[,1]) 
  knot1=as.matrix(knot1[2:(a1-1),]) 
  aa=rep(1,N) 
  data1=matrix(ncol=m,nrow=N) 
  data2=data[,2:M] #data x saja  
  nk1=nrow(knot1) 
  GCV=as.matrix(rep(NA,nk1),ncol=1);colnames(GCV)<-"GCV" 
  MSE=as.matrix(rep(NA,nk1),ncol=1);colnames(MSE)<-"MSE" 
  SSE=rep(NA,nk1) 
  SSR=rep(NA,nk1)  
  Rsq=as.matrix(rep(NA,nk1),ncol=1);colnames(Rsq)<-"Rsq" 
  knotke=matrix(c(1:nk1),ncol=1);colnames(knotke)<-"knot_ke" 
  for (i in 1:nk1) 
  { 
    for (j in 1:m)  
    {  
      for (k in 1:N)  
      {  
        if (data[k,(j+para+1)]<knot1[i,j]) data1[k,j]=0 else  
          data1[k,j]=data[k,(j+para+1)]-knot1[i,j] 
      } 
    } 
    mx=as.matrix(cbind(aa,data2,data1)) 
    C=pinv(t(mx)%*%mx) 
    B=C%*%(t(mx)%*%data[,1]) 
    yhat=mx%*%B 
    res=data[,1]-yhat 
    SSE[i]=sum((res)^2) 
    SSR[i]=sum((yhat-mean(data[,1]))^2) 
    MSE[i]=SSE[i]/(N) 
    Rsq[i]=(SSR[i]/(SSR[i]+SSE[i]))*100 
    A=mx%*%C%*%t(mx) 
    A1=(F-A) 
    A2=(sum(diag(A1))/N)^2 
    GCV[i]=MSE[i]/A2 
  } 
  dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot1)) 
  dataG=dataAll[order(GCV),-2] 
  write.csv(dataAll, file="C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\DATA ALL knot 1.csv")
  cat("==============================================","\n") 
  cat("HASIL GCV terkecil dengan 1 knot","\n") 
  cat("==============================================","\n") 
  print(((dataG[1,1:6]))) 
  cat("Nilai GCV 10 terkecil pertama","\n") 
  print(dataG[1:10,]) 
  mingcv=dataG[1,1] 
  knotgcv=as.matrix(knot1[dataG[1,4],]) 
  knotgcv1=matrix(knotgcv,nrow=1) 
  datagcv1=matrix(ncol=m,nrow=N)  
  for (j in 1:m)  
  {  
    for (k in 1:N)  
    {  
      if (data[k,(j+para+1)]<knotgcv[j,1]) datagcv1[k,j]=0 else  
        datagcv1[k,j]=data[k,(j+para+1)]-knotgcv[j,1]  
    }  
  } 
  mxgcv=as.matrix(cbind(aa,data2,datagcv1))  
  mxgcv=mxgcv[,c(2:6)] 
  list(knotgcv=knotgcv1,mingcv=mingcv,mxgcv=mxgcv) 
  cat("\n")  
  cat("==============================================","\n") 
  cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 1","\n") 
  cat("==============================================","\n") 
  print(B) 
  cat("\n") 
} 
GCV1(data)
## ============================================== 
## HASIL GCV terkecil dengan 1 knot 
## ============================================== 
##       GCV   knot_ke                                         
##  1.508187 45.000000 76.286735 14.630612 69.183673 75.922653 
## Nilai GCV 10 terkecil pertama 
##            GCV knot_ke                                     
##  [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
## 
## ============================================== 
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 1 
## ============================================== 
##                [,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

#UJI SIGNIFIKANSI SATU TITIK KNOT#

uji=function(alpha,para) 
{ 
  # --- baca & siapkan titik knot terbaik (GCV terkecil) ---
  knot = read.csv("C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\DATA ALL knot 1.csv", header=TRUE)
  knot = as.matrix(knot)
  knot = knot[,-1]                 # buang kolom index baris
  knot = knot[order(knot[,1]), ]   # urutkan berdasarkan GCV (kolom 1) terkecil
  baris_knot = knot[1, 4:7]        # ambil knot X1-X4 pada baris GCV terkecil
  
  data=as.matrix(data) 
  ybar=mean(data[,1]) 
  m=para+2 
  n=nrow(data) 
  q=ncol(data) 
  dataA=cbind(data[,m],data[,m+1], 
              data[,m+2],data[,m+3]) 
  dataA=as.matrix(dataA) 
  satu=rep(1,n) 
  n1=4                              # jumlah variabel X (X1-X4)
  data.knot=matrix(ncol=n1,nrow=n) 
  
  # --- hitung fungsi basis knot (loop i untuk variabel, j untuk observasi) ---
  for (i in 1:n1) 
  { 
    for (j in 1:n) 
    { 
      if(dataA[j,i] < baris_knot[i]) data.knot[j,i]=0 
      else data.knot[j,i] = dataA[j,i] - baris_knot[i] 
    } 
  } 
  
  mx=cbind(satu,data[,2],data.knot[,1:1],data[,3], 
           data.knot[,2:2],data[,4],data.knot[,3:3], 
           data[,5],data.knot[,4:4]) 
  mx=as.matrix(mx) 
  B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1] 
  n1=nrow(B) 
  yhat=mx%*%B 
  ybar=mean(data[,1]) 
  res=data[,1]-yhat 
  SSE=sum((data[,1]-yhat)^2) 
  SSR=sum((yhat-ybar)^2) 
  MSE=SSE/(n) 
  MSR=SSR/(n1-1) 
  SST=sum((data[,1]-ybar)^2) 
  Rsq=(SSR/(SSR+SSE))*100 
  
  #-------------------------------------------------------# 
  # UJI SIMULTAN 
  #-------------------------------------------------------# 
  Fhit=MSR/MSE 
  pvalue=pf(Fhit,(n1-1),(n-n1),lower.tail=FALSE) 
  if(pvalue<=alpha) 
  { 
    cat('---------------------------------------','\n') 
    cat('Kesimpulan hasil uji simultan','\n') 
    cat('---------------------------------------','\n') 
    cat('Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan','\n') 
    cat('','\n') 
  } 
  else                                # <-- perbaikan: "Else" jadi "else"
  { 
    cat('---------------------------------------','\n') 
    cat('Kesimpulan hasil uji simultan','\n') 
    cat('---------------------------------------','\n') 
    cat('Gagal Tolak Ho yakni semua variabel bebas tidak berpengaruh signifikan','\n') 
    cat('','\n') 
  } 
  
  #------------------------------------------------------# 
  # UJI PARSIAL 
  #------------------------------------------------------# 
  thit=rep(NA,n1) 
  pval=rep(NA,n1) 
  SE=sqrt(diag(MSE*(pinv(t(mx)%*%mx)))) 
  cat('---------------------------------------------','\n') 
  cat('Kesimpulan hasil uji parsial','\n') 
  cat('---------------------------------------------','\n') 
  for (i in 1:n1) 
  { 
    thit[i]=B[i,1]/SE[i] 
    pval[i]=2*(pt(abs(thit[i]),(n-n1),lower.tail=FALSE)) 
    if (pval[i]<=alpha) cat('Tolak Ho yakni variabel bebas signifikan dengan pvalue',pval[i],'\n') 
    else cat('Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue',pval[i],'\n') 
  } 
  thit=as.matrix(thit) 
  cat('=============================================','\n') 
  cat('nilai t hitung','\n') 
  cat('=============================================','\n') 
  print(thit) 
  cat('Analysis of Variance','\n') 
  cat('=============================================','\n') 
  cat('Sumber df SS MS Fhit','\n') 
  cat('Regresi ',(n1-1),' ',SSR,' ',MSR,' ',Fhit,'\n') 
  cat('Error ',n-n1,' ',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='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI RESIDUAL knot1.csv') 
  write.csv(mx, file='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI MX knot1.csv') 
  write.csv(yhat, file='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI yhat knot1.csv') 
}

uji(0.1, 0)
## --------------------------------------- 
## Kesimpulan hasil uji simultan 
## --------------------------------------- 
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan 
##  
## --------------------------------------------- 
## Kesimpulan hasil uji parsial 
## --------------------------------------------- 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.04471309 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 6.825617e-06 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.9079879 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.006735878 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 2.184894e-06 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.5049048 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.07036832 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.09074499 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 4.259298e-07 
## ============================================= 
## nilai t hitung 
## ============================================= 
##             [,1]
##  [1,]  2.0624803
##  [2,] -5.0617196
##  [3,]  0.1162015
##  [4,]  2.8350060
##  [5,] -5.3953305
##  [6,] -0.6719509
##  [7,]  1.8516205
##  [8,] -1.7269469
##  [9,]  5.8681386
## Analysis of Variance 
## ============================================= 
## Sumber df SS MS Fhit 
## Regresi  8   167.6864   20.9608   19.73027 
## Error  47   59.49258   1.062368 
## Total  55   227.1789 
## ============================================= 
## s= 1.030712  Rsq= 73.81246 
## pvalue(F)= 1.239e-12

#PEMILIHAN DUA TITIK KNOT#

GCV2=function(data)  
{  
  para=0 
  data=as.matrix(data)  
  N=nrow(data) 
  M=ncol(data) 
  m=ncol(data)-para-1 
  dataA=data[,(para+2):M]  
  dataA=as.matrix(dataA)  
  F=diag(N)  
  nk=50 
  knot1=matrix(ncol=m,nrow=nk)  
  for (i in (1:m)) 
  { 
    a=seq(min(dataA[,i]),max(dataA[,i]),length.out=nk)  
    knot1[,i]=t(as.matrix(a)) 
  } 
  a1=length(knot1[,1]) 
  knot1=as.matrix(knot1[2:(a1-1),]) 
  a2=nk-2  
  z=(a2*(a2-1)/2) 
  knot2=cbind(rep(NA,(z+1))) 
  for (i in (1:m)) 
  { 
    knot=rbind(rep(NA,2)) 
    for ( j in 1:(a2-1)) 
    { 
      for (k in (j+1):a2)  
      { xx=cbind(knot1[j,i],knot1[k,i]) 
      knot=rbind(knot,xx) 
      } 
    } 
    knot2=cbind(knot2,knot)  
  } 
  knot2=knot2[2:(z+1),2:(2*m+1)] 
  a3=nrow(knot2) 
  aa=rep(1,N) 
  data1=matrix(ncol=2*m,nrow=N) 
  data2=data[,2:M] 
  nk1=nrow(knot2) 
  GCV=as.matrix(rep(NA,nk1),ncol=1);colnames(GCV)<-"GCV"  
  MSE=as.matrix(rep(NA,nk1),ncol=1);colnames(MSE)<-"MSE" 
  SSE=rep(NA,nk1) 
  SSR=rep(NA,nk1)  
  Rsq=as.matrix(rep(NA,nk1),ncol=1);colnames(Rsq)<-"Rsq"  
  knotke=matrix(c(1:nk1),ncol=1);colnames(knotke)<-"knot_ke" 
  for (i in 1:a3)  
  {  
    for (j in 1:(2*m))  
    {  
      if (mod(j,2)==1) b=floor(j/2)+1 else b=j/2  
      for (k in 1:N)  
      { 
        if (data[k,(b+para+1)]<knot2[i,j]) data1[k,j]=0 else  
          data1[k,j]=data[k,(b+para+1)]-knot2[i,j] 
      } 
    } 
    mx=as.matrix(cbind(aa,data2,data1)) 
    C=pinv(t(mx)%*%mx) 
    B=C%*%(t(mx)%*%data[,1]) 
    yhat=mx%*%B 
    res=data[,1]-yhat 
    SSE[i]=sum((res)^2) 
    SSR[i]=sum((yhat-mean(data[,1]))^2) 
    MSE[i]=SSE[i]/(N) 
    Rsq[i]=(SSR[i]/(SSR[i]+SSE[i]))*100 
    A=mx%*%C%*%t(mx) 
    A1=(F-A) 
    A2=(sum(diag(A1))/N)^2 
    GCV[i]=MSE[i]/A2
  } 
  dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot2)) 
  dataG=dataAll[order(GCV),-2] 
  write.csv(dataAll,file="C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\DATA ALL knot 2.csv") 
  cat("==============================================","\n") 
  cat("HASIL GCV terkecil dengan 2 knot","\n")  
  cat("==============================================","\n") 
  print(((dataG[1,1:10])))  
  cat("Nilai GCV 10 terkecil pertama","\n") 
  print(dataG[1:10,]) 
  
  # --- PERBAIKAN: hitung ulang B dari titik knot dengan GCV TERKECIL ---
  # (sebelumnya B yang dicetak = B dari kombinasi knot TERAKHIR yang dicoba, bukan yang terbaik)
  baris_terbaik = dataG[1, "knot_ke"] 
  knotgcv = knot2[baris_terbaik, ] 
  datagcv1=matrix(ncol=2*m,nrow=N)  
  for (j in 1:(2*m)) 
  {  
    if (mod(j,2)==1) b=floor(j/2)+1 else b=j/2  
    for (k in 1:N)  
    { 
      if (data[k,(b+para+1)]<knotgcv[j]) datagcv1[k,j]=0 else 
        datagcv1[k,j]=data[k,(b+para+1)]-knotgcv[j] 
    } 
  } 
  mxgcv=as.matrix(cbind(aa,data2,datagcv1)) 
  Cgcv=pinv(t(mxgcv)%*%mxgcv) 
  B=Cgcv%*%(t(mxgcv)%*%data[,1]) 
  
  cat("\n")  
  cat("==============================================","\n") 
  cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 2","\n") 
  cat("==============================================","\n") 
  print(B) 
  cat("\n") 
} 
GCV2(data)
## ============================================== 
## HASIL GCV terkecil dengan 2 knot 
## ============================================== 
##        GCV    knot_ke                                                        
##   1.316068 135.000000  61.132449  76.286735  11.382041  14.630612   7.692245 
##                                  
##  69.183673  25.025510  75.922653 
## Nilai GCV 10 terkecil pertama 
##            GCV knot_ke                                                       
##  [1,] 1.316068     135 61.13245 76.28673 11.38204 14.63061  7.692245 69.18367
##  [2,] 1.317894     179 61.49327 76.28673 11.45939 14.63061  9.156327 69.18367
##  [3,] 1.321114     178 61.49327 75.92592 11.45939 14.55327  9.156327 67.71959
##  [4,] 1.321402     134 61.13245 75.92592 11.38204 14.55327  7.692245 67.71959
##  [5,] 1.358738     221 61.85408 75.92592 11.53673 14.55327 10.620408 67.71959
##  [6,] 1.359630     222 61.85408 76.28673 11.53673 14.63061 10.620408 69.18367
##  [7,] 1.359989     138 61.13245 77.36918 11.38204 14.86265  7.692245 73.57592
##  [8,] 1.388799     136 61.13245 76.64755 11.38204 14.70796  7.692245 70.64776
##  [9,] 1.393180     180 61.49327 76.64755 11.45939 14.70796  9.156327 70.64776
## [10,] 1.396087     177 61.49327 75.56510 11.45939 14.47592  9.156327 66.25551
##                        
##  [1,] 25.02551 75.92265
##  [2,] 26.23735 75.92265
##  [3,] 26.23735 74.71082
##  [4,] 25.02551 74.71082
##  [5,] 27.44918 74.71082
##  [6,] 27.44918 75.92265
##  [7,] 25.02551 79.55816
##  [8,] 25.02551 77.13449
##  [9,] 26.23735 77.13449
## [10,] 26.23735 73.49898
## 
## ============================================== 
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 2 
## ============================================== 
##              [,1]
##  [1,]  95.3587175
##  [2,]  -2.9213119
##  [3,]   6.1469387
##  [4,]   0.2059982
##  [5,]   0.7283636
##  [6,]   2.7449916
##  [7,]   0.7954129
##  [8,]  -5.3197766
##  [9,] -18.0705567
## [10,]  -0.2144538
## [11,]   0.2853691
## [12,]  -0.7619877
## [13,]   0.9034632

#UJI SIGNIFIKANSI DUA TITIK KNOT#

uji=function(alpha,para) 
{ 
  knot = read.csv("C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\DATA ALL knot 2.csv", header=TRUE)
  knot = as.matrix(knot) 
  knot = knot[,-1]                        # buang kolom index baris 
  knot = knot[order(knot[,1]), ]          # urutkan berdasarkan GCV terkecil 
  baris_knot = as.numeric(knot[1, 4:11])  # 8 nilai knot terbaik (X1a,X1b,X2a,X2b,X3a,X3b,X4a,X4b) 
  
  data=as.matrix(data) 
  ybar=mean(data[,1]) 
  m=para+2 
  n=nrow(data) 
  q=ncol(data) 
  dataA=cbind(data[,m],data[,m],data[,m+1],data[,m+1], 
              data[,m+2],data[,m+2],data[,m+3],data[,m+3]) 
  dataA=as.matrix(dataA) 
  satu=rep(1,n) 
  n1=length(baris_knot)                   # = 8 
  data.knot=matrix(ncol=n1,nrow=n) 
  for (i in 1:n1) 
  { 
    for (j in 1:n) 
    { 
      if(dataA[j,i]<baris_knot[i]) data.knot[j,i]=0 
      else data.knot[j,i]=dataA[j,i]-baris_knot[i] 
    } 
  } 
  mx=cbind(satu,data[,2],data.knot[,1:2],data[,3], 
           data.knot[,3:4],data[,4],data.knot[,5:6], 
           data[,5],data.knot[,7:8]) 
  mx=as.matrix(mx) 
  B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1] 
  n1=nrow(B) 
  yhat=mx%*%B 
  ybar=mean(data[,1]) 
  res=data[,1]-yhat 
  SSE=sum((data[,1]-yhat)^2)
  SSR=sum((yhat-ybar)^2) 
  MSE=SSE/(n) 
  MSR=SSR/(n1-1) 
  SST=sum((data[,1]-ybar)^2) 
  Rsq=(SSR/(SSR+SSE))*100 
  Fhit=MSR/MSE 
  pvalue=pf(Fhit,(n1-1),(n-n1),lower.tail=FALSE) 
  if(pvalue<=alpha) 
  { 
    cat('---------------------------------------','\n') 
    cat('Kesimpulan hasil uji simultan','\n') 
    cat('---------------------------------------','\n') 
    cat('Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan','\n') 
    cat('','\n') 
  } 
  else 
  { 
    cat('---------------------------------------','\n') 
    cat('Kesimpulan hasil uji simultan','\n') 
    cat('---------------------------------------','\n') 
    cat('Gagal Tolak Ho yakni semua variabel bebas tidak berpengaruh signifikan','\n') 
    cat('','\n') 
  } 
  thit=rep(NA,n1) 
  pval=rep(NA,n1) 
  SE=sqrt(diag(MSE*(pinv(t(mx)%*%mx)))) 
  cat('---------------------------------------------','\n') 
  cat('Kesimpulan hasil uji parsial','\n') 
  cat('---------------------------------------------','\n') 
  for (i in 1:n1) 
  { 
    thit[i]=B[i,1]/SE[i] 
    pval[i]=2*(pt(abs(thit[i]),(n-n1),lower.tail=FALSE)) 
    if (pval[i]<=alpha) cat('Tolak Ho yakni variabel bebas signifikan dengan pvalue',pval[i],'\n') 
    else cat('Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue',pval[i],'\n') 
  } 
  thit=as.matrix(thit) 
  cat('=============================================','\n') 
  cat('nilai t hitung','\n') 
  cat('=============================================','\n') 
  print(thit) 
  cat('Analysis of Variance','\n') 
  cat('=============================================','\n') 
  cat('Sumber df SS MS Fhit','\n') 
  cat('Regresi ',(n1-1),' ',SSR,' ',MSR,' ',Fhit,'\n') 
  cat('Error ',n-n1,' ',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='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI RESIDUAL knot2.csv')
  write.csv(mx,file='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI MX knot2.csv') 
  write.csv(yhat,file='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI yhat knot2.csv') 
} 
uji(0.1, 0)
## --------------------------------------- 
## Kesimpulan hasil uji simultan 
## --------------------------------------- 
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan 
##  
## --------------------------------------------- 
## Kesimpulan hasil uji parsial 
## --------------------------------------------- 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.2064143 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.002583234 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.004789845 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.2422198 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.1561285 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.2229634 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 1.532736e-05 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.1557346 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.1498621 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.1649558 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.02812549 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.02406025 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 3.439945e-06 
## ============================================= 
## nilai t hitung 
## ============================================= 
##            [,1]
##  [1,]  1.282856
##  [2,] -3.199948
##  [3,]  2.975107
##  [4,]  1.185782
##  [5,]  1.443491
##  [6,] -1.236557
##  [7,] -4.872993
##  [8,]  1.444899
##  [9,] -1.466226
## [10,]  1.412651
## [11,]  2.272317
## [12,] -2.338846
## [13,]  5.327626
## Analysis of Variance 
## ============================================= 
## Sumber df SS MS Fhit 
## Regresi  12   183.7252   15.31043   19.73097 
## Error  43   43.45373   0.7759595 
## Total  55   227.1789 
## ============================================= 
## s= 0.8808856  Rsq= 80.87246 
## pvalue(F)= 1.075619e-13
#PEMILIHAN TIGA TITIK KNOT# 
GCV3=function(data) 
{ 
  para=0 
  data=as.matrix(data) 
  N=nrow(data) 
  M=ncol(data) 
  m=ncol(data)-para-1 
  dataA=data[,(para+2):M] 
  dataA=as.matrix(dataA) 
  F=diag(N) 
  nk=50 
  knot1=matrix(ncol=m,nrow=nk) 
  for (i in (1:m)) 
  { 
    a=seq(min(dataA[,i]),max(dataA[,i]),length.out=nk) 
    knot1[,i]=t(as.matrix(a)) 
  } 
  a1=length(knot1[,1]) 
  knot1=as.matrix(knot1[2:(a1-1),]) 
  a2=nk-2 
  z=(a2*(a2-1)*(a2-2)/6) 
  knot2=cbind(rep(NA,(z+1))) 
  for (i in (1:m)) 
  { 
    knot=rbind(rep(NA,3)) 
    for ( j in 1:(a2-2)) 
    { 
      for (k in (j+1):(a2-1)) 
      { 
        for (g in (k+1):a2) 
        { 
          xx=cbind(knot1[j,i],knot1[k,i],knot1[g,i]) 
          knot=rbind(knot,xx) 
        } 
      } 
    } 
    knot2=cbind(knot2,knot) 
  } 
  knot2=knot2[2:(z+1),2:(3*m+1)] 
  a3=nrow(knot2) 
  aa=rep(1,N) 
  data1=matrix(ncol=3*m,nrow=N) 
  data2=data[,2:M] 
  nk1=nrow(knot2) 
  GCV=as.matrix(rep(NA,nk1),ncol=1);colnames(GCV)<-"GCV" 
  MSE=as.matrix(rep(NA,nk1),ncol=1);colnames(MSE)<-"MSE" 
  SSE=rep(NA,nk1) 
  SSR=rep(NA,nk1) 
  Rsq=as.matrix(rep(NA,nk1),ncol=1);colnames(Rsq)<-"Rsq" 
  knotke=matrix(c(1:nk1),ncol=1);colnames(knotke)<-"knot_ke" 
  for (i in 1:a3) 
  { 
    for (j in 1:(3*m)) 
    { 
      b=ceiling(j/3) 
      for (k in 1:N) 
      { 
        if (data[k,(b+para+1)]<knot2[i,j]) data1[k,j]=0 else  
          data1[k,j]=data[k,(b+para+1)]-knot2[i,j] 
      } 
    } 
    mx=as.matrix(cbind(aa,data2,data1)) 
    C=pinv(t(mx)%*%mx) 
    B=C%*%(t(mx)%*%data[,1]) 
    yhat=mx%*%B 
    res=data[,1]-yhat 
    SSE[i]=sum((res)^2) 
    SSR[i]=sum((yhat-mean(data[,1]))^2) 
    MSE[i]=SSE[i]/(N) 
    Rsq[i]=(SSR[i]/(SSR[i]+SSE[i]))*100 
    A=mx%*%C%*%t(mx) 
    A1=(F-A) 
    A2=(sum(diag(A1))/N)^2 
    GCV[i]=MSE[i]/A2 
  } 
  dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot2)) 
  dataG=dataAll[order(GCV),-2] 
  write.csv(dataAll,file="C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\DATA ALL knot 3.csv") 
  cat("==============================================","\n") 
  cat("HASIL GCV terkecil dengan 3 knot","\n") 
  cat("==============================================","\n") 
  print(((dataG[1,1:14]))) 
  cat("Nilai GCV 10 terkecil pertama","\n") 
  print(dataG[1:10,]) 
  
  # --- PERBAIKAN: hitung ulang B dari titik knot dengan GCV TERKECIL --- 
  baris_terbaik = dataG[1, "knot_ke"] 
  knotgcv = knot2[baris_terbaik, ] 
  datagcv1=matrix(ncol=3*m,nrow=N) 
  for (j in 1:(3*m)) 
  { 
    b=ceiling(j/3) 
    for (k in 1:N) 
    { 
      if (data[k,(b+para+1)]<knotgcv[j]) datagcv1[k,j]=0 else 
        datagcv1[k,j]=data[k,(b+para+1)]-knotgcv[j] 
    } 
  } 
  mxgcv=as.matrix(cbind(aa,data2,datagcv1)) 
  Cgcv=pinv(t(mxgcv)%*%mxgcv) 
  B=Cgcv%*%(t(mxgcv)%*%data[,1]) 
  
  cat("\n")  
  cat("==============================================","\n") 
  cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 3","\n") 
  cat("==============================================","\n") 
  print(B) 
  cat("\n") 
} 
GCV3(data)
## ============================================== 
## HASIL GCV terkecil dengan 3 knot 
## ============================================== 
##         GCV     knot_ke                                                 
##    1.444246 2200.000000   61.132449   61.854082   76.286735   11.382041 
##                                                                         
##   11.536735   14.630612    7.692245   10.620408   69.183673   25.025510 
##                         
##   27.449184   75.922653 
## Nilai GCV 10 terkecil pertama 
##            GCV knot_ke                                                      
##  [1,] 1.444246    2200 61.13245 61.85408 76.28673 11.38204 11.53673 14.63061
##  [2,] 1.445793    3146 61.49327 61.85408 76.28673 11.45939 11.53673 14.63061
##  [3,] 1.447104    2157 61.13245 61.49327 76.28673 11.38204 11.45939 14.63061
##  [4,] 1.452451    3106 61.13245 77.00837 77.36918 11.38204 14.78531 14.86265
##  [5,] 1.453903    2199 61.13245 61.85408 75.92592 11.38204 11.53673 14.55327
##  [6,] 1.454858    3145 61.49327 61.85408 75.92592 11.45939 11.53673 14.55327
##  [7,] 1.458180    2156 61.13245 61.49327 75.92592 11.38204 11.45939 14.55327
##  [8,] 1.458931    2242 61.13245 62.21490 76.28673 11.38204 11.61408 14.63061
##  [9,] 1.464584    3037 61.13245 73.03939 76.28673 11.38204 13.93449 14.63061
## [10,] 1.465580    3101 61.13245 76.28673 76.64755 11.38204 14.63061 14.70796
##                                                             
##  [1,] 7.692245 10.620408 69.18367 25.02551 27.44918 75.92265
##  [2,] 9.156327 10.620408 69.18367 26.23735 27.44918 75.92265
##  [3,] 7.692245  9.156327 69.18367 25.02551 26.23735 75.92265
##  [4,] 7.692245 72.111837 73.57592 25.02551 78.34633 79.55816
##  [5,] 7.692245 10.620408 67.71959 25.02551 27.44918 74.71082
##  [6,] 9.156327 10.620408 67.71959 26.23735 27.44918 74.71082
##  [7,] 7.692245  9.156327 67.71959 25.02551 26.23735 74.71082
##  [8,] 7.692245 12.084490 69.18367 25.02551 28.66102 75.92265
##  [9,] 7.692245 56.006939 69.18367 25.02551 65.01612 75.92265
## [10,] 7.692245 69.183673 70.64776 25.02551 75.92265 77.13449
## 
## ============================================== 
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 3 
## ============================================== 
##               [,1]
##  [1,]  -4.87954772
##  [2,]  -0.57070096
##  [3,]   1.29861827
##  [4,]   0.20162854
##  [5,]   1.23897094
##  [6,]  -2.60712939
##  [7,]   3.03208076
##  [8,]   0.76491110
##  [9,]  13.80910095
## [10,] -14.12188194
## [11,] -19.59971920
## [12,]  -0.22960609
## [13,]   0.02361872
## [14,]   0.25980650
## [15,]  -1.94377923
## [16,]   0.67438516
## [17,]   0.92445503

#UJI SIGNIFIKANSI TIGA TITIK KNOT#

uji=function(alpha,para) 
{ 
  knot = read.csv("C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\DATA ALL knot 3.csv", header=TRUE)
  knot = as.matrix(knot) 
  knot = knot[,-1] 
  knot = knot[order(knot[,1]), ] 
  baris_knot = as.numeric(knot[1, 4:15])   # 12 nilai knot terbaik (3*m = 12) 
  
  data=as.matrix(data) 
  ybar=mean(data[,1]) 
  m=para+2 
  n=nrow(data) 
  q=ncol(data) 
  dataA=cbind(data[,m],data[,m],data[,m], 
              data[,m+1],data[,m+1],data[,m+1], 
              data[,m+2],data[,m+2],data[,m+2], 
              data[,m+3],data[,m+3],data[,m+3]) 
  dataA=as.matrix(dataA) 
  satu=rep(1,n) 
  n1=length(baris_knot)              # = 12 
  data.knot=matrix(ncol=n1,nrow=n) 
  for (i in 1:n1) 
  { 
    for (j in 1:n) 
    { 
      if(dataA[j,i]<baris_knot[i]) data.knot[j,i]=0 
      else data.knot[j,i]=dataA[j,i]-baris_knot[i] 
    } 
  } 
  mx=cbind(satu,data[,2],data.knot[,1:3],data[,3], 
           data.knot[,4:6],data[,4],data.knot[,7:9], 
           data[,5],data.knot[,10:12]) 
  mx=as.matrix(mx) 
  B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1] 
  n1=nrow(B) 
  yhat=mx%*%B 
  ybar=mean(data[,1]) 
  res=data[,1]-yhat 
  SSE=sum((data[,1]-yhat)^2) 
  SSR=sum((yhat-ybar)^2) 
  MSE=SSE/(n) 
  MSR=SSR/(n1-1) 
  SST=sum((data[,1]-ybar)^2) 
  Rsq=(SSR/(SSR+SSE))*100
  Fhit=MSR/MSE 
  pvalue=pf(Fhit,(n1-1),(n-n1),lower.tail=FALSE) 
  if(pvalue<=alpha) 
  { 
    cat('---------------------------------------','\n') 
    cat('Kesimpulan hasil uji simultan','\n') 
    cat('---------------------------------------','\n') 
    cat('Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan','\n') 
    cat('','\n') 
  } 
  else 
  { 
    cat('---------------------------------------','\n') 
    cat('Kesimpulan hasil uji simultan','\n') 
    cat('---------------------------------------','\n') 
    cat('Gagal Tolak Ho yakni semua variabel bebas tidak berpengaruh signifikan','\n') 
    cat('','\n') 
  } 
  thit=rep(NA,n1) 
  pval=rep(NA,n1) 
  SE=sqrt(diag(MSE*(pinv(t(mx)%*%mx)))) 
  cat('---------------------------------------------','\n') 
  cat('Kesimpulan hasil uji parsial','\n') 
  cat('---------------------------------------------','\n') 
  for (i in 1:n1) 
  { 
    thit[i]=B[i,1]/SE[i] 
    pval[i]=2*(pt(abs(thit[i]),(n-n1),lower.tail=FALSE)) 
    if (pval[i]<=alpha) cat('Tolak Ho yakni variabel bebas signifikan dengan pvalue',pval[i],'\n') 
    else cat('Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue',pval[i],'\n') 
  } 
  thit=as.matrix(thit) 
  cat('=============================================','\n') 
  cat('nilai t hitung','\n') 
  cat('=============================================','\n') 
  print(thit) 
  cat('Analysis of Variance','\n') 
  cat('=============================================','\n') 
  cat('Sumber df SS MS Fhit','\n') 
  cat('Regresi ',(n1-1),' ',SSR,' ',MSR,' ',Fhit,'\n') 
  cat('Error ',n-n1,' ',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='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI RESIDUAL knot3.csv')
  write.csv(mx,file='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI mx knot3.csv') 
  write.csv(yhat,file='C:\\File Lama\\BACKUP\\Kodingan dan segalanya\\HTML R\\OUTPUT UJI yhat knot3.csv') 
} 
uji(0.1, 0)
## --------------------------------------- 
## Kesimpulan hasil uji simultan 
## --------------------------------------- 
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang signifikan 
##  
## --------------------------------------------- 
## Kesimpulan hasil uji parsial 
## --------------------------------------------- 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.6087605 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.7734961 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.5228609 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.1758342 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.3014451 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.8976796 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.6190703 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.4479088 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 1.040868e-05 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.4567152 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.6607985 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.9315849 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.2034532 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 0.04980339 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.1138228 
## Gagal tolak Ho yakni variabel tidak signifikan dengan pvalue 0.2996611 
## Tolak Ho yakni variabel bebas signifikan dengan pvalue 2.587599e-06 
## ============================================= 
## nilai t hitung 
## ============================================= 
##              [,1]
##  [1,] -0.51601043
##  [2,] -0.28981099
##  [3,] -0.64475471
##  [4,]  1.37873375
##  [5,]  1.04722022
##  [6,]  0.12943381
##  [7,]  0.50116758
##  [8,] -0.76664478
##  [9,] -5.05857856
## [10,]  0.75174871
## [11,] -0.44218412
## [12,]  0.08640707
## [13,]  1.29348236
## [14,]  2.02451772
## [15,] -1.61753944
## [16,]  1.05114669
## [17,]  5.49666947
## Analysis of Variance 
## ============================================= 
## Sumber df SS MS Fhit 
## Regresi  16   185.9148   11.61967   15.76917 
## Error  39   41.26416   0.7368601 
## Total  55   227.1789 
## ============================================= 
## s= 0.8584055  Rsq= 81.83627 
## pvalue(F)= 2.762142e-12