#=========================================================#
#ANALISIS STATISTIKA DESKRIPTIF, DETEKSI MULTIKOLINIERITAS, &
#MODEL REGRESI NONPARAMTERIK SPLINE TRUNCATED #=========================================================#
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)
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
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))
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
#scatterplot#
library(ggplot2)
library(gridExtra)
options(scipen = 999)

buat_gg <- function(xvar, xlab) {
  xlab <- paste(strwrap(xlab, width = 28), collapse = "\n")   # label x dipotong jadi 2 baris
  ggplot(data, aes(x = .data[[xvar]], y = Y)) +
    geom_point(shape = 21, color = "steelblue4", fill = "steelblue",
               size = 2.5, alpha = 0.7) +
    geom_smooth(method = "loess", formula = y ~ x,
                color = "firebrick", fill = "pink", alpha = 0.4, se = TRUE) +
    labs(x = xlab, y = "Tingkat Pengangguran\nTerbuka") +
    theme_bw(base_size = 10) +
    theme(axis.title = element_text(size = 9),
          plot.margin = margin(8, 10, 8, 8))
}

sp1 <- buat_gg("X1", "Tingkat Partisipasi Angkatan Kerja (X1)")
sp2 <- buat_gg("X2", "Harapan Lama Sekolah (X2)")
sp3 <- buat_gg("X3", "PDRB Atas Dasar Harga Berlaku (X3)")
sp4 <- buat_gg("X4", "Indeks Pembangunan Manusia (X4)")

grid.arrange(sp1, sp2, sp3, sp4, ncol = 2)

#scatterplot#
library(ggplot2)
options(scipen = 999)

buat_gg <- function(xvar, xlab) {
  ggplot(data, aes(x = .data[[xvar]], y = Y)) +
    geom_point(shape = 21, color = "steelblue4", fill = "steelblue",
               size = 3, alpha = 0.7) +
    geom_smooth(method = "loess", formula = y ~ x,
                color = "firebrick", fill = "pink", alpha = 0.4, se = TRUE) +
    labs(x = xlab, y = "Tingkat Pengangguran Terbuka") +
    theme_bw(base_size = 12)
}

# Scatterplot X1
buat_gg("X1", "Tingkat Partisipasi Angkatan Kerja (X1)")

# Scatterplot X2
buat_gg("X2", "Harapan Lama Sekolah (X2)")

# Scatterplot X3
buat_gg("X3", "PDRB Atas Dasar Harga Berlaku (X3)")

# Scatterplot X4
buat_gg("X4", "Indeks Pembangunan Manusia (X4)")

#======================#
#PEMILIHAN 1 TITIK KNOT#
#======================#
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
GCV1=function(data)
{
  library(Matrix)
  library(pracma)
  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="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/dataAll 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,2],])
  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))
  C=pinv(t(mxgcv)%*%mxgcv)
  B=C%*%(t(mxgcv)%*%data[,1])
  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)
