#=========================================================#
#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)
folder="D:/Lenovo/Documents/nuget/Materi Kuliah/Semester 5/Pengantar Regresi Nonparametrik/Tugas 6_Regresi Spline Truncated (DATASET 2)/"
data=read.table(file.choose(),header=TRUE)
data=data[,c("Y","X1","X2","X3","X4","X5","X6")] #Y dijadikan kolom pertama
stat.desc(data)
## Y X1 X2 X3 X4
## nbr.val 5.140000e+02 514.0000000 5.140000e+02 5.140000e+02 514.0000000
## nbr.null 0.000000e+00 0.0000000 1.000000e+00 0.000000e+00 1.0000000
## nbr.na 0.000000e+00 0.0000000 0.000000e+00 0.000000e+00 0.0000000
## min 1.414000e+01 2.2700000 0.000000e+00 1.280000e+00 0.0000000
## max 9.637000e+01 40.0100000 9.962000e+01 1.274000e+01 86.8200000
## range 8.223000e+01 37.7400000 9.962000e+01 1.146000e+01 86.8200000
## sum 3.860677e+04 5910.9500000 1.474862e+04 4.549390e+03 2072.5100000
## median 7.902500e+01 9.6150000 2.500500e+01 8.800000e+00 1.2600000
## mean 7.511045e+01 11.4999027 2.869381e+01 8.850953e+00 4.0321206
## SE.mean 6.364823e-01 0.3155689 8.844743e-01 6.816275e-02 0.3855870
## CI.mean.0.95 1.250432e+00 0.6199663 1.737637e+00 1.339125e-01 0.7575239
## var 2.082264e+02 51.1860318 4.020996e+02 2.388127e+00 76.4201684
## std.dev 1.443005e+01 7.1544414 2.005242e+01 1.545357e+00 8.7418630
## coef.var 1.921178e-01 0.6221306 6.988412e-01 1.745977e-01 2.1680559
## X5 X6
## nbr.val 5.140000e+02 5.140000e+02
## nbr.null 0.000000e+00 1.500000e+01
## nbr.na 0.000000e+00 0.000000e+00
## min 5.572000e+01 0.000000e+00
## max 7.793000e+01 5.020000e+01
## range 2.221000e+01 5.020000e+01
## sum 3.608291e+04 1.148700e+04
## median 7.045000e+01 2.220000e+01
## mean 7.020021e+01 2.234825e+01
## SE.mean 1.493169e-01 4.045886e-01
## CI.mean.0.95 2.933478e-01 7.948544e-01
## var 1.145990e+01 8.413767e+01
## std.dev 3.385248e+00 9.172659e+00
## coef.var 4.822276e-02 4.104419e-01
model=(lm(formula=Y~X1+X2+X3+X4+X5+X6,data=data))
summary(model)
##
## Call:
## lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.041 -3.610 2.047 5.431 21.434
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 26.47106 11.16110 2.372 0.0181 *
## X1 -0.61512 0.07038 -8.740 < 2e-16 ***
## X2 -0.21398 0.02453 -8.721 < 2e-16 ***
## X3 -0.21664 0.29665 -0.730 0.4656
## X4 -0.42017 0.04933 -8.518 < 2e-16 ***
## X5 0.92267 0.14182 6.506 1.85e-10 ***
## X6 0.03101 0.04248 0.730 0.4657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.209 on 507 degrees of freedom
## Multiple R-squared: 0.6801, Adjusted R-squared: 0.6763
## F-statistic: 179.7 on 6 and 507 DF, p-value: < 2.2e-16
vif(model)
## X1 X2 X3 X4 X5 X6
## 1.930060 1.842417 1.599757 1.415399 1.754443 1.155668
#scatterplot#
library(ggplot2)
library(gridExtra)
options(scipen = 999)
buat_gg <- function(xvar) {
ggplot(data, aes(x = .data[[xvar]], y = Y)) +
geom_point(shape = 21, color = "steelblue4", fill = "steelblue",
size = 2, alpha = 0.6) +
geom_smooth(method = "loess", formula = y ~ x,
color = "firebrick", fill = "pink", alpha = 0.4, se = TRUE) +
labs(x = xvar, y = "Y") +
theme_bw(base_size = 10) +
theme(plot.margin = margin(8, 10, 8, 8))
}
sp1 <- buat_gg("X1")
sp2 <- buat_gg("X2")
sp3 <- buat_gg("X3")
sp4 <- buat_gg("X4")
sp5 <- buat_gg("X5")
sp6 <- buat_gg("X6")
grid.arrange(sp1, sp2, sp3, sp4, sp5, sp6, ncol = 3)
options(scipen = 0) #dikembalikan agar p-value kecil tampil dalam notasi ilmiah
#======================#
#PEMILIHAN 1 TITIK KNOT#
#======================#
data
## Y X1 X2 X3 X4 X5 X6
## 1 73.84 12.10 34.05 8.50 2.22 64.81 40.2
## 2 76.87 12.45 44.19 10.13 2.20 68.65 32.9
## 3 77.19 13.39 38.66 8.99 1.68 69.07 29.7
## 4 69.08 14.38 26.30 10.04 2.16 69.18 23.6
## 5 78.27 17.86 23.36 10.27 1.46 68.34 33.4
## 6 84.72 13.38 10.06 10.87 1.38 70.14 30.1
## 7 76.69 18.78 31.46 9.29 0.70 67.29 29.5
## 8 77.12 16.64 37.08 9.37 0.56 69.15 25.2
## 9 77.54 17.92 17.76 9.91 1.68 65.61 30.7
## 10 46.10 19.15 40.31 9.02 1.74 67.78 34.1
## 11 80.83 12.12 27.66 9.67 0.46 71.66 32.9
## 12 78.05 15.43 24.17 9.20 1.70 65.48 27.9
## 13 77.06 18.82 20.48 8.39 6.33 65.92 15.4
## 14 79.69 12.42 26.92 8.98 4.54 67.55 34.0
## 15 76.60 17.25 24.00 9.14 2.19 69.60 31.6
## 16 81.16 12.51 14.56 9.58 1.27 70.04 35.9
## 17 46.98 18.31 61.03 10.30 1.15 69.64 32.2
## 18 78.95 18.40 32.95 9.85 0.46 70.60 29.4
## 19 87.96 7.04 0.11 12.74 0.01 72.02 21.7
## 20 73.80 14.59 3.95 10.94 0.20 70.98 25.6
## 21 83.03 10.73 1.81 11.31 0.08 72.06 20.7
## 22 78.05 10.53 6.24 11.27 0.17 69.78 25.6
## 23 41.12 16.41 55.11 8.78 1.40 64.41 29.6
## 24 73.81 11.50 35.19 8.96 1.42 67.90 23.8
## 25 78.14 8.54 51.69 10.19 2.22 69.57 27.4
## 26 77.29 7.01 55.31 9.57 3.22 65.56 15.6
## 27 72.41 15.10 64.15 7.06 0.87 70.34 20.3
## 28 79.92 9.23 26.92 8.97 1.85 69.64 16.9
## 29 82.28 7.98 41.28 10.28 1.51 72.28 24.7
## 30 86.20 3.44 14.42 10.41 0.51 72.31 33.8
## 31 82.42 7.87 41.91 9.78 1.78 72.20 17.7
## 32 76.65 8.21 24.26 9.18 1.82 69.09 11.0
## 33 76.96 7.99 41.36 9.69 1.39 70.77 20.2
## 34 82.50 7.47 36.67 10.08 2.02 70.11 32.6
## 35 84.33 8.04 25.69 10.50 2.06 71.24 28.0
## 36 70.54 8.86 55.14 9.37 3.08 63.47 20.7
## 37 65.87 16.39 63.23 6.99 1.20 69.58 31.8
## 38 75.32 7.54 56.82 9.90 3.57 66.89 28.9
## 39 81.16 8.69 49.20 9.93 3.18 70.27 18.4
## 40 80.35 11.66 47.18 9.20 2.85 72.24 22.4
## 41 85.35 7.44 33.30 9.14 1.51 69.59 14.4
## 42 78.37 11.38 35.54 9.13 0.72 67.88 17.7
## 43 75.14 8.79 42.52 9.74 3.26 67.82 21.8
## 44 76.75 7.89 39.91 9.86 2.67 67.71 17.7
## 45 57.85 8.06 24.70 9.32 2.79 69.54 16.0
## 46 80.30 9.08 37.52 9.28 2.18 70.25 9.6
## 47 73.17 21.79 46.34 7.01 1.41 70.24 20.3
## 48 65.35 22.81 85.31 7.19 0.69 69.96 28.9
## 49 90.25 8.00 2.92 11.59 0.02 73.93 5.8
## 50 89.82 7.24 8.31 11.56 0.04 74.75 7.7
## 51 81.03 11.42 4.27 10.45 0.01 70.18 10.6
## 52 74.59 12.21 1.10 10.13 0.07 64.28 5.7
## 53 88.22 4.79 5.42 11.07 0.05 73.13 19.4
## 54 77.31 9.49 28.45 10.88 0.04 71.63 10.4
## 55 71.32 6.85 37.86 11.11 0.14 70.20 26.1
## 56 66.12 14.78 20.47 8.52 0.25 72.09 18.9
## 57 86.71 7.34 15.76 9.13 2.78 71.52 27.0
## 58 82.56 7.13 32.61 8.57 3.29 69.56 25.4
## 59 82.05 5.88 27.30 9.23 3.62 67.02 28.5
## 60 89.40 4.16 15.74 9.56 1.37 70.84 18.5
## 61 85.08 6.34 23.54 8.86 0.92 69.70 19.4
## 62 86.92 6.60 20.57 9.63 1.95 73.23 20.1
## 63 82.59 6.80 33.84 8.37 3.12 70.30 28.6
## 64 82.09 6.80 33.26 8.37 2.69 68.29 29.4
## 65 53.86 13.72 72.89 7.85 8.11 65.10 33.7
## 66 84.98 5.56 26.93 9.23 2.33 72.24 17.7
## 67 85.60 6.45 26.75 9.36 6.80 68.71 14.7
## 68 80.91 6.92 34.36 9.69 2.10 68.53 29.7
## 69 89.63 4.17 4.19 11.60 0.10 74.16 24.2
## 70 92.90 3.05 2.11 11.27 0.05 74.39 16.3
## 71 83.45 2.27 6.49 10.45 0.44 70.69 19.5
## 72 87.05 5.24 12.33 12.23 0.05 73.23 15.8
## 73 92.86 4.11 2.64 11.77 0.02 75.13 20.1
## 74 89.33 5.44 3.37 11.25 0.09 74.43 19.8
## 75 84.08 4.20 14.32 10.93 0.07 70.95 17.7
## 76 59.69 7.04 24.67 9.44 4.21 71.42 7.6
## 77 60.71 6.06 16.34 8.61 5.19 70.69 12.7
## 78 65.75 6.31 36.88 9.51 4.00 71.79 17.9
## 79 69.30 5.64 73.12 7.77 5.86 68.62 18.8
## 80 71.59 8.15 13.63 9.07 7.97 72.00 10.1
## 81 58.88 9.72 15.70 9.26 4.42 70.84 15.9
## 82 68.67 7.07 37.28 8.89 3.59 70.99 16.6
## 83 77.06 5.23 9.90 10.01 4.77 71.64 10.4
## 84 73.80 8.07 28.94 9.01 4.29 69.13 23.0
## 85 71.15 22.98 89.58 8.47 5.21 68.41 19.6
## 86 92.20 3.16 7.83 11.82 0.08 73.02 8.7
## 87 76.30 3.21 13.10 10.35 1.18 71.67 14.9
## 88 83.76 7.54 32.30 8.25 3.39 70.55 8.7
## 89 75.01 8.90 50.04 8.61 3.48 71.75 14.9
## 90 65.64 8.54 49.63 8.20 4.47 69.70 4.8
## 91 74.13 9.45 28.95 8.42 4.41 71.22 10.1
## 92 68.51 4.43 29.18 9.09 3.41 71.84 12.0
## 93 69.60 9.79 53.10 8.59 4.23 68.67 14.1
## 94 69.99 10.85 49.09 8.15 4.79 66.98 23.7
## 95 73.39 5.29 26.74 8.53 2.66 68.43 13.7
## 96 71.36 6.46 28.31 8.40 5.11 70.49 22.7
## 97 85.88 8.24 9.34 11.34 0.04 73.29 13.5
## 98 87.06 3.00 8.00 10.04 0.32 72.86 4.1
## 99 73.01 11.46 27.60 9.40 1.75 68.85 15.7
## 100 79.54 13.15 31.03 7.51 6.17 69.31 32.5
## 101 75.62 10.93 30.27 8.12 2.57 69.72 25.9
## 102 75.19 15.00 40.77 8.73 1.94 66.87 7.8
## 103 77.66 14.13 39.52 7.85 3.73 68.95 21.9
## 104 79.48 14.90 28.89 8.14 6.17 69.51 16.5
## 105 79.17 9.58 55.26 7.60 5.17 69.76 20.4
## 106 85.17 9.99 37.40 8.40 1.46 69.78 9.3
## 107 75.37 10.36 49.77 8.40 2.53 67.88 23.0
## 108 77.55 13.28 28.50 8.35 1.34 66.33 22.9
## 109 67.42 11.80 71.69 7.83 2.38 65.73 32.6
## 110 78.79 10.91 31.52 7.53 1.98 68.96 15.4
## 111 59.58 18.26 44.73 7.70 6.83 66.37 33.1
## 112 83.62 10.22 2.15 11.09 0.04 71.99 18.9
## 113 57.79 8.88 41.12 9.68 0.87 67.70 23.3
## 114 66.10 12.65 33.15 9.51 0.26 70.25 17.5
## 115 68.11 11.23 40.09 10.31 0.24 71.27 15.4
## 116 74.91 17.51 59.02 9.53 1.05 68.40 24.0
## 117 73.11 14.79 55.84 9.31 1.16 69.22 28.6
## 118 72.62 11.29 43.03 8.22 3.15 68.77 21.6
## 119 74.75 17.83 46.75 8.25 3.87 67.27 14.3
## 120 71.61 18.00 57.84 8.20 2.06 68.35 26.0
## 121 76.57 10.76 36.02 8.87 4.18 67.30 27.1
## 122 77.06 11.15 44.85 8.49 3.00 63.93 15.7
## 123 72.43 14.12 55.68 8.41 0.95 68.46 22.1
## 124 62.14 9.40 60.54 7.57 1.55 68.52 23.2
## 125 78.72 14.71 12.26 11.90 0.05 70.74 6.7
## 126 84.46 12.79 20.90 8.07 0.92 69.97 10.3
## 127 84.93 10.65 31.53 8.00 1.41 70.29 16.7
## 128 76.82 17.17 60.90 8.55 1.23 69.83 23.5
## 129 76.18 11.17 52.76 8.57 1.90 68.13 24.6
## 130 88.78 8.04 22.06 7.89 2.34 70.42 9.8
## 131 76.40 10.52 53.84 7.35 2.13 69.23 17.1
## 132 85.22 13.80 30.49 8.34 2.07 71.25 14.2
## 133 79.61 11.02 44.87 7.83 2.51 69.91 22.7
## 134 80.85 12.89 43.85 8.21 0.94 69.74 10.0
## 135 87.35 9.14 11.49 8.56 0.37 71.05 15.8
## 136 88.18 6.73 22.71 7.20 3.02 68.75 5.0
## 137 83.74 7.25 55.59 8.05 1.35 70.42 10.5
## 138 75.77 13.49 59.30 8.75 4.24 64.29 16.1
## 139 84.64 7.77 8.71 10.94 0.03 71.91 13.4
## 140 85.78 7.28 9.50 10.80 0.04 72.12 7.1
## 141 74.35 4.32 11.86 8.79 1.69 71.60 23.2
## 142 62.34 6.46 8.85 9.16 2.66 71.54 20.8
## 143 84.87 3.11 23.99 7.49 5.04 68.98 20.7
## 144 60.38 5.29 14.30 7.48 2.89 72.11 18.2
## 145 60.25 2.71 20.73 7.87 3.40 70.43 20.6
## 146 62.02 6.73 14.77 9.09 3.85 72.59 17.3
## 147 87.27 4.27 7.03 10.39 0.07 73.96 20.7
## 148 60.42 5.90 20.42 9.82 1.26 71.09 21.6
## 149 62.68 5.95 12.19 9.20 0.84 72.01 17.9
## 150 60.51 5.25 8.26 9.28 2.67 66.22 16.1
## 151 53.02 11.26 23.63 7.92 3.31 63.50 20.5
## 152 55.34 6.95 45.37 8.25 1.16 68.10 15.2
## 153 88.83 5.02 2.72 11.12 0.20 73.94 16.1
## 154 83.22 7.95 12.38 10.90 0.10 72.85 15.2
## 155 58.05 13.13 1.54 9.01 0.05 69.61 18.6
## 156 91.63 4.68 2.40 11.33 0.00 74.81 19.1
## 157 87.42 6.78 0.86 10.65 0.02 73.59 19.8
## 158 91.35 4.09 2.09 10.94 0.01 74.16 17.1
## 159 90.52 3.10 13.99 11.63 0.01 74.78 16.6
## 160 91.81 4.20 8.10 11.77 0.01 75.12 16.8
## 161 70.01 7.27 31.89 8.67 0.29 71.92 27.6
## 162 78.59 7.01 44.67 7.62 0.87 71.83 27.0
## 163 80.60 10.22 32.56 7.13 0.84 70.79 11.4
## 164 80.63 6.40 22.97 9.28 0.31 74.27 29.2
## 165 78.99 9.77 46.60 8.04 0.63 72.07 24.1
## 166 78.47 10.28 35.61 8.24 0.98 70.19 20.7
## 167 83.20 7.42 21.45 8.32 0.64 72.57 25.4
## 168 81.43 12.12 29.18 7.94 0.39 74.29 23.4
## 169 82.77 11.20 15.84 7.97 0.18 72.76 22.9
## 170 82.58 11.21 26.26 7.74 0.43 71.05 24.1
## 171 86.88 9.36 17.60 8.89 0.61 73.19 14.4
## 172 86.52 12.13 6.41 7.09 0.57 72.46 18.4
## 173 85.65 9.52 27.91 7.67 0.65 73.24 18.7
## 174 80.19 8.46 29.17 8.04 0.30 71.74 24.0
## 175 89.04 7.87 11.45 8.13 0.27 72.90 17.1
## 176 87.48 4.93 6.36 9.79 0.13 74.30 23.2
## 177 75.79 10.52 24.46 8.16 0.48 73.10 25.1
## 178 83.86 8.98 22.31 8.16 0.98 72.17 23.9
## 179 85.53 6.67 9.88 10.36 0.02 74.45 18.2
## 180 80.90 7.50 12.32 10.26 0.02 73.11 26.9
## 181 91.55 3.96 4.68 10.91 0.01 75.04 16.3
## 182 83.78 9.16 3.34 10.35 0.01 73.08 19.9
## 183 93.90 4.10 7.14 11.69 0.02 75.79 10.3
## 184 89.81 2.38 19.52 11.42 0.03 75.23 14.3
## 185 88.68 4.66 6.98 11.12 0.01 74.80 24.5
## 186 74.37 11.53 16.47 9.64 0.06 72.92 27.1
## 187 78.92 6.14 9.82 8.77 0.12 71.80 23.6
## 188 84.69 10.99 27.64 7.46 0.54 74.25 18.5
## 189 80.74 12.53 30.51 8.04 0.22 73.98 20.9
## 190 78.47 14.99 28.37 7.65 0.30 73.37 26.0
## 191 79.40 14.90 23.50 7.15 0.51 74.47 19.9
## 192 81.21 16.34 33.23 7.74 0.33 73.83 21.9
## 193 86.22 11.33 23.10 8.34 0.49 75.21 20.6
## 194 77.87 15.58 20.14 7.07 0.59 72.17 29.2
## 195 80.20 10.96 26.04 7.80 0.48 74.20 25.8
## 196 87.70 9.81 19.18 7.98 0.37 76.23 21.5
## 197 85.64 12.28 28.11 9.08 0.17 77.07 24.5
## 198 90.93 7.58 17.61 9.91 0.12 77.86 24.3
## 199 88.91 10.94 21.08 7.66 0.72 76.56 19.5
## 200 89.67 9.79 12.04 9.05 0.32 77.72 22.2
## 201 89.37 12.87 8.69 7.83 0.29 75.97 18.4
## 202 87.29 11.72 9.41 7.35 0.53 75.04 20.2
## 203 88.53 11.49 8.14 7.36 0.84 74.71 21.2
## 204 86.85 14.17 7.08 7.86 0.51 74.77 19.5
## 205 90.88 9.31 9.89 7.96 0.38 76.39 18.5
## 206 88.81 7.24 16.34 9.32 0.13 76.86 15.7
## 207 87.00 6.61 23.56 8.39 0.36 76.04 18.9
## 208 90.47 12.01 4.64 8.47 0.38 75.60 9.5
## 209 85.57 7.17 21.70 8.46 0.38 75.95 18.8
## 210 78.91 9.26 12.34 7.87 0.43 75.77 25.1
## 211 84.83 9.39 19.24 7.96 0.35 74.58 22.4
## 212 84.77 8.92 20.88 7.38 0.46 74.85 24.7
## 213 80.65 9.67 26.99 7.70 0.41 73.87 28.6
## 214 81.81 15.03 21.40 6.84 0.40 73.85 15.3
## 215 83.71 7.30 18.31 7.55 0.27 72.00 21.5
## 216 78.79 15.78 25.43 6.68 0.43 69.96 21.6
## 217 91.41 6.11 6.68 11.08 0.01 77.22 15.4
## 218 87.67 8.44 10.65 10.91 0.01 77.63 16.0
## 219 91.35 4.66 10.74 11.30 0.03 77.93 16.9
## 220 92.74 4.23 4.00 10.57 0.03 77.90 15.7
## 221 78.27 6.81 21.88 9.53 0.03 74.60 28.2
## 222 85.73 7.68 1.64 9.90 0.02 74.77 22.3
## 223 85.73 15.64 18.44 9.17 0.33 75.29 21.2
## 224 82.92 11.96 18.86 9.97 0.15 73.94 20.5
## 225 82.97 15.60 41.29 7.32 0.78 74.24 22.2
## 226 84.91 7.52 17.85 11.10 0.07 75.08 12.4
## 227 84.21 6.49 26.33 12.14 0.01 74.91 16.8
## 228 82.94 13.65 28.82 7.95 1.04 72.86 20.9
## 229 88.74 9.53 17.71 7.63 0.49 73.55 16.1
## 230 83.17 10.63 36.23 7.80 0.68 74.64 15.8
## 231 87.35 6.53 25.05 8.72 0.29 74.91 21.5
## 232 86.43 8.69 18.29 7.95 0.71 74.34 20.3
## 233 81.19 10.72 31.55 8.37 0.46 73.27 16.8
## 234 81.09 9.45 28.02 7.97 0.53 73.26 19.5
## 235 79.58 8.93 33.64 7.36 0.62 70.96 29.9
## 236 78.97 9.51 32.29 6.66 0.58 70.03 29.7
## 237 85.77 7.34 17.57 7.82 0.92 71.38 21.9
## 238 78.44 13.34 35.14 6.40 0.69 67.60 17.0
## 239 80.22 11.90 36.32 6.90 0.79 69.94 4.1
## 240 74.22 17.19 34.59 6.49 0.71 68.12 35.4
## 241 80.31 9.24 20.80 7.56 0.42 70.81 27.9
## 242 84.49 5.00 12.26 10.82 0.09 74.69 8.4
## 243 87.32 9.80 11.52 9.09 0.30 73.25 16.2
## 244 86.43 9.15 16.41 8.65 0.29 73.22 18.0
## 245 84.88 10.89 26.30 8.27 0.52 72.28 17.1
## 246 86.22 11.04 20.54 7.81 0.65 72.28 16.5
## 247 87.98 9.80 17.53 8.48 0.38 73.29 19.2
## 248 87.53 14.40 11.10 7.81 0.56 73.20 14.0
## 249 87.63 12.18 10.93 7.44 0.67 72.57 14.1
## 250 85.85 14.91 7.08 7.46 0.79 72.36 17.8
## 251 90.30 12.42 4.97 8.29 0.40 73.22 9.4
## 252 90.90 10.96 0.79 10.20 0.27 73.30 15.4
## 253 76.02 19.35 34.55 6.04 0.49 70.79 10.2
## 254 78.97 21.76 11.31 5.46 0.55 68.64 14.1
## 255 79.08 13.85 10.80 7.28 0.32 68.31 25.1
## 256 79.77 18.70 21.19 5.70 0.66 72.47 16.7
## 257 76.67 7.15 40.23 10.77 0.02 74.67 18.6
## 258 80.80 7.30 27.37 10.87 0.03 74.66 17.7
## 259 88.76 4.26 11.17 11.19 0.02 74.13 17.3
## 260 78.62 6.48 8.81 9.60 0.05 70.99 31.8
## 261 84.16 6.60 8.66 9.73 0.04 72.31 11.7
## 262 89.30 5.77 9.94 11.32 0.01 74.10 11.0
## 263 92.29 4.74 1.25 11.53 0.02 73.44 12.8
## 264 93.06 4.65 0.00 10.71 0.02 74.75 1.6
## 265 81.31 3.31 20.60 10.05 0.19 73.29 23.1
## 266 75.63 9.27 53.73 7.26 1.14 65.58 28.6
## 267 72.76 8.68 58.52 6.95 1.01 68.13 35.5
## 268 77.60 6.93 13.62 9.16 0.11 70.65 26.4
## 269 81.39 4.85 27.36 8.14 0.50 65.60 23.9
## 270 85.92 5.89 7.38 10.76 0.02 72.24 17.6
## 271 80.12 3.98 11.17 10.32 0.08 67.39 22.0
## 272 71.95 6.20 25.44 8.72 0.10 68.98 22.3
## 273 88.67 2.57 20.29 11.44 0.02 73.11 9.2
## 274 87.64 4.96 18.29 8.40 0.77 73.20 8.7
## 275 92.51 4.70 15.09 9.00 0.35 74.48 6.3
## 276 92.90 2.30 7.11 10.88 0.09 75.88 4.9
## 277 92.77 4.47 5.08 9.66 0.13 74.52 6.3
## 278 87.47 5.61 13.91 8.09 0.21 72.28 4.9
## 279 77.19 5.28 37.96 7.67 0.34 71.33 10.2
## 280 83.97 6.56 24.50 6.22 0.60 71.25 6.4
## 281 83.26 5.85 16.00 7.35 0.55 72.70 6.2
## 282 96.37 2.68 0.49 11.37 0.01 75.69 10.8
## 283 75.74 13.67 28.81 7.24 0.52 68.08 19.9
## 284 79.50 12.93 22.48 6.97 0.37 67.12 17.6
## 285 73.32 15.63 35.08 7.59 0.35 66.95 27.6
## 286 85.82 13.91 7.04 8.83 2.41 68.54 25.7
## 287 81.24 12.62 37.61 9.07 1.26 67.76 12.4
## 288 78.26 14.39 25.58 8.83 1.41 67.27 36.7
## 289 86.47 12.95 9.76 8.88 1.83 69.20 10.5
## 290 68.39 25.80 37.84 6.36 0.87 68.17 33.5
## 291 80.15 8.62 6.16 9.45 0.01 72.55 25.6
## 292 75.53 8.67 16.57 10.82 0.19 71.19 31.8
## 293 70.45 21.78 35.11 7.85 3.69 65.64 38.4
## 294 67.02 25.18 32.28 7.49 2.37 66.89 50.1
## 295 73.61 21.85 24.02 8.48 1.63 67.61 42.7
## 296 71.21 14.30 36.62 7.65 0.96 65.63 48.1
## 297 63.81 19.97 55.04 8.26 1.87 62.35 39.3
## 298 68.61 11.77 27.31 7.91 1.05 65.96 37.2
## 299 77.22 12.56 20.42 7.32 0.90 68.30 33.3
## 300 72.93 22.86 21.14 8.51 1.11 66.12 27.5
## 301 82.72 12.06 24.86 9.01 1.45 68.71 21.3
## 302 75.69 19.69 17.70 7.83 0.59 67.63 36.8
## 303 70.98 28.08 43.49 8.11 4.48 65.82 26.3
## 304 71.31 27.17 26.10 7.33 0.89 67.57 42.5
## 305 70.11 24.78 37.63 7.68 1.47 67.87 35.1
## 306 71.27 27.05 26.04 8.03 1.62 65.60 39.8
## 307 78.00 16.82 32.38 8.20 1.91 68.00 36.2
## 308 75.20 12.33 43.78 8.29 1.56 67.91 24.9
## 309 65.96 31.78 57.42 7.40 2.89 68.87 39.5
## 310 62.53 27.48 63.65 6.95 1.14 68.99 44.3
## 311 66.92 25.06 48.86 8.22 1.22 68.49 43.7
## 312 56.66 28.37 53.35 7.47 0.92 61.06 36.9
## 313 71.78 14.42 43.40 7.42 0.90 65.67 47.7
## 314 75.94 8.61 21.61 11.36 0.04 70.52 29.9
## 315 70.83 7.08 92.51 6.71 3.67 69.76 30.8
## 316 73.85 5.21 81.54 7.70 2.55 71.74 27.2
## 317 78.51 4.79 43.57 7.70 9.04 71.77 22.1
## 318 77.17 9.25 45.82 8.05 12.71 71.45 0.0
## 319 76.17 8.18 38.35 7.89 12.07 72.41 24.8
## 320 74.43 8.16 41.91 8.16 22.35 72.88 16.7
## 321 81.24 6.28 40.03 8.15 4.09 74.20 32.7
## 322 76.96 9.97 54.82 7.60 7.32 73.77 31.0
## 323 73.84 5.90 67.93 7.18 7.16 72.74 12.2
## 324 59.56 11.12 46.91 7.63 10.67 73.32 35.3
## 325 76.53 9.13 60.21 6.58 7.70 69.22 22.9
## 326 68.90 4.23 85.45 7.50 5.49 71.26 25.4
## 327 72.55 4.45 57.92 10.34 0.03 73.87 16.7
## 328 68.27 4.70 47.34 8.50 0.30 72.81 20.1
## 329 62.27 4.18 13.52 8.94 8.01 71.28 17.9
## 330 73.03 5.69 26.11 7.92 8.45 70.39 35.5
## 331 82.66 5.21 36.50 8.50 11.12 69.23 16.2
## 332 67.09 4.72 34.12 9.40 6.39 67.75 23.9
## 333 81.65 5.35 19.30 8.95 9.96 71.68 15.3
## 334 77.99 4.99 34.52 9.25 21.48 66.43 34.0
## 335 65.03 7.12 30.94 7.78 18.29 69.63 25.8
## 336 61.78 3.96 23.05 8.51 8.22 72.03 29.1
## 337 59.34 3.12 24.79 8.63 12.59 69.84 13.2
## 338 56.38 5.47 43.53 9.48 10.31 70.96 12.9
## 339 81.88 4.58 44.92 8.58 14.03 68.60 24.0
## 340 51.29 6.44 49.56 8.38 23.69 69.94 21.0
## 341 81.16 6.63 33.78 9.96 4.63 68.91 21.7
## 342 80.69 3.44 8.72 11.69 1.13 73.70 28.0
## 343 80.59 3.73 28.97 8.09 2.79 70.11 41.7
## 344 75.34 4.32 22.59 7.66 8.06 69.77 20.1
## 345 81.54 2.44 28.92 8.18 2.01 68.01 30.1
## 346 82.85 4.60 35.70 7.90 2.63 66.80 15.9
## 347 88.58 3.19 24.05 8.13 2.40 71.16 14.4
## 348 83.19 4.01 30.53 8.19 1.31 66.90 25.4
## 349 81.30 5.84 40.01 8.46 1.25 66.86 13.0
## 350 81.20 6.25 27.66 8.06 0.95 64.97 36.0
## 351 85.94 5.77 20.34 9.21 3.12 71.28 18.1
## 352 83.15 4.12 11.33 8.53 2.95 70.94 25.1
## 353 82.70 5.22 12.79 8.36 2.35 68.40 33.4
## 354 85.62 4.63 0.01 10.05 0.02 71.89 26.5
## 355 86.34 3.92 13.93 11.05 0.18 72.62 12.4
## 356 80.04 9.11 16.49 9.35 7.27 72.99 22.4
## 357 87.17 7.61 2.51 9.15 10.69 72.75 17.6
## 358 81.91 5.54 6.56 9.54 14.55 72.41 23.0
## 359 60.87 9.72 21.44 8.90 10.26 73.19 22.0
## 360 61.64 9.06 12.47 9.40 13.67 73.57 29.0
## 361 88.11 6.97 4.31 8.65 3.81 71.83 24.6
## 362 56.53 11.38 12.76 9.10 47.25 72.46 0.0
## 363 91.23 2.31 2.06 10.67 0.13 74.89 21.6
## 364 89.68 4.81 0.89 11.41 0.16 74.68 24.4
## 365 88.89 4.11 0.01 11.12 0.13 74.67 27.4
## 366 76.61 8.99 18.80 9.06 13.78 72.78 22.6
## 367 74.67 6.54 27.49 9.22 42.89 71.51 20.5
## 368 78.80 5.53 30.10 8.53 13.26 71.42 15.8
## 369 61.25 4.62 16.79 8.66 10.43 71.53 15.1
## 370 87.53 6.10 3.61 10.57 0.17 74.08 14.8
## 371 84.02 7.37 28.53 8.81 3.02 70.07 24.8
## 372 79.61 6.87 20.89 10.61 0.71 71.82 23.1
## 373 60.47 11.01 27.43 9.23 0.65 70.78 19.0
## 374 79.36 8.46 29.79 10.01 1.41 70.91 19.3
## 375 72.82 8.89 21.34 9.76 2.02 70.71 26.4
## 376 78.79 6.65 10.85 10.49 0.78 71.99 10.9
## 377 75.82 11.84 13.73 9.18 1.40 70.86 15.0
## 378 80.12 7.90 26.84 8.92 2.91 68.33 27.8
## 379 53.11 8.76 49.95 9.68 0.44 71.59 24.9
## 380 75.22 5.80 23.25 8.61 1.14 68.51 28.4
## 381 71.81 12.04 27.48 8.65 3.16 64.95 33.0
## 382 87.56 5.79 2.91 11.25 0.03 72.50 21.8
## 383 82.00 6.56 4.42 10.11 0.29 71.68 19.5
## 384 83.92 5.60 19.78 11.24 0.19 72.84 10.5
## 385 80.21 5.03 17.60 10.08 0.11 71.34 20.5
## 386 85.14 6.94 30.04 8.83 3.99 71.02 29.1
## 387 84.60 15.16 15.74 9.53 5.38 71.33 26.5
## 388 73.97 16.25 39.70 8.20 4.74 67.86 34.1
## 389 79.39 12.85 31.23 8.87 2.60 66.78 29.0
## 390 79.97 13.36 17.55 9.58 4.27 69.73 30.0
## 391 81.44 12.31 14.99 8.86 2.97 69.37 26.0
## 392 66.50 12.90 14.24 8.79 2.79 66.87 27.7
## 393 78.91 14.91 30.86 8.34 3.29 64.62 28.5
## 394 73.78 16.74 13.32 8.93 4.41 66.41 21.3
## 395 80.06 12.83 43.58 9.39 4.20 70.35 26.4
## 396 51.91 14.15 25.14 8.61 0.96 66.08 25.6
## 397 81.21 12.85 26.66 8.88 8.65 69.97 24.7
## 398 82.43 6.56 13.53 11.43 0.10 71.45 22.1
## 399 72.03 12.27 14.34 8.62 0.77 69.04 31.3
## 400 83.73 7.22 22.48 8.38 0.51 68.84 33.7
## 401 86.59 9.18 12.43 7.53 0.28 71.11 15.8
## 402 81.67 13.06 11.07 7.75 0.39 66.99 36.3
## 403 84.63 8.29 10.38 8.00 0.27 67.90 35.4
## 404 87.53 7.42 20.29 8.84 0.74 70.89 21.1
## 405 82.31 8.55 22.62 8.19 0.48 67.92 33.5
## 406 85.66 10.53 22.10 7.85 1.26 67.85 26.0
## 407 83.73 9.65 21.50 8.15 0.74 69.45 34.7
## 408 81.00 13.40 24.18 8.87 0.40 67.36 30.0
## 409 86.67 8.46 19.02 8.91 0.99 69.55 22.1
## 410 87.90 7.48 22.95 8.58 0.86 70.50 24.0
## 411 86.65 6.73 20.04 7.98 1.47 68.10 27.4
## 412 87.92 5.14 24.51 8.43 0.99 70.74 26.4
## 413 86.94 8.90 22.62 8.82 0.91 70.49 17.6
## 414 80.31 12.69 26.74 8.96 1.23 71.34 34.9
## 415 82.23 12.71 34.16 9.19 1.23 71.00 32.1
## 416 76.88 12.48 48.81 8.63 1.49 73.99 36.9
## 417 85.66 12.66 20.67 8.52 4.12 69.36 15.5
## 418 87.73 6.93 20.83 8.83 3.09 71.19 26.0
## 419 80.25 12.12 40.06 9.03 0.84 73.83 28.7
## 420 87.95 5.07 2.39 11.33 0.01 72.60 25.6
## 421 82.76 5.34 6.50 10.82 0.07 71.78 26.7
## 422 83.05 7.69 3.66 11.29 0.14 71.39 25.5
## 423 83.69 11.80 17.27 9.82 1.94 71.42 23.8
## 424 85.89 13.02 17.39 9.64 2.07 70.44 27.8
## 425 71.93 14.07 27.59 8.24 0.83 70.56 29.2
## 426 72.59 13.77 17.89 8.81 1.64 68.56 37.2
## 427 84.41 11.26 10.60 8.71 3.16 71.02 33.6
## 428 85.33 10.73 5.66 8.39 3.27 69.34 30.4
## 429 62.54 14.81 30.22 8.62 0.50 70.72 31.9
## 430 70.26 13.57 9.23 9.22 2.58 70.47 31.8
## 431 79.71 13.48 12.53 9.54 4.72 69.68 24.7
## 432 81.40 14.06 11.93 8.94 2.36 71.02 33.9
## 433 85.54 14.04 19.30 8.83 4.62 73.01 31.3
## 434 54.01 15.90 21.11 9.33 1.80 68.52 31.3
## 435 79.06 14.03 34.15 8.58 1.15 70.46 24.4
## 436 67.07 15.43 10.97 6.88 1.07 67.90 36.8
## 437 70.01 14.76 48.86 7.89 0.56 67.86 37.1
## 438 89.67 4.59 3.79 12.08 0.08 74.07 25.7
## 439 80.49 7.53 5.69 11.10 0.24 71.51 29.7
## 440 79.50 17.48 16.35 8.63 1.06 68.09 34.7
## 441 81.91 18.38 16.23 8.21 1.97 69.92 16.0
## 442 80.32 15.51 11.21 8.91 1.84 69.24 27.1
## 443 81.71 17.64 3.40 8.28 4.97 64.95 18.4
## 444 78.88 17.03 13.44 7.85 2.24 66.67 30.5
## 445 86.50 5.64 3.53 10.58 0.04 73.25 23.6
## 446 57.34 4.79 17.14 8.51 4.10 67.33 27.9
## 447 81.47 7.57 24.96 8.31 2.25 68.66 32.8
## 448 73.19 14.31 70.51 8.58 2.49 71.45 37.6
## 449 75.87 16.08 30.23 8.01 0.90 63.20 28.1
## 450 65.80 14.54 9.18 9.66 0.77 62.56 30.5
## 451 78.29 7.32 20.96 7.88 2.99 69.37 27.9
## 452 71.92 17.84 22.59 10.13 4.12 67.04 29.4
## 453 49.43 21.79 32.96 9.74 1.27 65.78 34.0
## 454 49.94 24.47 24.47 10.35 4.93 63.99 25.1
## 455 80.21 16.53 28.66 9.34 4.57 67.03 20.3
## 456 70.57 21.08 39.59 8.88 4.38 60.46 27.5
## 457 61.29 22.39 45.87 9.78 4.36 62.61 31.4
## 458 43.15 24.21 30.55 9.43 11.36 63.57 40.6
## 459 65.95 28.78 39.73 9.31 5.56 63.39 29.9
## 460 53.32 15.28 36.60 8.92 4.08 66.96 35.5
## 461 84.56 5.25 13.20 12.05 0.09 71.37 20.7
## 462 59.11 20.68 12.48 10.70 0.35 66.65 32.0
## 463 53.07 8.74 42.60 8.65 3.39 66.91 26.1
## 464 52.96 11.44 35.58 9.48 3.12 65.06 29.5
## 465 62.34 4.62 27.66 8.87 2.80 70.11 24.3
## 466 52.69 5.68 38.45 8.41 4.27 66.50 30.4
## 467 51.76 8.17 47.20 9.34 4.27 63.94 18.8
## 468 79.63 12.47 36.14 8.88 8.88 69.90 19.0
## 469 60.71 5.38 33.98 7.91 4.00 68.10 11.7
## 470 48.28 7.31 45.09 8.62 5.23 62.79 30.6
## 471 88.24 3.39 6.54 11.74 0.12 71.70 21.1
## 472 64.68 6.35 37.75 10.17 1.44 70.10 21.3
## 473 84.07 10.01 23.56 10.06 32.29 67.60 23.8
## 474 30.40 34.71 77.16 6.19 3.68 60.50 33.2
## 475 66.89 11.45 38.66 9.99 13.14 67.78 26.9
## 476 66.97 23.35 21.23 10.45 10.20 68.73 22.2
## 477 46.80 25.89 27.44 9.16 3.34 69.63 39.9
## 478 49.23 23.53 39.30 10.17 1.82 68.74 33.7
## 479 24.19 35.60 88.21 4.51 29.20 65.90 0.0
## 480 22.27 35.39 99.01 3.97 7.96 67.22 0.0
## 481 83.64 13.55 20.61 9.96 9.57 72.83 24.7
## 482 63.07 13.21 40.22 10.03 37.92 67.07 26.5
## 483 54.21 15.68 75.39 8.41 18.50 67.24 34.4
## 484 25.07 29.79 79.27 3.52 39.40 64.99 0.0
## 485 25.62 36.08 66.40 4.13 45.09 66.42 0.0
## 486 22.83 31.57 97.09 2.38 10.42 66.33 0.0
## 487 47.38 29.16 59.99 8.86 35.11 66.92 26.2
## 488 43.19 19.80 39.30 9.16 35.64 60.97 31.4
## 489 35.44 25.72 81.80 6.96 45.27 65.90 28.6
## 490 45.99 24.36 83.74 6.30 40.28 59.19 0.0
## 491 32.69 36.99 63.52 9.27 1.64 66.66 37.7
## 492 23.09 29.63 98.98 8.59 86.82 58.54 40.0
## 493 21.51 35.42 99.00 5.00 32.55 64.14 33.0
## 494 37.80 30.97 50.98 6.78 17.30 65.95 0.0
## 495 21.33 36.94 99.62 4.28 8.67 66.51 39.5
## 496 17.10 37.09 98.15 4.64 26.28 55.72 0.0
## 497 30.12 36.44 21.79 1.43 53.48 66.40 0.0
## 498 31.37 29.20 66.66 4.30 12.95 66.39 0.0
## 499 14.14 40.01 92.00 1.28 60.62 66.12 0.0
## 500 21.35 38.66 92.00 2.23 20.93 65.93 50.2
## 501 78.35 10.50 9.55 11.68 0.35 71.07 21.3
## 502 50.88 26.88 52.18 9.07 6.37 66.89 27.3
## 503 71.81 18.73 43.21 9.64 2.07 69.55 0.0
## 504 48.20 21.38 50.37 10.78 10.15 69.02 30.5
## 505 40.39 18.11 79.81 8.23 10.30 67.01 31.3
## 506 54.42 16.76 49.35 8.55 9.53 65.41 30.9
## 507 43.12 28.24 42.97 9.90 24.62 61.76 19.6
## 508 36.76 28.90 72.04 8.44 10.80 60.86 19.7
## 509 48.19 14.57 43.79 9.86 24.48 65.62 25.7
## 510 45.92 31.23 47.00 8.33 56.66 61.14 31.8
## 511 41.27 30.28 61.00 9.97 15.17 65.82 0.0
## 512 74.91 27.80 18.33 8.48 5.89 68.09 20.4
## 513 35.53 32.29 58.55 6.35 20.75 67.71 34.7
## 514 76.32 14.41 5.16 11.63 0.17 71.92 31.0
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)
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),])
colnames(knot1)=paste0("k1_x",1:m)
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)
{
data1[,j]=pmax(data[,(j+para+1)]-knot1[i,j],0)
}
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
A2=((N-sum((mx%*%C)*mx))/N)^2 #sama dengan (sum(diag(I-A))/N)^2
GCV[i]=MSE[i]/A2
}
dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot1))
dataG=dataAll[order(GCV),-2]
write.csv(dataAll,file=paste0(folder,"dataAll knot 1.csv"))
cat("==============================================","\n")
cat("HASIL GCV terkecil dengan 1 knot","\n")
cat("==============================================","\n")
print(dataG[1,1:(2+m)])
cat("Nilai GCV 10 terkecil pertama","\n")
print(dataG[1:10,])
#estimasi parameter pada titik knot optimal
mingcv=dataG[1,1]
knotgcv=knot1[dataG[1,2],]
datagcv1=matrix(ncol=m,nrow=N)
for (j in 1:m)
{
datagcv1[,j]=pmax(data[,(j+para+1)]-knotgcv[j],0)
}
mxgcv=as.matrix(cbind(aa,data2,datagcv1))
C=pinv(t(mxgcv)%*%mxgcv)
B=C%*%(t(mxgcv)%*%data[,1])
rownames(B)=c("b0",paste0("x",1:m),paste0("(x",1:m,"-k1)+"))
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 1","\n")
cat("==============================================","\n")
print(B)
cat("\n")
invisible(list(knot=matrix(knotgcv,nrow=1),mingcv=mingcv,B=B))
}
hasil1=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 k1_x1 k1_x2 k1_x3 k1_x4 k1_x5 k1_x6
## 60.891404 15.000000 13.823061 30.495918 4.788163 26.577551 62.518980 15.367347
## Nilai GCV 10 terkecil pertama
## GCV knot_ke k1_x1 k1_x2 k1_x3 k1_x4 k1_x5 k1_x6
## [1,] 60.89140 15 13.82306 30.49592 4.788163 26.57755 62.51898 15.36735
## [2,] 60.90658 16 14.59327 32.52898 5.022041 28.34939 62.97224 16.39184
## [3,] 60.93708 17 15.36347 34.56204 5.255918 30.12122 63.42551 17.41633
## [4,] 61.00914 18 16.13367 36.59510 5.489796 31.89306 63.87878 18.44082
## [5,] 61.03068 19 16.90388 38.62816 5.723673 33.66490 64.33204 19.46531
## [6,] 61.03475 14 13.05286 28.46286 4.554286 24.80571 62.06571 14.34286
## [7,] 61.08612 20 17.67408 40.66122 5.957551 35.43673 64.78531 20.48980
## [8,] 61.22630 21 18.44429 42.69429 6.191429 37.20857 65.23857 21.51429
## [9,] 61.33258 22 19.21449 44.72735 6.425306 38.98041 65.69184 22.53878
## [10,] 61.42555 13 12.28265 26.42980 4.320408 23.03388 61.61245 13.31837
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 1
## ==============================================
## [,1]
## b0 55.235597520
## x1 -0.060206537
## x2 -0.139780329
## x3 4.602967603
## x4 -0.613837239
## x5 0.008921615
## x6 0.059789517
## (x1-k1)+ -0.866308131
## (x2-k1)+ -0.081699055
## (x3-k1)+ -5.057603076
## (x4-k1)+ 0.720183134
## (x5-k1)+ 1.020919184
## (x6-k1)+ -0.129832245
#=========================#
#satu knot Uji signifikansi#
#=========================#
library(pracma)
uji=function(alpha,para)
{
knot_raw <- read.csv(paste0(folder,"dataAll knot 1.csv"), header = FALSE, skip = 1)
colnames(knot_raw) <- c("no","GCV","Rsq","knot_ke",paste0("k",1:6))
knotgcv <- knot_raw[which.min(knot_raw$GCV), paste0("k",1:6)]
knot <- matrix(as.numeric(knotgcv), nrow = 1)
data=as.matrix(data)
knot=as.matrix(knot)
ybar=mean(data[,1])
n=nrow(data)
m=ncol(data)-para-1 #banyaknya variabel prediktor
k=1 #banyaknya titik knot
dataA=data[,rep((para+2):ncol(data),each=k)]
dataA=as.matrix(dataA)
satu=rep(1,n)
n1=ncol(knot)
data.knot=matrix(ncol=n1,nrow=n)
for (i in 1:n1)
{
data.knot[,i]=pmax(dataA[,i]-knot[1,i],0)
}
mx=satu
for (j in 1:m)
{
mx=cbind(mx,data[,(para+1+j)],data.knot[,((j-1)*k+1):(j*k)])
}
mx=as.matrix(mx)
B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1]
rownames(B)=c("b0",unlist(lapply(1:m,function(j) c(paste0("x",j),paste0("(x",j,"-k",1:k,")+")))))
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')
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('Parameter',rownames(B)[i],': Tolak Ho, signifikan, pvalue =',pval[i],'\n') else
cat('Parameter',rownames(B)[i],': Gagal tolak Ho, tidak signifikan, pvalue =',pval[i],'\n')
}
tabel=cbind(Estimasi=B[,1],SE=SE,thit=thit,pvalue=pval)
cat('=============================================','\n')
cat('Estimasi parameter, nilai t hitung, dan p-value','\n')
cat('=============================================','\n')
print(tabel)
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=paste0(folder,"output uji residual knot1.csv"))
write.csv(mx,file=paste0(folder,"output uji mx knot1.csv"))
write.csv(yhat,file=paste0(folder,"output uji yhat knot1.csv"))
invisible(list(B=B,yhat=yhat,SSE=SSE,SSR=SSR,SST=SST,Rsq=Rsq,Fhit=Fhit,pvalue=pvalue))
}
hasil_uji1=uji(0.05,0)
## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang
## signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Parameter b0 : Gagal tolak Ho, tidak signifikan, pvalue = 0.3883954
## Parameter x1 : Gagal tolak Ho, tidak signifikan, pvalue = 0.6170333
## Parameter (x1-k1)+ : Tolak Ho, signifikan, pvalue = 6.567427e-06
## Parameter x2 : Tolak Ho, signifikan, pvalue = 0.003580801
## Parameter (x2-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.223957
## Parameter x3 : Tolak Ho, signifikan, pvalue = 0.005646631
## Parameter (x3-k1)+ : Tolak Ho, signifikan, pvalue = 0.003879527
## Parameter x4 : Tolak Ho, signifikan, pvalue = 5.156252e-13
## Parameter (x4-k1)+ : Tolak Ho, signifikan, pvalue = 1.687891e-05
## Parameter x5 : Gagal tolak Ho, tidak signifikan, pvalue = 0.9927584
## Parameter (x5-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3135006
## Parameter x6 : Gagal tolak Ho, tidak signifikan, pvalue = 0.6565236
## Parameter (x6-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.4467439
## =============================================
## Estimasi parameter, nilai t hitung, dan p-value
## =============================================
## Estimasi SE thit pvalue
## b0 55.235597557 63.98324665 0.863282194 3.883954e-01
## x1 -0.060206537 0.12032368 -0.500371464 6.170333e-01
## (x1-k1)+ -0.866308131 0.19016176 -4.555637852 6.567427e-06
## x2 -0.139780329 0.04775948 -2.926755428 3.580801e-03
## (x2-k1)+ -0.081699055 0.06709957 -1.217579486 2.239570e-01
## x3 4.602967603 1.65596221 2.779633244 5.646631e-03
## (x3-k1)+ -5.057603076 1.74321280 -2.901311348 3.879527e-03
## x4 -0.613837239 0.08276072 -7.417012199 5.156252e-13
## (x4-k1)+ 0.720183134 0.16575456 4.344876881 1.687891e-05
## x5 0.008921615 0.98248945 0.009080621 9.927584e-01
## (x5-k1)+ 1.020919185 1.01189430 1.008918804 3.135006e-01
## x6 0.059789517 0.13436334 0.444983846 6.565236e-01
## (x6-k1)+ -0.129832245 0.17050559 -0.761454471 4.467439e-01
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 12 77085.11 6423.759 108.2327
## Error 501 29735.03 59.35135
## Total 513 106820.1
## =============================================
## s= 7.703983 Rsq= 72.16346
## pvalue(F)= 1.293317e-130
#################
#DUA TITIK KNOT#
#################
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)
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)]
colnames(knot2)=paste0("k",rep(1:2,m),"_x",rep(1:m,each=2))
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))
{
b=ceiling(j/2)
data1[,j]=pmax(data[,(b+para+1)]-knot2[i,j],0)
}
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
A2=((N-sum((mx%*%C)*mx))/N)^2 #sama dengan (sum(diag(I-A))/N)^2
GCV[i]=MSE[i]/A2
}
dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot2))
dataG=dataAll[order(GCV),-2]
write.csv(dataAll,file=paste0(folder,"dataAll knot 2.csv"))
cat("==============================================","\n")
cat("HASIL GCV terkecil dengan 2 knot","\n")
cat("==============================================","\n")
print(dataG[1,1:(2+2*m)])
cat("Nilai GCV 10 terkecil pertama","\n")
print(dataG[1:10,])
#estimasi parameter pada titik knot optimal
mingcv=dataG[1,1]
knotgcv=knot2[dataG[1,2],]
datagcv1=matrix(ncol=2*m,nrow=N)
for (j in 1:(2*m))
{
b=ceiling(j/2)
datagcv1[,j]=pmax(data[,(b+para+1)]-knotgcv[j],0)
}
mxgcv=as.matrix(cbind(aa,data2,datagcv1))
C=pinv(t(mxgcv)%*%mxgcv)
B=C%*%(t(mxgcv)%*%data[,1])
rownames(B)=c("b0",paste0("x",1:m),paste0("(x",rep(1:m,each=2),"-k",rep(1:2,m),")+"))
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 2","\n")
cat("==============================================","\n")
print(B)
cat("\n")
invisible(list(knot=matrix(knotgcv,nrow=1),mingcv=mingcv,B=B))
}
hasil2=GCV2(data)
## ==============================================
## HASIL GCV terkecil dengan 2 knot
## ==============================================
## GCV knot_ke k1_x1 k2_x1 k1_x2 k2_x2 k1_x3
## 59.653741 245.000000 6.891224 22.295306 12.198367 52.859592 2.683265
## k2_x3 k1_x4 k2_x4 k1_x5 k2_x5 k1_x6 k2_x6
## 7.360816 10.631020 46.067755 58.439592 67.504898 6.146939 26.636735
## Nilai GCV 10 terkecil pertama
## GCV knot_ke k1_x1 k2_x1 k1_x2 k2_x2 k1_x3 k2_x3
## [1,] 59.65374 245 6.891224 22.29531 12.19837 52.85959 2.683265 7.360816
## [2,] 59.71889 239 6.891224 17.67408 12.19837 40.66122 2.683265 5.957551
## [3,] 59.71968 246 6.891224 23.06551 12.19837 54.89265 2.683265 7.594694
## [4,] 59.72621 244 6.891224 21.52510 12.19837 50.82653 2.683265 7.126939
## [5,] 59.72924 238 6.891224 16.90388 12.19837 38.62816 2.683265 5.723673
## [6,] 59.73756 905 23.065510 29.22714 54.89265 71.15714 7.594694 9.465714
## [7,] 59.76477 237 6.891224 16.13367 12.19837 36.59510 2.683265 5.489796
## [8,] 59.76563 241 6.891224 19.21449 12.19837 44.72735 2.683265 6.425306
## [9,] 59.76975 240 6.891224 18.44429 12.19837 42.69429 2.683265 6.191429
## [10,] 59.77079 242 6.891224 19.98469 12.19837 46.76041 2.683265 6.659184
## k1_x4 k2_x4 k1_x5 k2_x5 k1_x6 k2_x6
## [1,] 10.63102 46.06776 58.43959 67.50490 6.146939 26.63673
## [2,] 10.63102 35.43673 58.43959 64.78531 6.146939 20.48980
## [3,] 10.63102 47.83959 58.43959 67.95816 6.146939 27.66122
## [4,] 10.63102 44.29592 58.43959 67.05163 6.146939 25.61224
## [5,] 10.63102 33.66490 58.43959 64.33204 6.146939 19.46531
## [6,] 47.83959 62.01429 67.95816 71.58429 27.661224 35.85714
## [7,] 10.63102 31.89306 58.43959 63.87878 6.146939 18.44082
## [8,] 10.63102 38.98041 58.43959 65.69184 6.146939 22.53878
## [9,] 10.63102 37.20857 58.43959 65.23857 6.146939 21.51429
## [10,] 10.63102 40.75224 58.43959 66.14510 6.146939 23.56327
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 2
## ==============================================
## [,1]
## b0 232.61514834
## x1 -0.21284620
## x2 -0.47339109
## x3 -3.37097474
## x4 -0.93861870
## x5 -2.77940374
## x6 0.78647431
## (x1-k1)+ -0.04625130
## (x1-k2)+ -1.12800600
## (x2-k1)+ 0.33226449
## (x2-k2)+ -0.04946566
## (x3-k1)+ 6.15885749
## (x3-k2)+ -4.00900376
## (x4-k1)+ 0.87768543
## (x4-k2)+ 0.02702941
## (x5-k1)+ 4.11508131
## (x5-k2)+ -0.43791153
## (x6-k1)+ -0.91535404
## (x6-k2)+ 0.05822887
#=========================#
#dua knot Uji signifikansi#
#=========================#
library(pracma)
uji=function(alpha,para)
{
knot_raw <- read.csv(paste0(folder,"dataAll knot 2.csv"), header = FALSE, skip = 1)
nama_knot <- paste0("k",rep(1:6,each=2),rep(c("a","b"),6))
colnames(knot_raw) <- c("no","GCV","Rsq","knot_ke",nama_knot)
knotgcv <- knot_raw[which.min(knot_raw$GCV), nama_knot]
knot <- matrix(as.numeric(knotgcv), nrow = 1)
data=as.matrix(data)
knot=as.matrix(knot)
ybar=mean(data[,1])
n=nrow(data)
m=ncol(data)-para-1 #banyaknya variabel prediktor
k=2 #banyaknya titik knot
dataA=data[,rep((para+2):ncol(data),each=k)]
dataA=as.matrix(dataA)
satu=rep(1,n)
n1=ncol(knot)
data.knot=matrix(ncol=n1,nrow=n)
for (i in 1:n1)
{
data.knot[,i]=pmax(dataA[,i]-knot[1,i],0)
}
mx=satu
for (j in 1:m)
{
mx=cbind(mx,data[,(para+1+j)],data.knot[,((j-1)*k+1):(j*k)])
}
mx=as.matrix(mx)
B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1]
rownames(B)=c("b0",unlist(lapply(1:m,function(j) c(paste0("x",j),paste0("(x",j,"-k",1:k,")+")))))
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')
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('Parameter',rownames(B)[i],': Tolak Ho, signifikan, pvalue =',pval[i],'\n') else
cat('Parameter',rownames(B)[i],': Gagal tolak Ho, tidak signifikan, pvalue =',pval[i],'\n')
}
tabel=cbind(Estimasi=B[,1],SE=SE,thit=thit,pvalue=pval)
cat('=============================================','\n')
cat('Estimasi parameter, nilai t hitung, dan p-value','\n')
cat('=============================================','\n')
print(tabel)
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=paste0(folder,"output uji residual knot2.csv"))
write.csv(mx,file=paste0(folder,"output uji mx knot2.csv"))
write.csv(yhat,file=paste0(folder,"output uji yhat knot2.csv"))
invisible(list(B=B,yhat=yhat,SSE=SSE,SSR=SSR,SST=SST,Rsq=Rsq,Fhit=Fhit,pvalue=pvalue))
}
hasil_uji2=uji(0.05,0)
## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang
## signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Parameter b0 : Gagal tolak Ho, tidak signifikan, pvalue = 0.2033832
## Parameter x1 : Gagal tolak Ho, tidak signifikan, pvalue = 0.5802154
## Parameter (x1-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.9150597
## Parameter (x1-k2)+ : Tolak Ho, signifikan, pvalue = 0.0002143497
## Parameter x2 : Tolak Ho, signifikan, pvalue = 0.001689131
## Parameter (x2-k1)+ : Tolak Ho, signifikan, pvalue = 0.04228767
## Parameter (x2-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.580082
## Parameter x3 : Gagal tolak Ho, tidak signifikan, pvalue = 0.53971
## Parameter (x3-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3171688
## Parameter (x3-k2)+ : Tolak Ho, signifikan, pvalue = 0.002221382
## Parameter x4 : Tolak Ho, signifikan, pvalue = 4.68481e-11
## Parameter (x4-k1)+ : Tolak Ho, signifikan, pvalue = 7.791085e-06
## Parameter (x4-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.9171649
## Parameter x5 : Gagal tolak Ho, tidak signifikan, pvalue = 0.3726539
## Parameter (x5-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2028383
## Parameter (x5-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.304175
## Parameter x6 : Gagal tolak Ho, tidak signifikan, pvalue = 0.1014708
## Parameter (x6-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.0736619
## Parameter (x6-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.7026649
## =============================================
## Estimasi parameter, nilai t hitung, dan p-value
## =============================================
## Estimasi SE thit pvalue
## b0 232.61514829 182.63580167 1.2736558 2.033832e-01
## x1 -0.21284620 0.38459113 -0.5534350 5.802154e-01
## (x1-k1)+ -0.04625130 0.43341679 -0.1067132 9.150597e-01
## (x1-k2)+ -1.12800600 0.30248414 -3.7291410 2.143497e-04
## x2 -0.47339109 0.14993211 -3.1573696 1.689131e-03
## (x2-k1)+ 0.33226449 0.16319906 2.0359461 4.228767e-02
## (x2-k2)+ -0.04946566 0.08934786 -0.5536301 5.800820e-01
## x3 -3.37097474 5.49306747 -0.6136780 5.397100e-01
## (x3-k1)+ 6.15885749 6.15083282 1.0013046 3.171688e-01
## (x3-k2)+ -4.00900376 1.30374304 -3.0749953 2.221382e-03
## x4 -0.93861870 0.13945275 -6.7307292 4.684810e-11
## (x4-k1)+ 0.87768543 0.19423483 4.5186820 7.791085e-06
## (x4-k2)+ 0.02702941 0.25975154 0.1040587 9.171649e-01
## x5 -2.77940374 3.11479854 -0.8923222 3.726539e-01
## (x5-k1)+ 4.11508131 3.22702037 1.2751953 2.028383e-01
## (x5-k2)+ -0.43791153 0.42573989 -1.0285894 3.041750e-01
## x6 0.78647431 0.47931898 1.6408161 1.014708e-01
## (x6-k1)+ -0.91535404 0.51065506 -1.7925095 7.366190e-02
## (x6-k2)+ 0.05822887 0.15245290 0.3819466 7.026649e-01
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 18 78383.06 4354.614 75.80013
## Error 495 28437.08 57.44864
## Total 513 106820.1
## =============================================
## s= 7.579488 Rsq= 73.37854
## pvalue(F)= 1.905753e-129
##################
#TIGA TITIK KNOT#
##################
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)
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)]
colnames(knot2)=paste0("k",rep(1:3,m),"_x",rep(1:m,each=3))
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)
data1[,j]=pmax(data[,(b+para+1)]-knot2[i,j],0)
}
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
A2=((N-sum((mx%*%C)*mx))/N)^2 #sama dengan (sum(diag(I-A))/N)^2
GCV[i]=MSE[i]/A2
}
dataAll=as.matrix(cbind(GCV,Rsq,knotke,knot2))
dataG=dataAll[order(GCV),-2]
write.csv(dataAll,file=paste0(folder,"dataAll knot 3.csv"))
cat("==============================================","\n")
cat("HASIL GCV terkecil dengan 3 knot","\n")
cat("==============================================","\n")
print(dataG[1,1:(2+3*m)])
cat("Nilai GCV 10 terkecil pertama","\n")
print(dataG[1:10,])
#estimasi parameter pada titik knot optimal
mingcv=dataG[1,1]
knotgcv=knot2[dataG[1,2],]
datagcv1=matrix(ncol=3*m,nrow=N)
for (j in 1:(3*m))
{
b=ceiling(j/3)
datagcv1[,j]=pmax(data[,(b+para+1)]-knotgcv[j],0)
}
mxgcv=as.matrix(cbind(aa,data2,datagcv1))
C=pinv(t(mxgcv)%*%mxgcv)
B=C%*%(t(mxgcv)%*%data[,1])
rownames(B)=c("b0",paste0("x",1:m),paste0("(x",rep(1:m,each=3),"-k",rep(1:3,m),")+"))
cat("\n")
cat("==============================================","\n")
cat("HASIL ESTIMASI PARAMETER TITIK KNOT KE 3","\n")
cat("==============================================","\n")
print(B)
cat("\n")
invisible(list(knot=matrix(knotgcv,nrow=1),mingcv=mingcv,B=B))
}
hasil3=GCV3(data)
## ==============================================
## HASIL GCV terkecil dengan 3 knot
## ==============================================
## GCV knot_ke k1_x1 k2_x1 k3_x1 k1_x2
## 57.871554 11127.000000 13.052857 24.605918 28.456939 28.462857
## k2_x2 k3_x2 k1_x3 k2_x3 k3_x3 k1_x4
## 58.958776 69.124082 4.554286 8.062449 9.231837 24.805714
## k2_x4 k3_x4 k1_x5 k2_x5 k3_x5 k1_x6
## 51.383265 60.242449 62.065714 68.864694 71.131020 14.342857
## k2_x6 k3_x6
## 29.710204 34.832653
## Nilai GCV 10 terkecil pertama
## GCV knot_ke k1_x1 k2_x1 k3_x1 k1_x2 k2_x2 k3_x2
## [1,] 57.87155 11127 13.05286 24.60592 28.45694 28.46286 58.95878 69.12408
## [2,] 57.90291 11144 13.05286 25.37612 27.68673 28.46286 60.99184 67.09102
## [3,] 57.91615 11126 13.05286 24.60592 27.68673 28.46286 58.95878 67.09102
## [4,] 57.97474 11655 13.82306 24.60592 28.45694 30.49592 58.95878 69.12408
## [5,] 58.00390 11672 13.82306 25.37612 27.68673 30.49592 60.99184 67.09102
## [6,] 58.01200 11654 13.82306 24.60592 27.68673 30.49592 58.95878 67.09102
## [7,] 58.01834 10566 12.28265 24.60592 28.45694 26.42980 58.95878 69.12408
## [8,] 58.05155 11145 13.05286 25.37612 28.45694 28.46286 60.99184 69.12408
## [9,] 58.07181 10583 12.28265 25.37612 27.68673 26.42980 60.99184 67.09102
## [10,] 58.08304 10565 12.28265 24.60592 27.68673 26.42980 58.95878 67.09102
## k1_x3 k2_x3 k3_x3 k1_x4 k2_x4 k3_x4 k1_x5 k2_x5
## [1,] 4.554286 8.062449 9.231837 24.80571 51.38327 60.24245 62.06571 68.86469
## [2,] 4.554286 8.296327 8.997959 24.80571 53.15510 58.47061 62.06571 69.31796
## [3,] 4.554286 8.062449 8.997959 24.80571 51.38327 58.47061 62.06571 68.86469
## [4,] 4.788163 8.062449 9.231837 26.57755 51.38327 60.24245 62.51898 68.86469
## [5,] 4.788163 8.296327 8.997959 26.57755 53.15510 58.47061 62.51898 69.31796
## [6,] 4.788163 8.062449 8.997959 26.57755 51.38327 58.47061 62.51898 68.86469
## [7,] 4.320408 8.062449 9.231837 23.03388 51.38327 60.24245 61.61245 68.86469
## [8,] 4.554286 8.296327 9.231837 24.80571 53.15510 60.24245 62.06571 69.31796
## [9,] 4.320408 8.296327 8.997959 23.03388 53.15510 58.47061 61.61245 69.31796
## [10,] 4.320408 8.062449 8.997959 23.03388 51.38327 58.47061 61.61245 68.86469
## k3_x5 k1_x6 k2_x6 k3_x6
## [1,] 71.13102 14.34286 29.71020 34.83265
## [2,] 70.67776 14.34286 30.73469 33.80816
## [3,] 70.67776 14.34286 29.71020 33.80816
## [4,] 71.13102 15.36735 29.71020 34.83265
## [5,] 70.67776 15.36735 30.73469 33.80816
## [6,] 70.67776 15.36735 29.71020 33.80816
## [7,] 71.13102 13.31837 29.71020 34.83265
## [8,] 71.13102 14.34286 30.73469 34.83265
## [9,] 70.67776 13.31837 30.73469 33.80816
## [10,] 70.67776 13.31837 29.71020 33.80816
##
## ==============================================
## HASIL ESTIMASI PARAMETER TITIK KNOT KE 3
## ==============================================
## [,1]
## b0 167.379164922
## x1 0.006721307
## x2 -0.162048384
## x3 3.088515434
## x4 -0.572073896
## x5 -1.828642791
## x6 0.121548312
## (x1-k1)+ -0.884363855
## (x1-k2)+ 0.860349570
## (x1-k3)+ -2.205508264
## (x2-k1)+ -0.027020113
## (x2-k2)+ -0.313063519
## (x2-k3)+ 0.513304370
## (x3-k1)+ -1.193279344
## (x3-k2)+ -6.518365931
## (x3-k3)+ 6.397219439
## (x4-k1)+ 0.756445243
## (x4-k2)+ 0.935416563
## (x4-k3)+ -2.197128367
## (x5-k1)+ 3.616134314
## (x5-k2)+ -2.564442076
## (x5-k3)+ 2.261039366
## (x6-k1)+ -0.235273051
## (x6-k2)+ -0.185090874
## (x6-k3)+ 0.665478503
#=========================#
#tiga knot Uji signifikansi#
#=========================#
library(pracma)
uji=function(alpha,para)
{
knot_raw <- read.csv(paste0(folder,"dataAll knot 3.csv"), header = FALSE, skip = 1)
nama_knot <- paste0("k",rep(1:6,each=3),rep(c("a","b","c"),6))
colnames(knot_raw) <- c("no","GCV","Rsq","knot_ke",nama_knot)
knotgcv <- knot_raw[which.min(knot_raw$GCV), nama_knot]
knot <- matrix(as.numeric(knotgcv), nrow = 1)
data=as.matrix(data)
knot=as.matrix(knot)
ybar=mean(data[,1])
n=nrow(data)
m=ncol(data)-para-1 #banyaknya variabel prediktor
k=3 #banyaknya titik knot
dataA=data[,rep((para+2):ncol(data),each=k)]
dataA=as.matrix(dataA)
satu=rep(1,n)
n1=ncol(knot)
data.knot=matrix(ncol=n1,nrow=n)
for (i in 1:n1)
{
data.knot[,i]=pmax(dataA[,i]-knot[1,i],0)
}
mx=satu
for (j in 1:m)
{
mx=cbind(mx,data[,(para+1+j)],data.knot[,((j-1)*k+1):(j*k)])
}
mx=as.matrix(mx)
B=(pinv(t(mx)%*%mx))%*%t(mx)%*%data[,1]
rownames(B)=c("b0",unlist(lapply(1:m,function(j) c(paste0("x",j),paste0("(x",j,"-k",1:k,")+")))))
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')
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('Parameter',rownames(B)[i],': Tolak Ho, signifikan, pvalue =',pval[i],'\n') else
cat('Parameter',rownames(B)[i],': Gagal tolak Ho, tidak signifikan, pvalue =',pval[i],'\n')
}
tabel=cbind(Estimasi=B[,1],SE=SE,thit=thit,pvalue=pval)
cat('=============================================','\n')
cat('Estimasi parameter, nilai t hitung, dan p-value','\n')
cat('=============================================','\n')
print(tabel)
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=paste0(folder,"output uji residual knot3.csv"))
write.csv(mx,file=paste0(folder,"output uji mx knot3.csv"))
write.csv(yhat,file=paste0(folder,"output uji yhat knot3.csv"))
invisible(list(B=B,yhat=yhat,SSE=SSE,SSR=SSR,SST=SST,Rsq=Rsq,Fhit=Fhit,pvalue=pvalue))
}
hasil_uji3=uji(0.05,0)
## ---------------------------------------
## Kesimpulan hasil uji simultan
## ---------------------------------------
## Tolak Ho yakni minimal terdapat 1 variabel bebas yang
## signifikan
##
## ---------------------------------------------
## Kesimpulan hasil uji parsial
## ---------------------------------------------
## Parameter b0 : Tolak Ho, signifikan, pvalue = 0.030933
## Parameter x1 : Gagal tolak Ho, tidak signifikan, pvalue = 0.9589731
## Parameter (x1-k1)+ : Tolak Ho, signifikan, pvalue = 0.0008273402
## Parameter (x1-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3383506
## Parameter (x1-k3)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.05676414
## Parameter x2 : Tolak Ho, signifikan, pvalue = 0.002971664
## Parameter (x2-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.7632598
## Parameter (x2-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.3196921
## Parameter (x2-k3)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2179738
## Parameter x3 : Gagal tolak Ho, tidak signifikan, pvalue = 0.1575939
## Parameter (x3-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.6503021
## Parameter (x3-k2)+ : Tolak Ho, signifikan, pvalue = 0.0001412823
## Parameter (x3-k3)+ : Tolak Ho, signifikan, pvalue = 9.938082e-06
## Parameter x4 : Tolak Ho, signifikan, pvalue = 1.323051e-09
## Parameter (x4-k1)+ : Tolak Ho, signifikan, pvalue = 0.002019186
## Parameter (x4-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.4230511
## Parameter (x4-k3)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.1188693
## Parameter x5 : Gagal tolak Ho, tidak signifikan, pvalue = 0.129938
## Parameter (x5-k1)+ : Tolak Ho, signifikan, pvalue = 0.00774821
## Parameter (x5-k2)+ : Tolak Ho, signifikan, pvalue = 0.0005331614
## Parameter (x5-k3)+ : Tolak Ho, signifikan, pvalue = 0.001533134
## Parameter x6 : Gagal tolak Ho, tidak signifikan, pvalue = 0.4278882
## Parameter (x6-k1)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.2545297
## Parameter (x6-k2)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.6005636
## Parameter (x6-k3)+ : Gagal tolak Ho, tidak signifikan, pvalue = 0.1704717
## =============================================
## Estimasi parameter, nilai t hitung, dan p-value
## =============================================
## Estimasi SE thit pvalue
## b0 167.379165037 77.34026676 2.16419172 3.093300e-02
## x1 0.006721307 0.13059047 0.05146859 9.589731e-01
## (x1-k1)+ -0.884363855 0.26285790 -3.36441804 8.273402e-04
## (x1-k2)+ 0.860349570 0.89772315 0.95836847 3.383506e-01
## (x1-k3)+ -2.205508264 1.15493272 -1.90964221 5.676414e-02
## x2 -0.162048384 0.05427585 -2.98564456 2.971664e-03
## (x2-k1)+ -0.027020113 0.08965733 -0.30137092 7.632598e-01
## (x2-k2)+ -0.313063519 0.31428791 -0.99610423 3.196921e-01
## (x2-k3)+ 0.513304370 0.41612923 1.23352155 2.179738e-01
## x3 3.088515432 2.18211160 1.41537923 1.575939e-01
## (x3-k1)+ -1.193279342 2.63055870 -0.45362202 6.503021e-01
## (x3-k2)+ -6.518365931 1.69914752 -3.83625662 1.412823e-04
## (x3-k3)+ 6.397219439 1.43264913 4.46530789 9.938082e-06
## x4 -0.572073896 0.09251360 -6.18367326 1.323051e-09
## (x4-k1)+ 0.756445243 0.24369214 3.10410199 2.019186e-03
## (x4-k2)+ 0.935416563 1.16662812 0.80181212 4.230511e-01
## (x4-k3)+ -2.197128367 1.40636153 -1.56227849 1.188693e-01
## x5 -1.828642793 1.20551281 -1.51690034 1.299380e-01
## (x5-k1)+ 3.616134316 1.35237864 2.67390671 7.748210e-03
## (x5-k2)+ -2.564442076 0.73549926 -3.48666848 5.331614e-04
## (x5-k3)+ 2.261039366 0.70961060 3.18631001 1.533134e-03
## x6 0.121548312 0.15318553 0.79347125 4.278882e-01
## (x6-k1)+ -0.235273051 0.20624303 -1.14075640 2.545297e-01
## (x6-k2)+ -0.185090874 0.35327308 -0.52393144 6.005636e-01
## (x6-k3)+ 0.665478503 0.48479173 1.37271011 1.704717e-01
## Analysis of Variance
## =============================================
## Sumber df SS MS Fhit
## Regresi 24 79897.37 3329.057 60.46587
## Error 489 26922.77 55.05679
## Total 513 106820.1
## =============================================
## s= 7.420026 Rsq= 74.79617
## pvalue(F)= 1.15666e-129