LIBRARY
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=T)
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
1. STATISTIKA DESKRIPTIF
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
2. DETEKSI MULTIKOLINIERITAS
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
3. SCATTERPLOT
library(ggplot2)
library(gridExtra)
p1 = ggplot(data, aes(x = X1, y = Y)) +
geom_point(color = "steelblue4", size = 2, alpha = 0.7) +
geom_smooth(color = "firebrick", fill = "lightpink", se = TRUE) +
labs(
x = "Tingkat Partisipasi Angkatan Kerja (X1)",
y = "Tingkat Pengangguran Terbuka (Y)"
) +
theme_bw()
p2 = ggplot(data, aes(x = X2, y = Y)) +
geom_point(color = "steelblue4", size = 2, alpha = 0.7) +
geom_smooth(color = "firebrick", fill = "lightpink", se = TRUE) +
labs(
x = "Harapan Lama Sekolah (X2)",
y = "Tingkat Pengangguran Terbuka (Y)"
) +
theme_bw()
p3 = ggplot(data, aes(x = X3, y = Y)) +
geom_point(color = "steelblue4", size = 2, alpha = 0.7) +
geom_smooth(color = "firebrick", fill = "lightpink", se = TRUE) +
labs(
x = "PDRB Atas Dasar Harga Berlaku (X3)",
y = "Tingkat Pengangguran Terbuka (Y)"
) +
theme_bw()
p4 = ggplot(data, aes(x = X4, y = Y)) +
geom_point(color = "steelblue4", size = 2, alpha = 0.7) +
geom_smooth(color = "firebrick", fill = "lightpink", se = TRUE) +
labs(
x = "Indeks Pembangunan Manusia (X4)",
y = "Tingkat Pengangguran Terbuka (Y)"
) +
theme_bw()
grid.arrange(p1, p2, p3, p4, ncol = 2)
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
4. PEMILIHAN 1 TITIK KNOT
data=read.table("C:/Users/acer/Downloads/Regnon/data skripsi GWS tanpa kolom kabupaten kota.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
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="C:/Users/acer/Downloads/Regnon/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,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")
idx_opt=dataG[1,2] # nomor baris knot dengan GCV terkecil
knot_opt=as.matrix(knot1[idx_opt,]) # 4 knot optimal (x1..x4)
data_opt=matrix(ncol=m,nrow=N)
for (j in 1:m)
{
for (k in 1:N)
{
if (data[k,(j+para+1)]<knot_opt[j,1]) data_opt[k,j]=0 else
data_opt[k,j]=data[k,(j+para+1)]-knot_opt[j,1]
}
}
mx_opt=as.matrix(cbind(aa,data2,data_opt))
B_opt=pinv(t(mx_opt)%*%mx_opt)%*%(t(mx_opt)%*%data[,1])
rownames(B_opt)=c("b0","x1","x2","x3","x4",
"(x1-k1)+","(x2-k1)+","(x3-k1)+","(x4-k1)+")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER 1 TITIK KNOT (OPTIMAL)","\n")
cat("==============================================","\n")
print(round(B_opt,7))
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,] 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
##
## ==============================================
## HASIL ESTIMASI PARAMETER 1 TITIK KNOT (OPTIMAL)
## ==============================================
## [,1]
## b0 10.1223931
## x1 -0.2006979
## x2 0.7628740
## x3 -0.0052449
## x4 -0.0173063
## (x1-k1)+ 0.0728504
## (x2-k1)+ -21.8437499
## (x3-k1)+ 0.4151389
## (x4-k1)+ 1.0531066
5. SATU KNOT UJI SIGNIFIKANSI
library(pracma)
uji=function(alpha,para)
{
alpha=0.1
para=0
data=read.table("C:/Users/acer/Downloads/Regnon/data skripsi GWS tanpa kolom kabupaten kota.txt",header=TRUE)
knotall=read.csv("C:/Users/acer/Downloads/Regnon/dataAll knot 1.csv",header=TRUE)
# Urutkan berdasarkan GCV terkecil
knotall=knotall[order(knotall$GCV),]
# Ambil 4 titik knot terbaik
knot=as.matrix(knotall[1,5:8])
data=as.matrix(data)
print(knot)
print(dim(knot))
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=4
data.knot=matrix(0,nrow=n,ncol=n1)
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="C:/Users/acer/Downloads/Regnon/output uji residual knot1.csv")
write.csv(mx,file="C:/Users/acer/Downloads/Regnon/output uji mx knot1.csv")
write.csv(yhat,file="C:/Users/acer/Downloads/Regnon/output uji yhat knot1.csv")
}
uji(alpha,0)
## X.1 X.2 X.3 X.4
## 45 76.28673 14.63061 69.18367 75.92265
## [1] 1 4
## ---------------------------------------
## 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.06500353
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 2.835411e-05
## 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 1.020449e-05
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.541134
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.09644074
## Gagal tolak Ho yakni variabel tidak
## signifikan dengan pvalue 0.1203336
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 2.335269e-06
## =============================================
## 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)= 2.451288e-11
6. DUA TITIK KNOT
data=read.table("C:/Users/acer/Downloads/Regnon/data skripsi GWS tanpa kolom kabupaten kota.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="C:/Users/acer/Downloads/Regnon/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,4],])
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))
mxgcv=mxgcv[,c(2:6)]
list(knotgcv=knotgcv1,mingcv=mingcv,mcgcv=mxgcv)
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 2","\n")
cat("==============================================","\n")
print(B)
cat("\n")
idx_opt=dataG[1,2]
knot_opt=as.matrix(knot2[idx_opt,])
data_opt=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)]<knot_opt[j,1]) data_opt[k,j]=0 else
data_opt[k,j]=data[k,(b+para+1)]-knot_opt[j,1]
}
}
mx_opt=as.matrix(cbind(aa,data2,data_opt))
B_opt=pinv(t(mx_opt)%*%mx_opt)%*%(t(mx_opt)%*%data[,1])
B_opt=B_opt[c(1:5, 6,8,10,12, 7,9,11,13),,drop=FALSE]
rownames(B_opt)=c("b0","x1","x2","x3","x4",
"(x1-k1)+","(x2-k1)+","(x3-k1)+","(x4-k1)+",
"(x1-k2)+","(x2-k2)+","(x3-k2)+","(x4-k2)+")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER 2 TITIK KNOT (OPTIMAL)","\n")
cat("==============================================","\n")
print(round(B_opt,7))
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,] 6.044955e+00
## [2,] -1.859007e-01
## [3,] 1.031486e+00
## [4,] -4.170952e-03
## [5,] -2.228411e-02
## [6,] 9.154406e+00
## [7,] -2.164657e+01
## [8,] -3.547727e+02
## [9,] 6.943980e+02
## [10,] 4.083286e-02
## [11,] 2.041643e-02
## [12,] 4.845712e+00
## [13,] -1.020433e+01
##
## ==============================================
## HASIL ESTIMASI PARAMETER 2 TITIK KNOT (OPTIMAL)
## ==============================================
## [,1]
## b0 95.3587176
## x1 -2.9213119
## x2 6.1469387
## x3 0.2059982
## x4 0.7283636
## (x1-k1)+ 2.7449916
## (x2-k1)+ -5.3197766
## (x3-k1)+ -0.2144538
## (x4-k1)+ -0.7619877
## (x1-k2)+ 0.7954129
## (x2-k2)+ -18.0705567
## (x3-k2)+ 0.2853691
## (x4-k2)+ 0.9034632
7. DUA KNOT UJI SIGNIFIKANSI
library(pracma)
uji=function(alpha,para)
{
alpha=0.1
para=0
data=read.table("C:/Users/acer/Downloads/Regnon/data skripsi GWS tanpa kolom kabupaten kota.txt",header=TRUE)
dataAll=read.csv("C:/Users/acer/Downloads/Regnon/dataAll knot 2.csv",header=TRUE)
# Urutkan berdasarkan GCV terkecil
dataAll=dataAll[order(dataAll$GCV),]
# Ambil 8 nilai knot dari baris dengan GCV terkecil
knot=as.numeric(dataAll[1,5:12])
# Ubah menjadi matriks 1 x 8
knot=matrix(knot,nrow=1)
print(knot)
print(dim(knot))
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=8
data.knot=matrix(0,nrow=n,ncol=n1)
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="C:/Users/acer/Downloads/Regnon/output uji residual knot2.csv")
write.csv(mx,file="C:/Users/acer/Downloads/Regnon/output uji mx knot2.csv")
write.csv(yhat,file="C:/Users/acer/Downloads/Regnon/output uji yhat knot2.csv")
}
uji(alpha,0)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 61.13245 76.28673 11.38204 14.63061 7.692245 69.18367 25.02551 75.92265
## [1] 1 8
## ---------------------------------------
## 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
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.05283963
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 0.0465475
## Tolak Ho yakni variabel bebas
## signifikan dengan pvalue 2.97441e-05
## =============================================
## 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)= 9.568368e-12
8. TIGA TITIK KNOT
data=read.table("C:/Users/acer/Downloads/Regnon/data skripsi GWS tanpa kolom kabupaten kota.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="C:/Users/acer/Downloads/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")
mingcv=dataG[1,1]
knotgcv=as.matrix(knot2[dataG[1,4],])
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))
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")
idx_opt=dataG[1,2] # nomor baris knot dengan GCV terkecil
knot_opt=as.matrix(knot2[idx_opt,]) # 12 knot: x1k1,x1k2,x1k3,x2k1,...
data_opt=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)]<knot_opt[j,1]) data_opt[k,j]=0 else
data_opt[k,j]=data[k,(b+para+1)]-knot_opt[j,1]
}
}
mx_opt=as.matrix(cbind(aa,data2,data_opt))
B_opt=pinv(t(mx_opt)%*%mx_opt)%*%(t(mx_opt)%*%data[,1])
# urutkan: knot ke-1 semua x, lalu knot ke-2, lalu knot ke-3
B_opt=B_opt[c(1:5, 6,9,12,15, 7,10,13,16, 8,11,14,17),,drop=FALSE]
rownames(B_opt)=c("b0","x1","x2","x3","x4",
"(x1-k1)+","(x2-k1)+","(x3-k1)+","(x4-k1)+",
"(x1-k2)+","(x2-k2)+","(x3-k2)+","(x4-k2)+",
"(x1-k3)+","(x2-k3)+","(x3-k3)+","(x4-k3)+")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER 3 TITIK KNOT (OPTIMAL)","\n")
cat("==============================================","\n")
print(round(B_opt,7))
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
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 3
## ==============================================
## [,1]
## [1,] 6.072786923
## [2,] -0.185927169
## [3,] 1.029244397
## [4,] -0.004126666
## [5,] -0.022270952
## [6,] 6.655420960
## [7,] -4.161967634
## [8,] -14.979356227
## [9,] -6.028540779
## [10,] 14.203131537
## [11,] 9.970360436
## [12,] 0.022331619
## [13,] 0.014887746
## [14,] 0.007443873
## [15,] -12.848704041
## [16,] 31.920041594
## [17,] -28.068467337
##
## ==============================================
## HASIL ESTIMASI PARAMETER 3 TITIK KNOT (OPTIMAL)
## ==============================================
## [,1]
## b0 -4.8795477
## x1 -0.5707010
## x2 1.2986183
## x3 0.2016285
## x4 1.2389709
## (x1-k1)+ -2.6071294
## (x2-k1)+ 13.8091009
## (x3-k1)+ -0.2296061
## (x4-k1)+ -1.9437792
## (x1-k2)+ 3.0320808
## (x2-k2)+ -14.1218819
## (x3-k2)+ 0.0236187
## (x4-k2)+ 0.6743852
## (x1-k3)+ 0.7649111
## (x2-k3)+ -19.5997192
## (x3-k3)+ 0.2598065
## (x4-k3)+ 0.9244550
9. TIGA KNOT UJI SIGNIFIKANSI
library(pracma)
uji=function(alpha,para)
{
alpha=0.1
para=0
data=read.table("C:/Users/acer/Downloads/Regnon/data skripsi GWS tanpa kolom kabupaten kota.txt",header=TRUE)
knotall=read.csv("C:/Users/acer/Downloads/dataAll knot 3.csv",header=TRUE)
# Mengurutkan berdasarkan GCV terkecil
knotall=knotall[order(knotall$GCV),]
# Mengambil 12 titik knot dari GCV terkecil
knot=as.matrix(knotall[1,5:16])
data=as.matrix(data)
print(knot)
print(dim(knot))
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="C:/Users/acer/Downloads/Regnon/output uji residual knot3.csv")
write.csv(mx,file="C:/Users/acer/Downloads/Regnon/output uji mx knot3.csv")
write.csv(yhat,file="C:/Users/acer/Downloads/Regnon/output uji yhat knot3.csv")
}
uji(alpha,0)
## X.1 X.2 X.3 X.4 X.5 X.6 X.7 X.8
## 2200 61.13245 61.85408 76.28673 11.38204 11.53673 14.63061 7.692245 10.62041
## X.9 X.10 X.11 X.12
## 2200 69.18367 25.02551 27.44918 75.92265
## [1] 1 12
## ---------------------------------------
## 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
## Tolak Ho yakni variabel bebas
## 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 4.564589e-05
## =============================================
## 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)= 7.291165e-10