## 
## Attaching package: 'pracma'
## The following objects are masked from 'package:Matrix':
## 
##     expm, lu, tril, triu
## The following object is masked from 'package:car':
## 
##     logit
## ============================================== 
## 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,]  10.122393144
##  [2,]  -0.200697922
##  [3,]   0.762874047
##  [4,]  -0.005244856
##  [5,]  -0.017306261
##  [6,]   0.072850412
##  [7,] -21.843749933
##  [8,]   0.415138942
##  [9,]   1.053106627
#=========================#
#satu knot Uji signifikansi#
#=========================#
library(pracma)
uji=function(alpha,para)
{
  data=read.table("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/data regresi spline truncated.txt",header=TRUE)
  knot_raw <- read.csv("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/dataAll knot 1.csv", header = FALSE, skip = 1)
  colnames(knot_raw) <- c("no","GCV","Rsq","knot_ke","k1","k2","k3","k4")
  knotgcv <- knot_raw[which.min(knot_raw$GCV), c("k1","k2","k3","k4")]
  knot <- matrix(as.numeric(knotgcv), nrow = 1)
  data=as.matrix(data)
  knot=as.matrix(knot)
  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=ncol(knot)
  data.knot=matrix(ncol=n1,nrow=n)
  for (i in 1:n1)
  {
    for (j in 1:n)
    {
      if(dataA[j,i]<knot[1,i])data.knot[j,i]=0
      else data.knot[j,i]=dataA[j,i]-knot[1,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-n1)
  MSR=SSR/(n1-1)
  SST=sum((data[,1]-ybar)^2)
  Rsq=(SSR/(SSR+SSE))*100
  #-------------------------------------------------------#
  #SINTAKS UJI SIMULTAN DENGAN TITIK KNOT #
  #-------------------------------------------------------#
  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')
  }
  #------------------------------------------------------#
  #SINTAKS UJI PARSIAL DENGAN TITIK KNOT #
  #------------------------------------------------------#
  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')
  thit=rep(NA,n1)
  pval=rep(NA,n1)
  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="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji residual knot1.csv")
  write.csv(mx,file="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji mx knot1.csv")
  write.csv(yhat,file="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji yhat knot1.csv")
}
uji(0.05,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.06500353 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.00002835411 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.9156743 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.01250812 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.00001020449 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.541134 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.09644074 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.1203336 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.000002335269 
## ============================================= 
## nilai t hitung 
## ============================================= 
##             [,1]
##  [1,]  1.8894905
##  [2,] -4.6371698
##  [3,]  0.1064552
##  [4,]  2.5972210
##  [5,] -4.9427992
##  [6,] -0.6155913
##  [7,]  1.6963165
##  [8,] -1.5820999
##  [9,]  5.3759508
## 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)= 0.00000000002451288
#DUA TITIK KNOT#
data=read.table("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/data regresi spline truncated.txt",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
GCV2=function(data)
{
  library(Matrix)
  library(pracma)
  para=0
  data=as.matrix(data)
  N=nrow(data)
  M=ncol(data)
  m=ncol(data)-para-1 #m = banyaknya var non parametrik
  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),])
  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] # data x saja
  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="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/dataAll 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,])
  mingcv=dataG[1,1]
  knotgcv=as.matrix(knot2[dataG[1,2],])
  knotgcv1=matrix(knotgcv,nrow=2)
  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,1]) datagcv1[k,j]=0 else
        datagcv1[k,j]=data[k,(b+para+1)]-knotgcv[j,1]
    }
  }
  mxgcv=as.matrix(cbind(aa,data2,datagcv1))
  C=pinv(t(mxgcv)%*%mxgcv)
  B=C%*%(t(mxgcv)%*%data[,1])
  mxgcv=mxgcv[,c(2:6)]
  list(knotgcv=knotgcv1,mingcv=mingcv,mxgcv=mxgcv)
  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.3587176
##  [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
#=========================#
#dua knot Uji signifikansi#
#=========================#
library(pracma)
uji=function(alpha,para)
{
  data=read.table("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/data regresi spline truncated.txt",header=TRUE)
  knot_raw <- read.csv("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/dataAll knot 2.csv", header = FALSE, skip = 1)
  colnames(knot_raw) <- c("no","GCV","Rsq","knot_ke",
                          "k1a","k1b","k2a","k2b","k3a","k3b","k4a","k4b")
  knotgcv <- knot_raw[which.min(knot_raw$GCV),
                      c("k1a","k1b","k2a","k2b","k3a","k3b","k4a","k4b")]
  knot <- matrix(as.numeric(knotgcv), nrow = 1) 
  data=as.matrix(data)
  knot=as.matrix(knot)
  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=ncol(knot)
  data.knot=matrix(ncol=n1,nrow=n)
  for (i in 1:n1)
  {
    for (j in 1:n)
    {
      if(dataA[j,i]<knot[1,i])data.knot[j,i]=0
      else data.knot[j,i]=dataA[j,i]-knot[1,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-n1)
  MSR=SSR/(n1-1)
  SST=sum((data[,1]-ybar)^2)
  Rsq=(SSR/(SSR+SSE))*100
  #-------------------------------------------------------#
  #SINTAKS UJI SIMULTAN DENGAN TITIK KNOT #
  #-------------------------------------------------------#
  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')
  }
  #------------------------------------------------------#
  #SINTAKS UJI PARSIAL DENGAN TITIK KNOT #
  #------------------------------------------------------#
  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')
  thit=rep(NA,n1)
  pval=rep(NA,n1)
  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="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji residual knot2.csv")
  write.csv(mx,file="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji mx knot2.csv")
  write.csv(yhat,file="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji yhat knot2.csv")
}
uji(0.05,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.2671925 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.007544198 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.01250391 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.3045796 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2127214 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2845971 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.0001058266 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2122837 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2057346 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2224808 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.05283963 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.0465475 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.0000297441 
## ============================================= 
## nilai t hitung 
## ============================================= 
##            [,1]
##  [1,]  1.124135
##  [2,] -2.804034
##  [3,]  2.607011
##  [4,]  1.039070
##  [5,]  1.264895
##  [6,] -1.083563
##  [7,] -4.270080
##  [8,]  1.266128
##  [9,] -1.284816
## [10,]  1.237870
## [11,]  1.991173
## [12,] -2.049471
## [13,]  4.668464
## 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)= 0.000000000009568368
#TIGA TITIK KNOT#
data=read.table("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/data regresi spline truncated.txt",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
GCV3=function(data)
{
  library(Matrix)
  library(pracma)
  para=0
  data=as.matrix(data)
  N=nrow(data)
  M=ncol(data)
  m=ncol(data)-para-1 #m = banyaknya var non parametrik
  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),])
  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] # data x saja
  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="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/dataAll 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,])
  mingcv=dataG[1,1]
  knotgcv=as.matrix(knot2[dataG[1,2],])
  knotgcv1=matrix(knotgcv,nrow=3)
  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,1]) datagcv1[k,j]=0 else
        datagcv1[k,j]=data[k,(b+para+1)]-knotgcv[j,1]
    }
  }
  mxgcv=as.matrix(cbind(aa,data2,datagcv1))
  C=pinv(t(mxgcv)%*%mxgcv)
  B=C%*%(t(mxgcv)%*%data[,1])
  mxgcv=mxgcv[,c(2:6)]
  list(knotgcv=knotgcv1,mingcv=mingcv,mxgcv=mxgcv)
  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.87954771
##  [2,]  -0.57070096
##  [3,]   1.29861828
##  [4,]   0.20162854
##  [5,]   1.23897094
##  [6,]  -2.60712938
##  [7,]   3.03208076
##  [8,]   0.76491110
##  [9,]  13.80910092
## [10,] -14.12188192
## [11,] -19.59971920
## [12,]  -0.22960609
## [13,]   0.02361872
## [14,]   0.25980650
## [15,]  -1.94377923
## [16,]   0.67438515
## [17,]   0.92445503
#=========================#
#tiga knot Uji signifikansi#
#=========================#
library(pracma)
uji=function(alpha,para)
{
  data=read.table("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/data regresi spline truncated.txt",header=TRUE)
  
  knot_raw <- read.csv("D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/dataAll knot 3.csv", header = FALSE, skip = 1)
  colnames(knot_raw) <- c("no","GCV","Rsq","knot_ke",
                          "k1a","k1b","k1c","k2a","k2b","k2c",
                          "k3a","k3b","k3c","k4a","k4b","k4c")
  knotgcv <- knot_raw[which.min(knot_raw$GCV),
                      c("k1a","k1b","k1c","k2a","k2b","k2c",
                        "k3a","k3b","k3c","k4a","k4b","k4c")]
  knot <- matrix(as.numeric(knotgcv), nrow = 1)
  data=as.matrix(data)
  knot=as.matrix(knot)
  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=ncol(knot)
  data.knot=matrix(ncol=n1,nrow=n)
  for (i in 1:n1)
  {
    for (j in 1:n)
    {
      if(dataA[j,i]<knot[1,i])data.knot[j,i]=0
      else data.knot[j,i]=dataA[j,i]-knot[1,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-n1)
  MSR=SSR/(n1-1)
  SST=sum((data[,1]-ybar)^2)
  Rsq=(SSR/(SSR+SSE))*100
  #-------------------------------------------------------#
  #SINTAKS UJI SIMULTAN DENGAN TITIK KNOT #
  #-------------------------------------------------------#
  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')
  }
  
  #------------------------------------------------------#
  #SINTAKS UJI PARSIAL DENGAN TITIK KNOT #
  #------------------------------------------------------#
  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')
  thit=rep(NA,n1)
  pval=rep(NA,n1)
  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="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji residual knot3.csv")
  write.csv(mx,file="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji mx knot3.csv")
  write.csv(yhat,file="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 5_Regresi Spline Truncated (DATA SKRIPSI)/output uji yhat knot3.csv")
}
uji(0.05,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.6691134 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.8101606 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.5935914 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2569091 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.3875066 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.9145372 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.6780688 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.5260547 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.0001404269 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.5340848 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.7141141 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.942884 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.2870221 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.09910266 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.1848397 
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.3857466 
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.00004564589 
## ============================================= 
## nilai t hitung 
## ============================================= 
##              [,1]
##  [1,] -0.43062256
##  [2,] -0.24185393
##  [3,] -0.53806261
##  [4,]  1.15058497
##  [5,]  0.87392932
##  [6,]  0.10801549
##  [7,]  0.41823585
##  [8,] -0.63978267
##  [9,] -4.22149996
## [10,]  0.62735156
## [11,] -0.36901280
## [12,]  0.07210868
## [13,]  1.07944073
## [14,]  1.68950652
## [15,] -1.34987380
## [16,]  0.87720604
## [17,]  4.58709688
## 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)= 0.0000000007291165