PSD 8 & 9
PERTEMUAN 8
Untuk tugas bagian pertemuan 8 ini saya akan melakukan interpretasi pada data bangkitan set.seed 1401201089 saja, sedangkan pada data triceps tidak saya interpretasikan karena isinya kurang lebih sama.
Library
library("MultiKink")# DATA TRICEPS
library(tidyverse)
library(ggplot2)
library(dplyr)
library(purrr)
library(rsample)Data
Data yang digunakan yaitu data yang dibangkitkan dari fungsi set.seed 1401201089 (NIM saya). Dibangkitkan 1000 data yang menyebar normal dengan rata-rata dan ragam sama dengan 1. Model yang digunakan pada data tersebut yaitu Model Polynomial Ordo 2 dengan persamaan sebagai berikut: \[y=5+2(data.x)+3(data.x)^2+error\]
set.seed(1401201089)
data.x <- rnorm(1000,1,1)
err <- rnorm(1000)
y <- 5+2*data.x+3*data.x^2+err ## Model Polynomial Ordo 2
plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70))
Dari grafik di atas, dapat diketahui bahwa data tidak membentuk pola yang linier karena dilihat dari bentuk grafiknya yang cenderung melengkung dan tidak membentuk garis lurus (linier).
Regresi Linier
regresi linear merupakan suatu pendekatan untuk mengetahui hubungan antara satu atau lebih variabel dependen dan juga variabel independen.
lin.mod <-lm(y~data.x)
plot(data.x,y,xlim=c( -2,5), ylim=c( -10,70))
lines(data.x,lin.mod$fitted.values,col="red")
Dari grafik di atas, dapat diketahui bahwa garis regresi linier tidak mendekati pola data, sehingga tidak akan tepat apabila digunakan fungsi liner biasa pada data ini.
summary(lin.mod)##
## Call:
## lm(formula = y ~ data.x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.228 -2.561 -1.441 1.091 29.837
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7871 0.1924 24.88 <2e-16 ***
## data.x 8.1205 0.1362 59.62 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.255 on 998 degrees of freedom
## Multiple R-squared: 0.7808, Adjusted R-squared: 0.7806
## F-statistic: 3555 on 1 and 998 DF, p-value: < 2.2e-16Polynomial
pol.mod <- lm(y~data.x+I(data.x^2)) #ordo 2
ix <- sort(data.x,index.return=T)$ix
plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70))
lines(data.x[ix], pol.mod$fitted.values[ix],col="blue", cex=2)
Garis biru berupakan hasil fungsi polynomial orde 2. Dapat dilihat bahwa
garis tersebut mendekati observasi yang ada, sehingga errornya menjadi
lebih kecil daripada fungsi regresi linier di atas.
summary(pol.mod)##
## Call:
## lm(formula = y ~ data.x + I(data.x^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3268 -0.6134 -0.0023 0.6302 2.7398
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.98088 0.04339 114.78 <2e-16 ***
## data.x 2.03376 0.05413 37.57 <2e-16 ***
## I(data.x^2) 2.98285 0.02185 136.54 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9591 on 997 degrees of freedom
## Multiple R-squared: 0.9889, Adjusted R-squared: 0.9888
## F-statistic: 4.43e+04 on 2 and 997 DF, p-value: < 2.2e-16Fungsi Tangga
#regresi fungsi tangga
range(data.x)## [1] -2.269277 3.974688data.x memiliki nilai minimum -2.269277 dan nilai maksimum 3.974688
c1 <- as.factor(ifelse(data.x<=0,1,0))
c2 <- as.factor(ifelse(data.x<=2 & data.x>0,1,0))
c3 <- as.factor(ifelse(data.x>2,1,0))
data.frame(y,c1,c2,c3)## y c1 c2 c3
## 1 39.499515 0 0 1
## 2 12.749125 0 1 0
## 3 6.016863 0 1 0
## 4 20.909510 0 1 0
## 5 12.983689 0 1 0
## 6 3.837443 0 1 0
## 7 6.145309 0 1 0
## 8 34.945835 0 0 1
## 9 11.216720 0 1 0
## 10 5.270016 0 1 0
## 11 7.576125 0 1 0
## 12 12.060699 0 1 0
## 13 6.101816 0 1 0
## 14 19.192795 0 1 0
## 15 5.587806 1 0 0
## 16 15.059727 0 1 0
## 17 6.276712 0 1 0
## 18 31.777766 0 0 1
## 19 9.415421 0 1 0
## 20 12.080452 0 1 0
## 21 13.890919 0 1 0
## 22 32.211739 0 0 1
## 23 15.610584 0 1 0
## 24 14.871444 0 1 0
## 25 49.832872 0 0 1
## 26 10.040268 0 1 0
## 27 3.677834 0 1 0
## 28 17.323060 0 1 0
## 29 6.695144 0 1 0
## 30 4.788230 1 0 0
## 31 31.482923 0 0 1
## 32 9.732201 0 1 0
## 33 7.345214 0 1 0
## 34 7.326607 0 1 0
## 35 27.964237 0 0 1
## 36 23.830006 0 0 1
## 37 4.989869 1 0 0
## 38 32.066483 0 0 1
## 39 8.181434 0 1 0
## 40 9.868413 0 1 0
## 41 19.459859 0 1 0
## 42 19.013075 0 1 0
## 43 3.984603 1 0 0
## 44 19.541820 0 1 0
## 45 8.081965 0 1 0
## 46 12.378170 0 1 0
## 47 3.661592 1 0 0
## 48 9.478298 0 1 0
## 49 11.865078 0 1 0
## 50 18.454943 0 1 0
## 51 10.550340 0 1 0
## 52 12.948436 0 1 0
## 53 24.624490 0 0 1
## 54 5.388519 0 1 0
## 55 5.701085 1 0 0
## 56 6.831478 1 0 0
## 57 5.999681 0 1 0
## 58 31.529093 0 0 1
## 59 33.058633 0 0 1
## 60 4.734353 1 0 0
## 61 7.502639 0 1 0
## 62 9.505205 0 1 0
## 63 7.468565 0 1 0
## 64 35.019255 0 0 1
## 65 15.542595 0 1 0
## 66 13.610536 0 1 0
## 67 17.459436 0 1 0
## 68 8.998782 0 1 0
## 69 5.305957 0 1 0
## 70 21.818425 0 0 1
## 71 4.560831 0 1 0
## 72 9.229415 0 1 0
## 73 32.524910 0 0 1
## 74 18.983873 0 1 0
## 75 9.153255 0 1 0
## 76 4.600770 1 0 0
## 77 8.579895 0 1 0
## 78 8.138211 1 0 0
## 79 15.741459 0 1 0
## 80 4.172438 0 1 0
## 81 13.258199 0 1 0
## 82 4.734818 1 0 0
## 83 7.448220 0 1 0
## 84 16.180029 0 1 0
## 85 7.602129 0 1 0
## 86 8.656882 0 1 0
## 87 35.337800 0 0 1
## 88 16.747389 0 1 0
## 89 5.125901 0 1 0
## 90 54.752018 0 0 1
## 91 21.878484 0 0 1
## 92 8.779497 0 1 0
## 93 6.010760 0 1 0
## 94 4.918885 0 1 0
## 95 11.061403 0 1 0
## 96 4.530402 1 0 0
## 97 5.599072 0 1 0
## 98 11.889344 0 1 0
## 99 24.991938 0 0 1
## 100 14.523728 0 1 0
## 101 22.653474 0 0 1
## 102 25.432802 0 0 1
## 103 12.032475 0 1 0
## 104 8.735108 0 1 0
## 105 4.936549 0 1 0
## 106 5.605776 1 0 0
## 107 21.893566 0 0 1
## 108 7.003312 0 1 0
## 109 12.938143 0 1 0
## 110 18.725358 0 1 0
## 111 13.069686 0 1 0
## 112 47.030774 0 0 1
## 113 11.770828 0 1 0
## 114 16.870363 0 1 0
## 115 13.604770 0 1 0
## 116 33.385939 0 0 1
## 117 4.914885 0 1 0
## 118 14.295488 0 1 0
## 119 9.753226 0 1 0
## 120 19.546576 0 1 0
## 121 7.952468 0 1 0
## 122 6.386436 1 0 0
## 123 17.129103 0 1 0
## 124 9.903268 0 1 0
## 125 14.716121 0 1 0
## 126 8.574384 0 1 0
## 127 6.121186 1 0 0
## 128 12.628375 0 1 0
## 129 18.035160 0 1 0
## 130 10.836767 0 1 0
## 131 29.376726 0 0 1
## 132 4.173904 1 0 0
## 133 19.478166 0 1 0
## 134 21.489350 0 0 1
## 135 14.026918 0 1 0
## 136 10.119741 0 1 0
## 137 28.809429 0 0 1
## 138 59.125920 0 0 1
## 139 5.070682 1 0 0
## 140 5.477479 0 1 0
## 141 17.920956 0 1 0
## 142 8.825429 0 1 0
## 143 8.755277 0 1 0
## 144 18.651082 0 1 0
## 145 11.469977 0 1 0
## 146 10.309696 0 1 0
## 147 31.438219 0 0 1
## 148 6.413382 0 1 0
## 149 9.015866 0 1 0
## 150 12.339588 0 1 0
## 151 16.689514 0 1 0
## 152 3.243224 1 0 0
## 153 7.918566 0 1 0
## 154 14.002696 0 1 0
## 155 14.720993 0 1 0
## 156 28.499829 0 0 1
## 157 8.928185 0 1 0
## 158 17.354999 0 1 0
## 159 8.584380 0 1 0
## 160 9.370546 0 1 0
## 161 12.098260 0 1 0
## 162 18.298912 0 1 0
## 163 5.612740 0 1 0
## 164 6.281431 0 1 0
## 165 4.027428 1 0 0
## 166 17.724560 0 1 0
## 167 8.242046 0 1 0
## 168 15.781609 0 1 0
## 169 7.211849 1 0 0
## 170 6.134113 0 1 0
## 171 11.868578 0 1 0
## 172 16.737163 0 1 0
## 173 13.488745 0 1 0
## 174 5.567144 1 0 0
## 175 5.456703 0 1 0
## 176 7.859555 0 1 0
## 177 7.415340 0 1 0
## 178 4.120318 1 0 0
## 179 8.810289 0 1 0
## 180 30.991906 0 0 1
## 181 6.176979 0 1 0
## 182 7.274394 1 0 0
## 183 11.555162 0 1 0
## 184 6.749226 0 1 0
## 185 14.617773 0 1 0
## 186 9.230112 0 1 0
## 187 3.778021 1 0 0
## 188 15.624169 0 1 0
## 189 5.875087 0 1 0
## 190 4.054607 1 0 0
## 191 33.124518 0 0 1
## 192 5.458182 1 0 0
## 193 13.914535 0 1 0
## 194 20.751006 0 1 0
## 195 18.045710 0 1 0
## 196 6.430649 0 1 0
## 197 10.258353 0 1 0
## 198 5.988092 0 1 0
## 199 6.008455 0 1 0
## 200 9.028511 0 1 0
## 201 10.675589 0 1 0
## 202 4.487455 1 0 0
## 203 23.990202 0 0 1
## 204 8.529949 0 1 0
## 205 27.828660 0 0 1
## 206 5.325857 0 1 0
## 207 5.598870 0 1 0
## 208 6.131010 0 1 0
## 209 23.424828 0 0 1
## 210 12.019492 0 1 0
## 211 4.827963 0 1 0
## 212 29.317585 0 0 1
## 213 3.866884 1 0 0
## 214 12.397518 0 1 0
## 215 6.921884 1 0 0
## 216 19.478883 0 1 0
## 217 5.440721 1 0 0
## 218 5.666796 1 0 0
## 219 8.427803 0 1 0
## 220 18.368152 0 1 0
## 221 6.733518 0 1 0
## 222 4.906534 1 0 0
## 223 40.779652 0 0 1
## 224 11.325396 0 1 0
## 225 8.103497 0 1 0
## 226 17.820196 0 1 0
## 227 19.728230 0 0 1
## 228 18.199179 0 1 0
## 229 11.663342 0 1 0
## 230 6.706007 0 1 0
## 231 9.402071 0 1 0
## 232 6.091208 1 0 0
## 233 21.720572 0 0 1
## 234 3.857362 1 0 0
## 235 3.937332 0 1 0
## 236 8.463425 0 1 0
## 237 19.262016 0 1 0
## 238 6.191226 0 1 0
## 239 7.539181 0 1 0
## 240 4.938916 0 1 0
## 241 4.781119 0 1 0
## 242 4.973809 0 1 0
## 243 7.746348 1 0 0
## 244 40.306514 0 0 1
## 245 4.105065 1 0 0
## 246 5.155292 1 0 0
## 247 14.619954 0 1 0
## 248 12.491501 0 1 0
## 249 6.232480 0 1 0
## 250 10.609233 0 1 0
## 251 19.484647 0 1 0
## 252 5.959519 0 1 0
## 253 9.100788 0 1 0
## 254 7.279067 0 1 0
## 255 6.330971 1 0 0
## 256 26.248618 0 0 1
## 257 13.523359 0 1 0
## 258 14.859530 0 1 0
## 259 11.123749 0 1 0
## 260 4.949847 0 1 0
## 261 14.609809 0 1 0
## 262 12.982827 0 1 0
## 263 6.893877 0 1 0
## 264 25.598143 0 0 1
## 265 10.199922 0 1 0
## 266 6.633484 1 0 0
## 267 13.732128 0 1 0
## 268 17.508419 0 1 0
## 269 7.003404 0 1 0
## 270 9.396874 0 1 0
## 271 8.200239 0 1 0
## 272 15.112267 0 1 0
## 273 21.012794 0 1 0
## 274 11.588666 0 1 0
## 275 15.176189 0 1 0
## 276 5.433403 1 0 0
## 277 19.041608 0 1 0
## 278 5.687769 1 0 0
## 279 4.274140 1 0 0
## 280 12.742207 0 1 0
## 281 32.536189 0 0 1
## 282 8.393627 0 1 0
## 283 10.459198 0 1 0
## 284 18.090353 0 1 0
## 285 6.194114 0 1 0
## 286 39.940484 0 0 1
## 287 7.621189 0 1 0
## 288 9.102056 0 1 0
## 289 4.682196 0 1 0
## 290 14.390928 0 1 0
## 291 13.284461 0 1 0
## 292 22.346089 0 0 1
## 293 15.913774 0 1 0
## 294 12.131900 0 1 0
## 295 4.933509 1 0 0
## 296 4.040284 1 0 0
## 297 3.912930 1 0 0
## 298 7.096302 0 1 0
## 299 3.779101 0 1 0
## 300 12.335089 0 1 0
## 301 9.181421 0 1 0
## 302 31.392249 0 0 1
## 303 13.260807 0 1 0
## 304 25.704232 0 0 1
## 305 5.995122 1 0 0
## 306 5.171751 0 1 0
## 307 10.597258 0 1 0
## 308 12.377886 0 1 0
## 309 16.324356 0 1 0
## 310 16.233471 0 1 0
## 311 3.522343 0 1 0
## 312 4.264001 1 0 0
## 313 13.582024 0 1 0
## 314 8.702469 0 1 0
## 315 29.263366 0 0 1
## 316 14.536128 0 1 0
## 317 20.348052 0 0 1
## 318 12.163597 0 1 0
## 319 9.202677 0 1 0
## 320 31.811065 0 0 1
## 321 7.022685 1 0 0
## 322 16.993332 0 1 0
## 323 19.275811 0 1 0
## 324 11.677275 0 1 0
## 325 32.142225 0 0 1
## 326 21.091392 0 1 0
## 327 10.608925 0 1 0
## 328 5.442323 0 1 0
## 329 6.737721 0 1 0
## 330 13.638793 0 1 0
## 331 7.018369 0 1 0
## 332 4.863600 1 0 0
## 333 9.793280 0 1 0
## 334 5.515575 1 0 0
## 335 19.527871 0 1 0
## 336 14.268116 0 1 0
## 337 6.852828 0 1 0
## 338 41.685272 0 0 1
## 339 4.961085 1 0 0
## 340 20.504749 0 0 1
## 341 5.993640 1 0 0
## 342 6.017448 0 1 0
## 343 33.690571 0 0 1
## 344 3.556195 1 0 0
## 345 24.683993 0 0 1
## 346 5.809703 0 1 0
## 347 9.326966 0 1 0
## 348 18.334509 0 1 0
## 349 17.510738 0 1 0
## 350 7.027679 1 0 0
## 351 4.271347 1 0 0
## 352 4.339711 0 1 0
## 353 10.375783 0 1 0
## 354 8.336044 0 1 0
## 355 6.376198 1 0 0
## 356 18.979789 0 1 0
## 357 13.231006 0 1 0
## 358 29.627797 0 0 1
## 359 21.738015 0 0 1
## 360 5.036217 1 0 0
## 361 14.021403 0 1 0
## 362 5.049024 0 1 0
## 363 12.190385 0 1 0
## 364 32.124246 0 0 1
## 365 4.511465 1 0 0
## 366 4.829100 0 1 0
## 367 5.930220 1 0 0
## 368 8.522877 0 1 0
## 369 24.755001 0 0 1
## 370 7.021806 0 1 0
## 371 7.267124 0 1 0
## 372 7.876300 0 1 0
## 373 6.085785 0 1 0
## 374 6.094197 0 1 0
## 375 6.489400 0 1 0
## 376 14.393220 0 1 0
## 377 14.908918 0 1 0
## 378 5.842871 0 1 0
## 379 6.061984 0 1 0
## 380 10.079401 1 0 0
## 381 17.784859 0 1 0
## 382 22.843261 0 0 1
## 383 14.439115 0 1 0
## 384 7.455201 0 1 0
## 385 12.405091 0 1 0
## 386 7.517362 0 1 0
## 387 19.250615 0 1 0
## 388 21.072621 0 0 1
## 389 10.692359 0 1 0
## 390 8.010914 0 1 0
## 391 7.176618 0 1 0
## 392 11.658447 0 1 0
## 393 7.733971 0 1 0
## 394 8.753829 0 1 0
## 395 18.468848 0 1 0
## 396 28.762294 0 0 1
## 397 4.499861 0 1 0
## 398 5.747477 1 0 0
## 399 10.553091 0 1 0
## 400 13.425967 0 1 0
## 401 15.613148 0 1 0
## 402 13.586466 0 1 0
## 403 25.964138 0 0 1
## 404 17.988211 0 1 0
## 405 8.151942 0 1 0
## 406 13.927106 0 1 0
## 407 10.296812 0 1 0
## 408 8.983471 0 1 0
## 409 13.788731 0 1 0
## 410 13.044337 0 1 0
## 411 19.476535 0 1 0
## 412 15.236596 0 1 0
## 413 4.236654 1 0 0
## 414 5.771360 0 1 0
## 415 16.994470 0 1 0
## 416 10.738946 0 1 0
## 417 13.600778 0 1 0
## 418 18.953287 0 1 0
## 419 5.776360 0 1 0
## 420 11.391127 0 1 0
## 421 7.271955 0 1 0
## 422 59.095288 0 0 1
## 423 5.931022 0 1 0
## 424 6.950260 0 1 0
## 425 10.621064 0 1 0
## 426 20.284817 0 1 0
## 427 6.974925 0 1 0
## 428 17.874730 0 1 0
## 429 6.949101 0 1 0
## 430 6.708219 0 1 0
## 431 27.446210 0 0 1
## 432 5.326149 1 0 0
## 433 26.842458 0 0 1
## 434 10.173912 0 1 0
## 435 25.338889 0 0 1
## 436 16.196702 1 0 0
## 437 8.618344 0 1 0
## 438 5.966481 0 1 0
## 439 12.109653 0 1 0
## 440 12.705767 0 1 0
## 441 11.415358 0 1 0
## 442 6.376131 1 0 0
## 443 17.292413 0 1 0
## 444 7.550905 0 1 0
## 445 13.257041 0 1 0
## 446 11.944531 0 1 0
## 447 3.252344 1 0 0
## 448 8.862922 0 1 0
## 449 20.963515 0 1 0
## 450 5.624471 1 0 0
## 451 4.817141 1 0 0
## 452 4.526401 1 0 0
## 453 9.024521 0 1 0
## 454 4.899938 0 1 0
## 455 16.515319 0 1 0
## 456 7.255999 0 1 0
## 457 6.584965 0 1 0
## 458 11.248456 0 1 0
## 459 8.713702 0 1 0
## 460 14.793042 0 1 0
## 461 4.117487 1 0 0
## 462 8.382449 0 1 0
## 463 29.246672 0 0 1
## 464 32.948462 0 0 1
## 465 12.756377 0 1 0
## 466 20.327120 0 1 0
## 467 33.553081 0 0 1
## 468 4.490854 0 1 0
## 469 20.359406 0 1 0
## 470 4.697233 0 1 0
## 471 10.931053 0 1 0
## 472 5.207160 1 0 0
## 473 5.489221 0 1 0
## 474 4.965134 0 1 0
## 475 5.774458 0 1 0
## 476 9.479559 0 1 0
## 477 7.224633 0 1 0
## 478 8.589174 0 1 0
## 479 6.328741 0 1 0
## 480 19.071543 0 1 0
## 481 8.201049 0 1 0
## 482 7.381175 0 1 0
## 483 22.486572 0 0 1
## 484 4.364924 1 0 0
## 485 5.854865 0 1 0
## 486 4.775599 1 0 0
## 487 5.106241 0 1 0
## 488 4.567316 0 1 0
## 489 11.009977 0 1 0
## 490 3.694699 1 0 0
## 491 15.023951 0 1 0
## 492 5.012009 0 1 0
## 493 17.515464 0 1 0
## 494 21.117951 0 1 0
## 495 4.472277 1 0 0
## 496 7.737537 0 1 0
## 497 16.750113 0 1 0
## 498 21.771398 0 0 1
## 499 9.781351 0 1 0
## 500 3.800707 1 0 0
## 501 9.307143 0 1 0
## 502 6.890080 0 1 0
## 503 15.978518 0 1 0
## 504 7.504390 0 1 0
## 505 4.829615 1 0 0
## 506 7.038775 0 1 0
## 507 37.066953 0 0 1
## 508 12.319610 0 1 0
## 509 3.513176 1 0 0
## 510 11.876785 0 1 0
## 511 7.079029 0 1 0
## 512 6.387699 0 1 0
## 513 26.213427 0 0 1
## 514 5.361934 0 1 0
## 515 31.913832 0 0 1
## 516 14.233342 0 1 0
## 517 19.430030 0 1 0
## 518 10.317624 0 1 0
## 519 10.814183 0 1 0
## 520 5.069451 0 1 0
## 521 7.715734 0 1 0
## 522 8.479671 0 1 0
## 523 9.429924 0 1 0
## 524 25.370231 0 0 1
## 525 6.859487 1 0 0
## 526 7.578863 0 1 0
## 527 22.150188 0 0 1
## 528 5.035815 0 1 0
## 529 13.395882 1 0 0
## 530 9.962596 0 1 0
## 531 15.253514 0 1 0
## 532 7.443501 0 1 0
## 533 8.051168 0 1 0
## 534 5.873838 0 1 0
## 535 5.896406 1 0 0
## 536 13.778115 0 1 0
## 537 8.608128 0 1 0
## 538 6.120142 0 1 0
## 539 25.147926 0 0 1
## 540 26.189757 0 0 1
## 541 15.898513 0 1 0
## 542 5.780057 1 0 0
## 543 18.079177 0 1 0
## 544 23.142875 0 0 1
## 545 6.664551 0 1 0
## 546 6.055575 0 1 0
## 547 6.627707 0 1 0
## 548 16.257137 0 1 0
## 549 6.522958 0 1 0
## 550 6.393694 1 0 0
## 551 20.604526 0 0 1
## 552 10.556927 0 1 0
## 553 4.922080 0 1 0
## 554 16.408212 0 1 0
## 555 5.764832 0 1 0
## 556 26.316357 0 0 1
## 557 4.486100 1 0 0
## 558 7.841130 0 1 0
## 559 4.984246 0 1 0
## 560 16.189537 0 1 0
## 561 11.727260 0 1 0
## 562 8.450329 0 1 0
## 563 10.967744 0 1 0
## 564 4.938071 1 0 0
## 565 22.480252 0 0 1
## 566 8.505660 0 1 0
## 567 9.946052 0 1 0
## 568 6.191816 0 1 0
## 569 6.620178 0 1 0
## 570 11.048745 0 1 0
## 571 27.930217 0 0 1
## 572 3.828709 0 1 0
## 573 5.771468 0 1 0
## 574 6.010009 0 1 0
## 575 15.415297 0 1 0
## 576 4.212023 0 1 0
## 577 25.943566 0 0 1
## 578 10.025855 0 1 0
## 579 5.676316 0 1 0
## 580 5.296645 0 1 0
## 581 8.551752 0 1 0
## 582 9.461462 0 1 0
## 583 10.463816 0 1 0
## 584 9.143854 0 1 0
## 585 14.314692 0 1 0
## 586 8.911230 0 1 0
## 587 17.049552 0 1 0
## 588 4.717885 1 0 0
## 589 5.962961 0 1 0
## 590 10.194583 0 1 0
## 591 18.154868 0 1 0
## 592 13.357120 0 1 0
## 593 18.778324 0 1 0
## 594 6.629229 0 1 0
## 595 6.323671 0 1 0
## 596 7.808920 0 1 0
## 597 4.582308 1 0 0
## 598 7.667718 0 1 0
## 599 5.792826 0 1 0
## 600 13.543554 0 1 0
## 601 4.797786 1 0 0
## 602 11.017875 0 1 0
## 603 8.814970 0 1 0
## 604 3.420027 1 0 0
## 605 47.388002 0 0 1
## 606 4.960279 1 0 0
## 607 9.139187 0 1 0
## 608 11.374338 0 1 0
## 609 18.089348 0 1 0
## 610 9.520971 0 1 0
## 611 12.633942 0 1 0
## 612 5.417920 0 1 0
## 613 20.895347 0 1 0
## 614 5.128721 0 1 0
## 615 6.394471 0 1 0
## 616 11.236696 0 1 0
## 617 18.159109 0 1 0
## 618 20.800108 0 1 0
## 619 7.263215 0 1 0
## 620 8.829737 0 1 0
## 621 5.817026 0 1 0
## 622 19.299474 0 1 0
## 623 4.806268 1 0 0
## 624 4.224639 0 1 0
## 625 7.672345 0 1 0
## 626 20.820447 0 1 0
## 627 8.615963 0 1 0
## 628 9.128849 0 1 0
## 629 5.291964 1 0 0
## 630 10.242622 0 1 0
## 631 5.096522 1 0 0
## 632 4.253981 0 1 0
## 633 4.526066 0 1 0
## 634 23.198147 0 0 1
## 635 24.708972 0 0 1
## 636 10.709886 0 1 0
## 637 40.518137 0 0 1
## 638 8.858591 0 1 0
## 639 15.642781 0 1 0
## 640 7.184234 0 1 0
## 641 29.172975 0 0 1
## 642 4.013314 1 0 0
## 643 16.132384 0 1 0
## 644 25.779494 0 0 1
## 645 30.749372 0 0 1
## 646 10.686053 0 1 0
## 647 4.779053 0 1 0
## 648 6.295036 0 1 0
## 649 11.738951 0 1 0
## 650 4.537598 1 0 0
## 651 4.283603 1 0 0
## 652 8.559437 0 1 0
## 653 13.298497 0 1 0
## 654 16.711513 0 1 0
## 655 5.082748 1 0 0
## 656 15.814903 0 1 0
## 657 4.846588 0 1 0
## 658 12.688076 0 1 0
## 659 39.994603 0 0 1
## 660 5.033172 1 0 0
## 661 15.405091 0 1 0
## 662 12.240178 0 1 0
## 663 9.481463 0 1 0
## 664 29.154630 0 0 1
## 665 37.043217 0 0 1
## 666 18.907626 0 1 0
## 667 32.812987 0 0 1
## 668 4.297686 0 1 0
## 669 5.333535 0 1 0
## 670 6.446952 0 1 0
## 671 21.282880 0 1 0
## 672 24.879188 0 0 1
## 673 10.499755 0 1 0
## 674 40.312927 0 0 1
## 675 12.111310 0 1 0
## 676 20.750416 0 1 0
## 677 32.826782 0 0 1
## 678 7.244367 0 1 0
## 679 10.052928 0 1 0
## 680 5.138219 0 1 0
## 681 15.406538 0 1 0
## 682 7.431255 0 1 0
## 683 11.391992 0 1 0
## 684 9.086483 0 1 0
## 685 15.387208 0 1 0
## 686 8.882860 0 1 0
## 687 11.279864 0 1 0
## 688 37.211835 0 0 1
## 689 15.252173 0 1 0
## 690 6.195815 0 1 0
## 691 16.489027 0 1 0
## 692 9.416736 0 1 0
## 693 6.898985 0 1 0
## 694 4.986166 1 0 0
## 695 5.229555 0 1 0
## 696 30.070393 0 0 1
## 697 15.075517 0 1 0
## 698 5.219333 0 1 0
## 699 10.632291 0 1 0
## 700 14.083849 0 1 0
## 701 24.045118 0 0 1
## 702 4.293447 0 1 0
## 703 7.556785 0 1 0
## 704 11.599063 0 1 0
## 705 7.897812 0 1 0
## 706 4.353300 1 0 0
## 707 26.879802 0 0 1
## 708 18.629699 0 1 0
## 709 5.849203 0 1 0
## 710 8.387038 0 1 0
## 711 9.412294 0 1 0
## 712 16.079605 0 1 0
## 713 6.777882 0 1 0
## 714 9.608504 0 1 0
## 715 5.672070 0 1 0
## 716 19.446155 0 1 0
## 717 6.570576 0 1 0
## 718 16.000362 0 1 0
## 719 6.065461 0 1 0
## 720 4.756558 1 0 0
## 721 27.218248 0 0 1
## 722 7.164145 0 1 0
## 723 6.687505 0 1 0
## 724 40.448613 0 0 1
## 725 36.702127 0 0 1
## 726 14.537559 0 1 0
## 727 18.816304 0 1 0
## 728 14.058001 0 1 0
## 729 17.809420 0 1 0
## 730 3.528607 1 0 0
## 731 5.928080 1 0 0
## 732 9.270316 0 1 0
## 733 12.482074 0 1 0
## 734 34.623844 0 0 1
## 735 9.060879 0 1 0
## 736 4.864503 0 1 0
## 737 4.499854 1 0 0
## 738 7.254173 0 1 0
## 739 6.111321 1 0 0
## 740 8.228830 1 0 0
## 741 35.624814 0 0 1
## 742 13.550380 0 1 0
## 743 7.204323 0 1 0
## 744 18.549076 0 1 0
## 745 6.720920 0 1 0
## 746 2.899034 0 1 0
## 747 12.779281 0 1 0
## 748 15.798453 0 1 0
## 749 3.438641 0 1 0
## 750 8.552490 0 1 0
## 751 9.405298 0 1 0
## 752 42.808224 0 0 1
## 753 9.457132 0 1 0
## 754 40.718166 0 0 1
## 755 9.039712 0 1 0
## 756 15.599785 0 1 0
## 757 2.903734 1 0 0
## 758 8.444956 0 1 0
## 759 12.886667 0 1 0
## 760 9.512222 0 1 0
## 761 14.170067 0 1 0
## 762 10.834174 0 1 0
## 763 11.597444 0 1 0
## 764 7.754950 0 1 0
## 765 5.632811 0 1 0
## 766 5.342089 1 0 0
## 767 24.067087 0 0 1
## 768 17.362931 0 1 0
## 769 21.240921 0 0 1
## 770 7.718403 0 1 0
## 771 19.560701 0 1 0
## 772 10.705238 0 1 0
## 773 4.118247 0 1 0
## 774 3.791250 1 0 0
## 775 20.905579 0 0 1
## 776 13.338733 0 1 0
## 777 18.324902 0 1 0
## 778 19.043431 0 1 0
## 779 23.724273 0 0 1
## 780 7.678348 0 1 0
## 781 12.547588 0 1 0
## 782 24.129432 0 0 1
## 783 13.315128 0 1 0
## 784 9.370972 0 1 0
## 785 26.682683 0 0 1
## 786 7.453065 0 1 0
## 787 11.147000 0 1 0
## 788 5.459133 0 1 0
## 789 10.562732 0 1 0
## 790 34.324648 0 0 1
## 791 8.410987 0 1 0
## 792 15.166243 0 1 0
## 793 7.167321 0 1 0
## 794 4.015857 0 1 0
## 795 8.146647 0 1 0
## 796 4.953792 1 0 0
## 797 5.866762 1 0 0
## 798 11.708156 0 1 0
## 799 4.724716 1 0 0
## 800 5.630216 0 1 0
## 801 18.716826 0 1 0
## 802 8.367368 0 1 0
## 803 45.157103 0 0 1
## 804 14.896583 0 1 0
## 805 4.728871 1 0 0
## 806 34.797241 0 0 1
## 807 8.504734 0 1 0
## 808 7.876150 0 1 0
## 809 6.731238 0 1 0
## 810 3.954798 0 1 0
## 811 11.317785 0 1 0
## 812 14.834882 0 1 0
## 813 15.197306 0 1 0
## 814 14.646447 0 1 0
## 815 30.964842 0 0 1
## 816 15.987976 0 1 0
## 817 6.040280 1 0 0
## 818 35.108554 0 0 1
## 819 3.530954 1 0 0
## 820 10.049612 1 0 0
## 821 12.161719 0 1 0
## 822 20.248516 0 1 0
## 823 10.946744 0 1 0
## 824 4.854908 1 0 0
## 825 26.886995 0 0 1
## 826 6.553675 0 1 0
## 827 5.488807 0 1 0
## 828 14.030327 0 1 0
## 829 5.554347 0 1 0
## 830 56.485565 0 0 1
## 831 8.232838 0 1 0
## 832 13.193736 0 1 0
## 833 29.540419 0 0 1
## 834 11.553571 0 1 0
## 835 27.847339 0 0 1
## 836 17.578296 0 1 0
## 837 24.154366 0 0 1
## 838 8.125342 0 1 0
## 839 10.064358 0 1 0
## 840 7.202656 0 1 0
## 841 14.694169 0 1 0
## 842 12.165906 0 1 0
## 843 10.670699 0 1 0
## 844 5.022667 0 1 0
## 845 11.260466 0 1 0
## 846 4.941506 1 0 0
## 847 7.601123 0 1 0
## 848 18.047522 0 1 0
## 849 9.379980 0 1 0
## 850 15.036262 0 1 0
## 851 19.627941 0 1 0
## 852 4.932316 0 1 0
## 853 11.465209 1 0 0
## 854 4.065589 0 1 0
## 855 23.783541 0 0 1
## 856 18.974884 0 1 0
## 857 7.309180 1 0 0
## 858 6.778905 0 1 0
## 859 20.174708 0 1 0
## 860 14.086368 0 1 0
## 861 7.439381 1 0 0
## 862 8.588694 0 1 0
## 863 9.714359 0 1 0
## 864 6.719785 0 1 0
## 865 7.660087 0 1 0
## 866 8.816867 0 1 0
## 867 10.857410 0 1 0
## 868 9.161216 0 1 0
## 869 4.951661 1 0 0
## 870 4.804959 1 0 0
## 871 5.155758 0 1 0
## 872 17.040609 0 1 0
## 873 9.974358 0 1 0
## 874 4.238603 1 0 0
## 875 7.872056 0 1 0
## 876 9.608057 0 1 0
## 877 14.177128 0 1 0
## 878 9.157281 1 0 0
## 879 23.533264 0 0 1
## 880 14.517945 0 1 0
## 881 29.790258 0 0 1
## 882 5.373406 0 1 0
## 883 21.044071 0 0 1
## 884 9.140217 0 1 0
## 885 6.157556 1 0 0
## 886 10.793414 0 1 0
## 887 9.486539 0 1 0
## 888 43.232189 0 0 1
## 889 6.494899 0 1 0
## 890 10.013589 0 1 0
## 891 10.140635 0 1 0
## 892 6.188418 0 1 0
## 893 17.303164 0 1 0
## 894 24.210949 0 0 1
## 895 22.443464 0 0 1
## 896 4.877480 1 0 0
## 897 9.640884 0 1 0
## 898 10.782703 0 1 0
## 899 22.053574 0 0 1
## 900 11.286052 0 1 0
## 901 15.772342 0 1 0
## 902 14.667931 0 1 0
## 903 16.229528 0 1 0
## 904 11.646333 0 1 0
## 905 11.839377 0 1 0
## 906 5.576124 1 0 0
## 907 8.065882 0 1 0
## 908 11.545457 0 1 0
## 909 6.031073 0 1 0
## 910 4.509250 1 0 0
## 911 4.616738 1 0 0
## 912 9.624244 0 1 0
## 913 7.793338 0 1 0
## 914 10.362696 0 1 0
## 915 11.659745 0 1 0
## 916 9.628667 0 1 0
## 917 15.327269 0 1 0
## 918 34.917748 0 0 1
## 919 8.690655 0 1 0
## 920 19.984439 0 1 0
## 921 13.476867 0 1 0
## 922 9.193601 0 1 0
## 923 9.640399 0 1 0
## 924 53.529935 0 0 1
## 925 12.840598 0 1 0
## 926 17.624812 0 1 0
## 927 26.123801 0 0 1
## 928 9.016406 0 1 0
## 929 14.341198 0 1 0
## 930 4.640169 1 0 0
## 931 5.407539 0 1 0
## 932 7.585694 0 1 0
## 933 8.971526 0 1 0
## 934 3.505992 1 0 0
## 935 2.481118 1 0 0
## 936 10.423638 0 1 0
## 937 10.901006 0 1 0
## 938 34.480665 0 0 1
## 939 5.096764 1 0 0
## 940 4.165246 1 0 0
## 941 39.702716 0 0 1
## 942 21.359875 0 1 0
## 943 9.947408 1 0 0
## 944 13.834518 0 1 0
## 945 20.805949 0 1 0
## 946 5.662576 1 0 0
## 947 13.494068 0 1 0
## 948 38.913437 0 0 1
## 949 5.190204 0 1 0
## 950 17.890777 0 1 0
## 951 8.820599 0 1 0
## 952 9.697169 0 1 0
## 953 33.916645 0 0 1
## 954 12.736998 0 1 0
## 955 13.536001 0 1 0
## 956 13.162205 0 1 0
## 957 9.509858 0 1 0
## 958 5.411333 0 1 0
## 959 7.818260 0 1 0
## 960 37.110161 0 0 1
## 961 8.707991 0 1 0
## 962 19.066135 0 1 0
## 963 5.680717 1 0 0
## 964 17.059499 0 1 0
## 965 18.124216 0 1 0
## 966 12.758513 0 1 0
## 967 24.593772 0 0 1
## 968 21.226613 0 0 1
## 969 37.884313 0 0 1
## 970 4.707693 1 0 0
## 971 3.640579 1 0 0
## 972 9.917166 0 1 0
## 973 17.154990 0 1 0
## 974 13.635604 0 1 0
## 975 6.885037 0 1 0
## 976 18.821687 0 1 0
## 977 5.614852 1 0 0
## 978 5.085337 0 1 0
## 979 6.670012 0 1 0
## 980 12.534300 0 1 0
## 981 29.839078 0 0 1
## 982 5.426319 0 1 0
## 983 7.125389 0 1 0
## 984 12.655627 0 1 0
## 985 14.651827 0 1 0
## 986 5.164873 1 0 0
## 987 13.518185 0 1 0
## 988 4.685202 1 0 0
## 989 28.102796 0 0 1
## 990 8.416940 0 1 0
## 991 7.318212 0 1 0
## 992 18.468263 0 1 0
## 993 6.402860 1 0 0
## 994 15.861102 0 1 0
## 995 5.856849 1 0 0
## 996 4.902764 0 1 0
## 997 33.210001 0 0 1
## 998 9.997875 0 1 0
## 999 23.931790 0 0 1
## 1000 18.663170 0 1 0c1 = data x yang lebih kecil dari sama dengan 0 diberi nilai 1, sisanya diberi nilai 0 c2 = data x yang nilainya di antara 0 hingga 2 diberi nilai 1, sisanya diberi nilai 0 c3 = data x yang lebih besar dari 2 diberi nilai 1, sisanya diberi nilai 0
step.mod <- lm(y~c1+c2+c3)
plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70))
lines(data.x,lin.mod$fitted.values,col="red")
lines(data.x[ix], step.mod$fitted.values[ix],col="green")
Garis hijau merupakan fungsi tangga dengan 3 selang, x<=0 (c1), 0<x<=2 (c2), dan x>2 (c3), sedangkan garis merah merupakan fungsi regresi linier. Apabila garis hijau dan garis merah dibandingkan secara eksploratif, maka dapat dilihat bahwa garis fungsi tangga lebih mengikuti pola dibandingkan garis fungsi regresi linier.
summary(step.mod)##
## Call:
## lm(formula = y ~ c1 + c2 + c3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8218 -3.5774 -0.6672 2.7147 28.5759
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.5501 0.4087 74.75 <2e-16 ***
## c11 -25.1776 0.5819 -43.27 <2e-16 ***
## c21 -19.7065 0.4497 -43.82 <2e-16 ***
## c31 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.989 on 997 degrees of freedom
## Multiple R-squared: 0.699, Adjusted R-squared: 0.6984
## F-statistic: 1158 on 2 and 997 DF, p-value: < 2.2e-16Komparasi Model
nilai_AIC <- rbind(AIC(lin.mod),
AIC(pol.mod),
AIC(step.mod))
nama_model <- c("Linear","Poly (ordo=2)","Tangga (breaks=3)")
data.frame(nama_model,nilai_AIC)## nama_model nilai_AIC
## 1 Linear 5738.001
## 2 Poly (ordo=2) 2759.435
## 3 Tangga (breaks=3) 6057.180Komparasi model dilakukan berdasarkan nilai AIC-nya. Semakin kecil nilai AIC, maka semakin baik model. Dari ketiga nilai AIC di atas, model polynomial ordo 2 memiliki nilai terkecil.
MSE = function(pred,actual){
mean((pred-actual)^2)
}
nilai_MSE <- rbind(MSE(predict(lin.mod),y),
MSE(predict(pol.mod),y),
MSE(predict(step.mod),y))
nama_model <- c("Linear","Poly (ordo=2)","Tangga (breaks=3)")
data.frame(nama_model,nilai_MSE)## nama_model nilai_MSE
## 1 Linear 18.067659
## 2 Poly (ordo=2) 0.917189
## 3 Tangga (breaks=3) 24.811385Selain menggunakan AIC, komparasi model juga dapat dilakukan berdasarkan nilai MSE-nya.Semakin kecil nilai MSE, semakin bagus modelnya. Sesuai dengan hasil AIC di atas, model polynomial ordo 2 juga memiliki nilai MSE terkecil.
Maka dari itu, baik secara eksploratif maupun berdasarkan nilai AIC dan MSE-nya, model plynomial ordo 2 merupakan model terbaik bagi data tersebut.
Contoh Data TRICEPS
Data yang digunakan untuk ilustrasi berasal dari studi antropometri terhadap 892 perempuan di bawah 50 tahun di tiga desa Gambia di Afrika Barat. Data terdiri dari 3 Kolom yaitu Age, Intriceps dan tricepts. Berikut adalah penjelasannya pada masing-masing kolom:
- age : umur responden
- Intriceps : logaritma dari ketebalan lipatan kulit triceps
- triceps: ketebalan lipatan kulit triceps
Lipatan kulit trisep diperlukan untuk menghitung lingkar otot lengan atas. Ketebalannya memberikan informasi tentang cadangan lemak tubuh, sedangkan massa otot yang dihitung memberikan informasi tentang cadangan protein
data("triceps", package="MultiKink")Kita gambarkan dalam bentuk scatter plot
ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
theme_bw() 
Regresi Linier
mod_linear = lm(triceps~age,data=triceps)
summary(mod_linear)##
## Call:
## lm(formula = triceps ~ age, data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9512 -2.3965 -0.5154 1.5822 25.1233
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.19717 0.21244 29.17 <2e-16 ***
## age 0.21584 0.01014 21.28 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.007 on 890 degrees of freedom
## Multiple R-squared: 0.3372, Adjusted R-squared: 0.3365
## F-statistic: 452.8 on 1 and 890 DF, p-value: < 2.2e-16ringkasan_linear <- summary(mod_linear)
ringkasan_linear$r.squared## [1] 0.3372092AIC(mod_linear)## [1] 5011.515tabel_data <- rbind(data.frame("y_aktual"=triceps$triceps,
"y_pred"=predict(mod_linear),
"residuals"=residuals(mod_linear)))
tabel_data## y_aktual y_pred residuals
## 1 3.4 8.798074 -5.398073843
## 2 4.0 8.336171 -4.336171174
## 3 4.2 8.364231 -4.164230898
## 4 4.2 8.677202 -4.477202306
## 5 4.4 8.381498 -3.981497986
## 6 4.4 8.759222 -4.359222149
## 7 4.8 8.552014 -3.752013362
## 8 4.8 8.362072 -3.562072044
## 9 4.8 8.575756 -3.775756155
## 10 5.0 8.355597 -3.355597021
## 11 5.0 8.649143 -3.649142581
## 12 5.0 8.653460 -3.653459528
## 13 5.2 8.573598 -3.373598063
## 14 5.2 8.353439 -3.153438738
## 15 5.2 8.346963 -3.146963525
## 16 5.2 8.772173 -3.572173068
## 17 5.2 8.627558 -3.427558658
## 18 5.2 8.372864 -3.172864585
## 19 5.4 8.353439 -2.953438452
## 20 5.4 8.806708 -3.406707529
## 21 5.4 8.441934 -3.041933794
## 22 5.4 8.355597 -2.955596925
## 23 5.4 8.573598 -3.173597777
## 24 5.4 8.754906 -3.354905408
## 25 5.4 8.580073 -3.180072991
## 26 5.4 8.441934 -3.041933794
## 27 5.5 8.657776 -3.157776268
## 28 5.6 8.716054 -3.116053905
## 29 5.6 8.703103 -3.103103271
## 30 5.6 8.364231 -2.764230803
## 31 5.6 8.340488 -2.740488216
## 32 5.6 8.495894 -2.895894580
## 33 5.8 8.817500 -3.017499594
## 34 5.8 8.465677 -2.665676493
## 35 5.8 8.478627 -2.678626920
## 36 5.8 8.463518 -2.663518019
## 37 5.8 8.830450 -3.030450022
## 38 6.0 8.782965 -2.782964831
## 39 6.0 8.493736 -2.493736217
## 40 6.0 8.700945 -2.700944909
## 41 6.0 8.422508 -2.422508249
## 42 6.0 8.485103 -2.485102530
## 43 6.0 8.638351 -2.638350627
## 44 6.0 8.435459 -2.435458676
## 45 6.2 8.703103 -2.503103367
## 46 6.2 8.653460 -2.453459719
## 47 6.2 8.450568 -2.250567768
## 48 6.2 8.614608 -2.414608025
## 49 6.4 8.785123 -2.385123209
## 50 6.4 8.666410 -2.266409860
## 51 6.4 8.383657 -1.983656459
## 52 6.6 8.459201 -1.859201359
## 53 6.6 8.666410 -2.066410051
## 54 6.6 8.657776 -2.057776364
## 55 6.6 8.670727 -2.070726997
## 56 6.8 8.763539 -1.963539000
## 57 6.8 8.644826 -1.844825650
## 58 6.8 8.830450 -2.030450022
## 59 7.0 8.495894 -1.495894484
## 60 7.0 8.463518 -1.463518210
## 61 7.0 8.467835 -1.467835156
## 62 7.0 8.651301 -1.651301055
## 63 7.0 8.605974 -1.605974147
## 64 7.0 8.782965 -1.782964831
## 65 7.0 8.599499 -1.599498933
## 66 7.0 8.504528 -1.504528171
## 67 7.2 8.461360 -1.261359928
## 68 7.2 8.731163 -1.531162901
## 69 7.4 8.815341 -1.415341216
## 70 7.4 8.413875 -1.013874466
## 71 8.0 8.316745 -0.316745327
## 72 8.0 8.787282 -0.787281778
## 73 8.4 8.651301 -0.251301436
## 74 8.6 8.683678 -0.083677153
## 75 8.8 8.789440 0.010559940
## 76 8.8 8.336171 0.463829017
## 77 9.0 8.754906 0.245094497
## 78 9.0 8.828292 0.171708261
## 79 9.8 8.834767 0.965233032
## 80 9.8 8.694469 1.105530702
## 81 10.0 8.623242 1.376758479
## 82 10.2 8.642667 1.557332442
## 83 10.2 8.547697 1.652302997
## 84 13.0 8.474310 4.525689630
## 85 14.2 8.452726 5.747273759
## 86 16.0 8.631875 7.368124792
## 87 9.4 14.077578 -4.677578494
## 88 8.2 6.840384 1.359616282
## 89 21.2 13.043693 8.156307429
## 90 9.6 7.116662 2.483338564
## 91 10.2 11.457252 -1.257252372
## 92 13.2 10.406100 2.793900193
## 93 13.6 14.986275 -1.386274771
## 94 12.2 10.026217 2.173782828
## 95 6.8 10.065069 -3.265068484
## 96 6.4 10.097445 -3.697444854
## 97 14.0 9.352789 4.647211215
## 98 6.0 7.876427 -1.876426986
## 99 11.2 10.108237 1.091762494
## 100 7.4 7.382148 0.017852251
## 101 8.2 11.560856 -3.360856615
## 102 13.4 11.850085 1.549914374
## 103 17.6 11.191766 6.408234639
## 104 15.6 15.633802 -0.033801906
## 105 4.2 7.766347 -3.566347514
## 106 14.0 10.123346 3.876653784
## 107 8.2 7.013057 1.186942389
## 108 8.2 10.932755 -2.732755326
## 109 21.8 12.139314 9.660685173
## 110 7.8 6.682819 1.117181603
## 111 6.0 6.531729 -0.531728912
## 112 9.2 9.149897 0.050102770
## 113 6.4 7.133929 -0.733929096
## 114 17.4 14.915047 2.484952846
## 115 8.6 7.371356 1.228644594
## 116 5.8 7.248326 -1.448325404
## 117 8.0 8.959956 -0.959955722
## 118 10.0 12.860228 -2.860227641
## 119 13.6 14.120747 -0.520746372
## 120 7.6 8.182923 -0.582923172
## 121 8.2 6.421649 1.778350508
## 122 8.0 6.892186 1.107814299
## 123 9.0 7.403732 1.596267836
## 124 23.4 12.048660 11.351339370
## 125 6.4 7.054067 -0.654067389
## 126 9.0 11.297528 -2.297528459
## 127 5.2 8.867143 -3.667143624
## 128 5.8 8.131121 -2.331120765
## 129 8.2 10.414734 -2.214733700
## 130 4.2 17.151174 -12.951174136
## 131 4.2 7.677852 -3.477852172
## 132 22.6 10.552872 12.047127882
## 133 6.8 8.107378 -1.307378177
## 134 9.4 9.106728 0.293271219
## 135 8.0 13.037218 -5.037218326
## 136 11.0 10.395308 0.604692338
## 137 11.8 10.473011 1.326989552
## 138 8.4 6.501511 1.898488636
## 139 15.8 11.161548 4.638452249
## 140 17.4 11.068736 6.331263966
## 141 6.0 6.570581 -0.570580555
## 142 4.4 8.012408 -3.612407511
## 143 12.6 12.609850 -0.009849722
## 144 10.6 13.943756 -3.343755687
## 145 9.2 6.365530 2.834469525
## 146 7.8 6.365530 1.434469906
## 147 6.0 6.510145 -0.510144695
## 148 7.2 7.669218 -0.469218485
## 149 16.2 11.146439 5.053561722
## 150 5.2 7.496544 -2.296544541
## 151 6.6 7.129612 -0.529612443
## 152 6.0 8.267102 -2.267101679
## 153 6.2 6.678502 -0.478501935
## 154 13.6 10.453585 3.146415590
## 155 8.0 10.978082 -2.978081837
## 156 9.0 6.454026 2.545974322
## 157 11.2 11.105429 0.094570936
## 158 10.2 12.579632 -2.379632493
## 159 6.2 7.844051 -1.644050799
## 160 5.8 7.310920 -1.510919685
## 161 8.4 9.283719 -0.883719671
## 162 14.6 13.691221 0.908779500
## 163 6.6 7.295811 -0.695811071
## 164 9.6 9.303145 0.296855245
## 165 5.2 7.073493 -1.873493471
## 166 10.4 9.307462 1.092537741
## 167 7.4 9.244868 -1.844867500
## 168 10.2 6.777789 3.422210563
## 169 6.8 7.969239 -1.169238981
## 170 24.8 17.367016 7.432982913
## 171 19.0 14.122906 4.877094362
## 172 3.6 8.992332 -5.392332092
## 173 8.2 7.082127 1.117872842
## 174 7.0 7.338979 -0.338979410
## 175 6.2 9.020392 -2.820391721
## 176 4.6 7.798724 -3.198723796
## 177 9.4 6.833908 2.566091356
## 178 6.6 9.326888 -2.726887819
## 179 5.6 8.858510 -3.258509842
## 180 6.4 7.721020 -1.321020320
## 181 21.6 13.995559 7.604441780
## 182 15.2 17.367016 -2.167016514
## 183 12.4 12.109096 0.290903767
## 184 5.8 7.772823 -1.972822449
## 185 13.6 13.693379 -0.093378561
## 186 19.4 14.241619 5.158380837
## 187 7.0 9.020392 -2.020391530
## 188 6.8 9.119679 -2.319678842
## 189 7.0 9.132630 -2.132629666
## 190 10.4 9.430492 0.969507652
## 191 7.0 6.404382 0.595618086
## 192 14.2 13.330764 0.869236128
## 193 8.8 9.093778 -0.293777781
## 194 9.2 6.941829 2.258170357
## 195 6.4 9.074352 -2.674352029
## 196 6.8 8.148388 -1.348388138
## 197 5.0 8.938371 -3.938371402
## 198 7.2 8.269260 -1.069260342
## 199 7.6 9.708929 -2.108928928
## 200 8.2 9.197382 -0.997382405
## 201 8.6 7.302286 1.297714193
## 202 5.4 6.954780 -1.554779887
## 203 6.8 9.255660 -2.455659565
## 204 5.0 8.247676 -3.247675832
## 205 8.4 6.546838 1.853161729
## 206 24.8 12.773890 12.026108876
## 207 5.4 7.686486 -2.286485573
## 208 6.8 6.831750 -0.031749650
## 209 8.0 7.246167 0.753832776
## 210 6.2 9.272927 -3.072927320
## 211 6.4 6.963414 -0.563413574
## 212 7.0 6.516620 0.483380040
## 213 6.6 9.698136 -3.098136562
## 214 7.2 6.590006 0.609993433
## 215 9.0 9.756414 -0.756414008
## 216 7.0 9.596691 -2.596690697
## 217 7.8 9.683028 -1.883027376
## 218 11.2 13.546606 -2.346606250
## 219 12.2 6.980681 5.219318715
## 220 16.0 9.706770 6.293229640
## 221 11.0 15.849645 -4.849644666
## 222 13.4 14.409976 -1.009975955
## 223 6.0 9.639859 -3.639859132
## 224 8.2 12.270978 -4.070977826
## 225 18.2 14.772591 3.427409928
## 226 5.8 8.118170 -2.318170131
## 227 7.6 7.453376 0.146623989
## 228 6.6 6.339629 0.260370693
## 229 7.0 7.347613 -0.347613097
## 230 9.8 9.836276 -0.036275678
## 231 6.0 7.733971 -1.733970946
## 232 17.2 13.848786 3.351215044
## 233 9.8 6.715195 3.084805226
## 234 11.0 9.434809 1.565191087
## 235 6.2 9.432650 -3.232650631
## 236 12.0 12.914188 -0.914188236
## 237 10.4 6.766997 3.633002481
## 238 9.2 15.849645 -6.649644857
## 239 7.6 6.421649 1.178350604
## 240 15.4 14.483362 0.916637604
## 241 6.0 8.157022 -2.157022016
## 242 8.4 6.497194 1.902805480
## 243 7.4 15.849645 -8.449644571
## 244 19.0 10.278753 8.721247420
## 245 7.4 9.654968 -2.254968143
## 246 9.6 13.652369 -4.052368806
## 247 16.2 11.539272 4.660728659
## 248 12.2 10.600358 1.599642134
## 249 10.8 11.776699 -0.976698612
## 250 4.8 8.085794 -3.285793857
## 251 14.0 10.090970 3.909030058
## 252 9.2 9.542730 -0.342730293
## 253 5.4 7.936863 -2.536862802
## 254 5.6 7.278544 -1.678543697
## 255 7.4 9.572948 -2.172948014
## 256 5.0 9.156372 -4.156372253
## 257 10.0 9.566473 0.433527310
## 258 7.2 8.970748 -1.770748073
## 259 7.8 6.466976 1.333023982
## 260 9.2 11.373073 -2.173073564
## 261 7.2 14.552432 -7.352431701
## 262 5.4 7.995140 -2.595140137
## 263 7.8 7.060543 0.739457441
## 264 6.0 7.528921 -1.528920728
## 265 12.0 9.182273 2.817726686
## 266 8.6 6.309411 2.290589113
## 267 10.4 10.026217 0.373782637
## 268 12.6 10.729863 1.870137197
## 269 8.0 10.276594 -2.276594107
## 270 8.6 7.166305 1.433694916
## 271 5.0 13.853103 -8.853102665
## 272 8.2 9.037659 -0.837659095
## 273 9.2 7.205157 1.994842649
## 274 3.2 9.777998 -6.577998280
## 275 16.8 12.074561 4.725438133
## 276 6.0 7.867793 -1.867793196
## 277 12.2 12.473870 -0.273869777
## 278 9.0 14.768274 -5.768273889
## 279 16.2 14.125064 2.074937063
## 280 9.4 10.738497 -1.338497459
## 281 7.0 8.893044 -1.893044494
## 282 9.0 14.770433 -5.770432774
## 283 11.0 14.768274 -3.768273889
## 284 15.6 14.768274 0.831726493
## 285 9.0 9.691661 -0.691661459
## 286 15.4 10.280911 5.119088565
## 287 11.0 12.184641 -1.184640766
## 288 12.0 13.907063 -1.907063260
## 289 8.2 6.387115 1.812885282
## 290 6.4 9.223283 -2.823283386
## 291 29.2 13.963182 15.236818847
## 292 7.2 8.888728 -1.688727944
## 293 24.8 12.398325 12.401674566
## 294 10.6 9.711087 0.888913075
## 295 13.0 13.395516 -0.395516230
## 296 7.8 7.114503 0.685496846
## 297 13.0 13.786191 -0.786191231
## 298 10.0 14.392709 -4.392708611
## 299 7.0 7.300128 -0.300127819
## 300 7.8 15.202118 -7.402117340
## 301 5.2 7.207316 -2.007315721
## 302 14.0 14.321480 -0.321480231
## 303 15.0 15.417960 -0.417959909
## 304 4.0 9.175798 -5.175798100
## 305 11.8 7.153355 4.646645255
## 306 8.0 9.639859 -1.639859132
## 307 10.2 6.797215 3.402784768
## 308 14.2 9.182273 5.017726495
## 309 12.0 15.847486 -3.847485781
## 310 11.4 8.912470 2.487529278
## 311 12.0 11.081686 0.918313920
## 312 6.6 6.943988 -0.343987969
## 313 8.2 7.205157 0.994842649
## 314 20.2 15.199959 5.000042117
## 315 14.6 16.065487 -1.465486663
## 316 14.2 16.274854 -2.074853783
## 317 10.2 6.413016 3.786984195
## 318 10.0 16.495013 -6.495012917
## 319 20.2 15.851803 4.348198035
## 320 11.0 15.849645 -4.849644666
## 321 8.4 6.855492 1.544507139
## 322 14.0 12.180324 1.819676180
## 323 10.2 6.859809 3.340190486
## 324 9.8 6.803690 2.996309884
## 325 5.2 7.518129 -2.318128758
## 326 19.0 14.701363 4.298636722
## 327 6.2 9.171481 -2.971481345
## 328 8.2 6.883552 1.316447796
## 329 6.6 12.031393 -5.431392970
## 330 7.2 6.341788 0.858212176
## 331 6.4 9.708929 -3.308928738
## 332 15.2 10.529130 4.670870103
## 333 6.6 7.151197 -0.551196661
## 334 13.2 9.965781 3.234218842
## 335 27.0 16.935332 10.064668433
## 336 9.0 16.931015 -7.931014620
## 337 13.8 16.926698 -3.126697483
## 338 19.0 15.204276 3.795724408
## 339 9.4 14.511422 -5.111422136
## 340 9.8 6.572739 3.227261214
## 341 30.4 13.958865 16.441134649
## 342 8.4 7.129612 1.270387271
## 343 8.4 6.486402 1.913597588
## 344 14.2 9.331205 4.868795139
## 345 6.4 8.191557 -1.791556668
## 346 9.2 6.438917 2.761083109
## 347 6.2 7.218108 -1.018107881
## 348 9.0 7.390782 1.609218366
## 349 7.0 12.402642 -5.402641618
## 350 8.2 7.712387 0.487613081
## 351 6.0 7.824625 -1.824624761
## 352 8.2 7.077810 1.122189686
## 353 8.2 6.898661 1.301338843
## 354 9.8 15.847486 -6.047485591
## 355 10.0 15.847486 -5.847485781
## 356 6.2 7.518129 -1.318128758
## 357 7.6 9.346314 -1.746313666
## 358 5.8 9.799583 -3.999582458
## 359 11.4 12.180324 -0.780324201
## 360 6.8 8.923263 -2.123262310
## 361 17.0 15.417960 1.582040091
## 362 20.6 16.497172 4.102828580
## 363 12.4 10.900379 1.499620757
## 364 7.6 8.286528 -0.686527621
## 365 7.0 9.393799 -2.393798952
## 366 11.4 10.166514 1.233485174
## 367 9.2 6.479927 2.720073070
## 368 6.4 9.136946 -2.736946311
## 369 19.8 11.204717 8.595282655
## 370 8.8 6.788581 2.011418836
## 371 10.2 13.697696 -3.497696080
## 372 27.4 13.475379 13.924621116
## 373 5.0 11.955848 -6.955847960
## 374 21.0 16.065487 4.934512955
## 375 11.8 17.142540 -5.342539862
## 376 5.6 9.534096 -3.934096511
## 377 4.6 13.978291 -9.378290912
## 378 12.4 11.964481 0.435518177
## 379 19.2 13.598409 5.601592170
## 380 8.0 12.173849 -4.173848812
## 381 7.6 6.760522 0.839478084
## 382 8.8 6.823116 1.976884037
## 383 8.4 6.883552 1.516447605
## 384 8.2 6.270560 1.929440164
## 385 6.2 7.187890 -0.987889977
## 386 7.8 7.321712 0.478288155
## 387 7.6 14.552432 -6.952431606
## 388 20.2 15.208593 4.991408224
## 389 21.5 13.296229 8.203771067
## 390 18.2 12.098304 6.101696866
## 391 17.2 12.070244 5.129756606
## 392 6.0 7.503020 -1.503019667
## 393 8.4 6.462659 1.937340253
## 394 11.8 9.911820 1.888179818
## 395 6.0 6.633175 -0.633174836
## 396 4.2 9.510354 -5.310354019
## 397 7.4 9.415383 -2.015382971
## 398 6.2 7.295811 -1.095811166
## 399 7.6 6.691452 0.908547630
## 400 11.2 11.355806 -0.155806191
## 401 7.0 8.157022 -1.157022016
## 402 6.8 6.831750 -0.031749650
## 403 18.2 15.417960 2.782040854
## 404 6.4 6.704403 -0.304402709
## 405 8.0 7.138246 0.861753965
## 406 9.0 6.413016 2.586984386
## 407 8.8 14.856769 -6.056769040
## 408 9.2 10.460060 -1.260060402
## 409 8.0 6.721670 1.278329770
## 410 5.8 7.807357 -2.007357197
## 411 20.4 11.770223 8.629776236
## 412 15.8 12.262344 3.537656037
## 413 18.0 9.713246 8.286754221
## 414 19.0 11.563015 7.436985103
## 415 9.2 11.595391 -2.395391362
## 416 16.0 11.647193 4.352806707
## 417 4.6 8.264943 -3.664943301
## 418 7.4 9.244868 -1.844867500
## 419 6.4 9.838434 -3.438434247
## 420 23.4 9.335521 14.064478208
## 421 6.2 7.556980 -1.356980452
## 422 17.8 13.261694 4.538305051
## 423 13.2 13.261694 -0.061694376
## 424 6.4 7.779298 -1.379297758
## 425 6.2 7.921754 -1.721753981
## 426 7.8 6.372006 1.427994628
## 427 8.8 9.283719 -0.483719099
## 428 8.0 7.015216 0.984784158
## 429 8.0 6.417332 1.582667542
## 430 11.2 6.725987 4.474012736
## 431 6.6 8.163497 -1.563497325
## 432 7.6 7.090761 0.509239251
## 433 5.0 9.262135 -4.262134969
## 434 14.8 9.631225 5.168774746
## 435 4.0 8.847718 -4.847717586
## 436 6.0 7.921754 -1.921753791
## 437 6.2 6.902978 -0.702978000
## 438 13.8 13.661003 0.138997111
## 439 36.8 14.125064 22.674935537
## 440 10.0 10.993191 -0.993190738
## 441 14.0 13.477537 0.522463436
## 442 6.0 7.608782 -1.608782383
## 443 7.6 13.557398 -5.957398109
## 444 5.8 7.153355 -1.353354745
## 445 10.8 6.555472 4.244528614
## 446 26.2 13.477537 12.722464199
## 447 7.6 6.946146 0.653853609
## 448 7.4 7.418841 -0.018840970
## 449 5.8 7.684327 -1.884327004
## 450 9.2 13.473220 -4.273219808
## 451 5.6 7.932546 -2.332546046
## 452 10.4 9.713246 0.686753839
## 453 20.2 13.261694 6.938306577
## 454 11.4 10.665111 0.734888983
## 455 9.0 12.219176 -3.219175514
## 456 7.0 8.146230 -1.146229856
## 457 7.8 7.092919 0.707081064
## 458 10.4 8.297319 2.102680139
## 459 5.2 7.675694 -2.475693699
## 460 8.0 6.741096 1.258903975
## 461 7.0 9.430492 -2.430491967
## 462 4.6 8.023200 -3.423199861
## 463 14.8 13.693379 1.106621248
## 464 9.8 6.626700 3.173300620
## 465 13.2 12.232126 0.967873868
## 466 6.0 6.477768 -0.477768317
## 467 13.4 12.542939 0.857060537
## 468 21.4 13.475379 7.924621116
## 469 8.4 11.431351 -3.031351296
## 470 16.4 12.396166 4.003833420
## 471 17.2 14.086212 3.113788757
## 472 7.2 6.857651 0.342348908
## 473 12.0 9.883761 2.116238956
## 474 27.2 16.719489 10.480511575
## 475 5.4 8.893044 -3.493044399
## 476 5.8 6.926720 -1.126720309
## 477 6.2 7.997299 -1.797298896
## 478 8.8 6.717353 2.082646804
## 479 7.0 7.572089 -0.572089162
## 480 5.0 7.779298 -2.779297854
## 481 14.8 9.549205 5.250794874
## 482 6.0 9.115362 -3.115362292
## 483 12.8 9.976573 2.823426858
## 484 8.2 7.000107 1.199892919
## 485 16.0 12.314146 3.685853724
## 486 8.2 6.972047 1.227952453
## 487 8.8 6.993632 1.806368566
## 488 7.2 6.393590 0.806410003
## 489 26.0 14.768274 11.231726111
## 490 7.0 9.035500 -2.035500431
## 491 7.4 7.233217 0.166783402
## 492 11.2 13.641577 -2.441577424
## 493 7.6 7.060543 0.539457155
## 494 11.4 14.804968 -3.404967903
## 495 8.4 16.928857 -8.528856940
## 496 7.0 7.440425 -0.440425385
## 497 6.2 11.398975 -5.198974831
## 498 7.8 8.100903 -0.300902758
## 499 5.8 7.852684 -2.052684105
## 500 14.0 10.902537 3.097463077
## 501 5.8 8.133279 -2.333279238
## 502 9.6 6.499353 3.100647821
## 503 23.2 12.044343 11.155657460
## 504 8.2 8.299478 -0.099478144
## 505 9.6 6.648284 2.951716593
## 506 6.6 15.847486 -9.247485877
## 507 10.6 10.801091 -0.201090771
## 508 9.0 6.758363 2.241636601
## 509 5.4 7.334663 -1.934662471
## 510 5.6 7.513812 -1.913811819
## 511 15.0 11.966640 3.033360086
## 512 10.2 6.454026 3.745974131
## 513 9.0 6.400065 2.599934929
## 514 12.7 12.609850 0.090149706
## 515 13.2 10.039167 3.160832401
## 516 8.6 6.931037 1.668963038
## 517 6.4 7.967081 -1.567080603
## 518 7.4 6.410857 0.989142903
## 519 8.2 11.893253 -3.693253664
## 520 6.6 9.581582 -2.981581892
## 521 4.2 8.154864 -3.954863734
## 522 9.8 15.631643 -5.831643212
## 523 17.5 13.924330 3.575669778
## 524 8.4 10.293861 -1.893861862
## 525 8.2 6.922404 1.277596204
## 526 7.0 6.525254 0.474746353
## 527 8.2 9.044134 -0.844134308
## 528 11.4 13.477537 -2.077536946
## 529 6.8 7.272068 -0.472068094
## 530 13.4 15.849645 -2.449645048
## 531 9.2 7.341138 1.858861926
## 532 14.4 12.178165 2.221834272
## 533 11.4 16.067645 -4.667645488
## 534 11.2 11.496103 -0.296103654
## 535 27.2 12.594741 14.605259560
## 536 8.4 6.833908 1.566091356
## 537 9.0 6.898661 2.101339034
## 538 12.6 15.199959 -2.599958264
## 539 10.6 6.479927 4.120073642
## 540 8.6 14.984116 -6.384115886
## 541 7.2 7.397257 -0.197257039
## 542 7.2 7.395098 -0.195098668
## 543 5.3 13.259536 -7.959535933
## 544 6.6 8.910312 -2.310311963
## 545 16.0 7.032483 8.967516784
## 546 8.0 10.062910 -2.062910202
## 547 7.8 7.423158 0.376842282
## 548 8.8 11.720580 -2.920579544
## 549 18.2 16.065487 2.134513718
## 550 10.0 12.583949 -2.583949249
## 551 9.2 12.374582 -3.174582481
## 552 11.4 6.516620 4.883379659
## 553 8.0 9.398116 -1.398115693
## 554 16.0 14.986275 1.013724848
## 555 15.6 10.065069 5.534931706
## 556 7.0 7.334663 -0.334662566
## 557 14.8 11.003983 3.796017087
## 558 4.0 11.601867 -7.601866591
## 559 24.0 15.202118 8.797882469
## 560 10.0 10.496753 -0.496753432
## 561 8.8 6.890027 1.909972912
## 562 6.8 8.979382 -2.179381378
## 563 7.0 7.136088 -0.136087562
## 564 5.4 9.145580 -3.745579998
## 565 5.4 7.287177 -1.887177193
## 566 23.6 10.688853 12.911147364
## 567 6.0 7.882902 -1.882902200
## 568 7.0 8.215300 -1.215299557
## 569 36.6 11.476678 25.123320858
## 570 9.8 7.198682 2.601318348
## 571 27.2 15.417960 11.782040854
## 572 11.4 6.786423 4.613576685
## 573 6.8 6.404382 0.395618276
## 574 10.0 6.989315 3.010685219
## 575 15.8 10.116871 5.683129394
## 576 9.0 10.179465 -1.179465284
## 577 5.4 8.003774 -2.603773824
## 578 33.4 14.986275 18.413726374
## 579 8.6 6.443234 2.156766837
## 580 9.0 7.071335 1.928665090
## 581 11.0 14.986275 -3.986275152
## 582 12.6 11.155073 1.444927860
## 583 8.0 6.853334 1.146665942
## 584 16.4 10.291703 6.108296611
## 585 10.2 6.419491 3.780508930
## 586 7.4 6.417332 0.982667638
## 587 6.8 8.087953 -1.287952330
## 588 8.0 6.432441 1.567558565
## 589 7.8 6.384956 1.415044085
## 590 9.2 15.415801 -6.215801215
## 591 6.8 6.313728 0.486272066
## 592 24.0 14.094846 9.905154102
## 593 15.0 14.584808 0.415191804
## 594 6.0 6.570581 -0.570580555
## 595 7.8 6.512303 1.287697074
## 596 7.0 9.970098 -2.970097913
## 597 8.4 9.037659 -0.637659285
## 598 8.0 7.710228 0.289771642
## 599 6.0 6.441075 -0.441075122
## 600 6.4 12.247235 -5.847235158
## 601 7.6 7.552663 0.047336487
## 602 11.8 10.062910 1.737089989
## 603 5.6 7.589357 -1.989356631
## 604 5.2 7.630367 -2.430366791
## 605 7.0 7.157672 -0.157671779
## 606 7.7 7.004424 0.695576076
## 607 8.8 6.445392 2.354608225
## 608 8.8 11.964481 -3.164481251
## 609 11.0 6.372006 4.627994437
## 610 11.4 14.213559 -2.813559423
## 611 6.0 8.046942 -2.046942354
## 612 5.6 6.982840 -1.382839611
## 613 8.4 6.333154 2.066845685
## 614 7.4 7.639000 -0.239000192
## 615 11.2 9.903187 1.296812918
## 616 7.0 10.375882 -3.375881815
## 617 9.0 6.415174 2.584825964
## 618 7.0 6.413016 0.586984386
## 619 8.6 6.410857 2.189143189
## 620 12.0 16.067645 -4.067645106
## 621 8.2 9.417542 -1.217541730
## 622 15.6 16.713014 -1.113013799
## 623 6.2 7.343296 -1.143296444
## 624 5.2 8.210983 -3.010982801
## 625 18.6 9.983049 8.616951629
## 626 9.2 12.024917 -2.824917646
## 627 4.4 7.941180 -3.541179542
## 628 4.8 7.615258 -2.815257509
## 629 8.4 9.262135 -0.862135351
## 630 10.0 10.714754 -0.714754284
## 631 6.2 7.533238 -1.333237762
## 632 5.4 6.754047 -1.354046460
## 633 6.8 6.374164 0.425836206
## 634 8.2 7.274227 0.925773154
## 635 9.2 8.994490 0.205509340
## 636 10.0 13.015634 -3.015634006
## 637 5.2 9.139105 -3.939105070
## 638 10.2 11.664461 -1.464460858
## 639 9.0 9.616117 -0.616116544
## 640 16.2 12.603375 3.596625667
## 641 9.8 6.721670 3.078329961
## 642 6.6 12.055135 -5.455135352
## 643 5.8 6.808007 -1.008006960
## 644 6.4 9.141263 -2.741263257
## 645 6.2 7.520287 -1.320287231
## 646 7.0 6.766997 0.233002862
## 647 11.2 6.529570 4.670429319
## 648 7.2 7.585040 -0.385039883
## 649 8.2 6.775631 1.424368985
## 650 9.0 6.782106 2.217893910
## 651 4.4 7.345455 -2.945454631
## 652 4.6 7.764189 -3.164189048
## 653 8.2 6.792898 1.407101611
## 654 7.2 7.794407 -0.594407048
## 655 14.2 11.273786 2.926213732
## 656 8.4 6.525254 1.874745972
## 657 8.0 6.553313 1.446686845
## 658 22.2 15.204276 6.995725171
## 659 7.6 6.546838 1.053162015
## 660 5.4 7.187890 -1.787889691
## 661 6.4 6.566264 -0.166263616
## 662 17.6 11.405450 6.194550733
## 663 7.4 8.139755 -0.739754547
## 664 6.8 7.507337 -0.707336320
## 665 5.2 16.065487 -10.865487235
## 666 6.2 7.960605 -1.760605675
## 667 5.2 8.232567 -3.032567122
## 668 9.8 6.529570 3.270429700
## 669 14.0 12.719930 1.280070234
## 670 8.0 7.326029 0.673971121
## 671 16.8 9.829800 6.970198787
## 672 9.2 9.311779 -0.111779014
## 673 8.0 9.898870 -1.898869945
## 674 17.6 10.447110 7.152890598
## 675 13.6 15.130889 -1.530888769
## 676 6.8 6.732462 0.067537852
## 677 9.0 10.021900 -1.021900035
## 678 11.0 11.970957 -0.970956861
## 679 9.0 9.413225 -0.413224593
## 680 14.0 11.299687 2.700313068
## 681 6.0 8.234725 -2.234725198
## 682 9.6 7.185731 2.414269069
## 683 9.2 10.887428 -1.687428213
## 684 7.2 9.719721 -2.519720978
## 685 11.6 6.905136 4.694864150
## 686 6.4 9.115362 -2.715362197
## 687 5.6 8.929738 -3.329737810
## 688 4.0 7.537554 -3.537554414
## 689 5.8 7.228900 -1.428899660
## 690 20.4 13.332923 7.067077053
## 691 19.8 14.125064 5.674935537
## 692 14.0 14.554590 -0.554590395
## 693 14.6 11.103270 3.496729981
## 694 10.2 9.199541 1.000459121
## 695 12.0 9.350630 2.649369688
## 696 7.0 6.326679 0.673321332
## 697 6.2 7.546188 -1.346188292
## 698 13.6 13.725756 -0.125755247
## 699 18.0 11.206875 6.793125357
## 700 11.0 10.701804 0.298196144
## 701 6.2 7.658426 -1.458426325
## 702 5.0 7.865635 -2.865634826
## 703 10.4 7.187890 3.212109833
## 704 9.2 9.754256 -0.554255725
## 705 7.0 7.526762 -0.526762254
## 706 5.6 6.348263 -0.748262993
## 707 15.2 17.360541 -2.160541507
## 708 14.4 12.186799 2.213200379
## 709 6.0 8.940530 -2.940529875
## 710 18.0 10.280911 7.719088946
## 711 17.0 13.475379 3.524621497
## 712 7.5 7.382148 0.117852156
## 713 12.8 15.847486 -3.047485591
## 714 7.0 6.810166 0.189834428
## 715 5.4 9.093778 -3.693777876
## 716 4.6 8.159180 -3.559180585
## 717 9.2 7.064860 2.135140216
## 718 12.0 10.436317 1.563682582
## 719 19.2 11.813392 7.386608740
## 720 11.0 9.184432 1.815568213
## 721 10.2 8.174289 2.025710419
## 722 6.0 7.444742 -1.444742229
## 723 9.0 11.401133 -2.401133113
## 724 7.8 7.477119 0.322881688
## 725 13.8 13.693379 0.106621248
## 726 8.2 15.417960 -7.217960100
## 727 8.6 7.496544 1.103456031
## 728 11.2 13.907063 -2.707063450
## 729 8.0 7.142563 0.857437122
## 730 13.0 8.888728 4.111272247
## 731 18.6 14.455302 4.144698106
## 732 8.2 7.392940 0.807059805
## 733 17.0 9.935563 7.064436834
## 734 8.2 8.299478 -0.099478144
## 735 6.4 7.768506 -1.368505701
## 736 12.4 14.915047 -2.515047154
## 737 6.6 6.676343 -0.076343418
## 738 15.0 14.986275 0.013724848
## 739 9.0 10.507545 -1.507545386
## 740 5.4 7.923912 -2.523912168
## 741 7.0 16.281329 -9.281329423
## 742 6.8 8.103061 -1.303061231
## 743 14.5 14.543798 -0.043797618
## 744 8.8 14.472570 -5.672569870
## 745 7.0 7.522445 -0.522445411
## 746 6.2 6.611591 -0.411590784
## 747 10.6 8.871460 1.728540002
## 748 6.8 6.527412 0.272588122
## 749 6.8 8.992332 -2.192331806
## 750 16.0 16.065487 -0.065487045
## 751 7.2 7.332504 -0.132504387
## 752 7.6 15.204276 -7.604275688
## 753 11.6 9.816850 1.783150359
## 754 11.0 13.924330 -2.924330222
## 755 8.0 7.202999 0.797001313
## 756 5.4 6.333154 -0.933153838
## 757 11.6 12.825692 -1.225692101
## 758 11.2 13.689062 -2.489062187
## 759 7.4 7.466326 -0.066326351
## 760 6.4 7.669218 -1.269218199
## 761 7.0 7.067018 -0.067018066
## 762 4.2 14.770433 -10.570432964
## 763 8.8 7.330346 1.469654468
## 764 6.2 6.669868 -0.469868248
## 765 6.6 7.423158 -0.823158004
## 766 8.4 6.590006 1.809993243
## 767 4.0 11.314796 -7.314795833
## 768 7.0 8.113854 -1.113853582
## 769 8.0 11.066577 -3.066577180
## 770 20.0 12.063769 7.936230850
## 771 19.0 15.202118 3.797882469
## 772 5.4 10.673744 -5.273744021
## 773 24.6 16.706538 7.893462032
## 774 5.0 8.005932 -3.005932392
## 775 11.0 9.801741 1.198258878
## 776 5.0 8.882252 -3.882252334
## 777 8.6 7.224583 1.375417375
## 778 12.3 16.713014 -4.413013990
## 779 6.0 7.835417 -1.835416922
## 780 10.6 9.864335 0.735665184
## 781 7.0 6.641809 0.358191477
## 782 13.6 14.768274 -1.168273507
## 783 7.6 12.137156 -4.537155686
## 784 13.6 14.697046 -1.097045951
## 785 7.8 6.993632 0.806368566
## 786 6.8 7.967081 -1.167080508
## 787 5.4 7.155513 -1.755513313
## 788 10.0 9.937722 0.062278361
## 789 5.6 7.524604 -1.924603979
## 790 5.6 9.348472 -3.748472140
## 791 28.0 11.621292 16.378707973
## 792 10.0 7.518129 2.481871433
## 793 9.4 6.922404 2.477596014
## 794 5.2 8.157022 -2.957022207
## 795 14.0 10.496753 3.503246568
## 796 9.8 14.554590 -4.754590204
## 797 6.4 6.507986 -0.107986178
## 798 7.4 7.582881 -0.182881227
## 799 6.0 7.561297 -1.561297105
## 800 5.0 17.144699 -12.144698937
## 801 8.4 8.295161 0.104838406
## 802 17.4 11.925630 5.474369460
## 803 9.4 9.497403 -0.097403576
## 804 9.2 12.035709 -2.835709600
## 805 17.8 11.491787 6.308212308
## 806 7.0 9.419700 -2.419700013
## 807 11.0 11.707629 -0.707629307
## 808 6.4 6.434600 -0.034599762
## 809 11.0 10.401783 0.598217330
## 810 7.0 7.701595 -0.701594569
## 811 10.2 9.855702 0.344298093
## 812 9.2 9.424017 -0.224016944
## 813 6.2 7.362722 -1.162722291
## 814 8.0 7.308762 0.691238494
## 815 5.8 7.483594 -1.683593629
## 816 24.0 10.239901 13.760099114
## 817 8.0 7.030325 0.969675206
## 818 7.8 6.486402 1.313598161
## 819 20.8 13.259536 7.540463113
## 820 8.6 11.010458 -2.410457730
## 821 11.6 10.578773 1.021227027
## 822 6.6 7.425316 -0.825316477
## 823 8.6 7.464168 1.135832409
## 824 12.0 13.697696 -1.697695889
## 825 5.8 7.440425 -1.640425195
## 826 7.8 7.088602 0.711397907
## 827 10.4 9.989524 0.410475858
## 828 14.2 15.849645 -1.649644857
## 829 5.0 9.387324 -4.387323532
## 830 7.8 11.170182 -3.370181643
## 831 7.0 6.853334 0.146665942
## 832 7.0 8.979382 -1.979381569
## 833 8.0 7.505178 0.494821963
## 834 6.8 6.533887 0.266112857
## 835 19.6 12.454444 7.145556642
## 836 18.0 14.986275 3.013724848
## 837 14.2 12.050819 2.149181087
## 838 15.2 12.374582 2.825417519
## 839 5.0 6.372006 -1.372005563
## 840 17.6 15.202118 2.397882851
## 841 15.2 13.186150 2.013850127
## 842 24.8 17.360541 7.439457921
## 843 8.8 14.498471 -5.698470725
## 844 6.0 9.048451 -3.048451064
## 845 8.0 9.721879 -1.721879260
## 846 6.0 7.328187 -1.328187250
## 847 9.0 7.541871 1.458128742
## 848 5.6 7.967081 -2.367080794
## 849 9.4 9.527621 -0.127621583
## 850 7.0 12.394008 -5.394007725
## 851 8.0 6.680660 1.319339834
## 852 7.4 9.264293 -1.864293347
## 853 16.4 10.934914 5.465086010
## 854 12.0 9.212491 2.787508679
## 855 5.8 7.787932 -1.987931350
## 856 14.8 9.521146 5.278854408
## 857 4.0 8.152705 -4.152705276
## 858 8.0 6.253292 1.746707748
## 859 4.4 8.072843 -3.672843319
## 860 29.6 11.532797 18.067203697
## 861 11.2 9.657127 1.542873098
## 862 11.0 6.898661 4.101339034
## 863 4.2 11.316954 -7.116954497
## 864 6.8 9.370056 -2.570055968
## 865 13.8 14.261045 -0.461044438
## 866 14.0 14.986275 -0.986275152
## 867 5.2 7.852684 -2.652684486
## 868 5.4 6.510145 -1.110144599
## 869 7.4 6.393590 1.006410290
## 870 14.4 13.963182 0.436817703
## 871 8.2 6.449709 1.750291000
## 872 6.2 7.846209 -1.646209169
## 873 7.0 8.947005 -1.947005089
## 874 22.6 12.635751 9.964249011
## 875 12.8 15.199959 -2.399958455
## 876 7.7 7.995140 -0.295140423
## 877 6.6 7.459851 -0.859851225
## 878 7.0 7.077810 -0.077810124
## 879 5.2 8.273577 -3.073577083
## 880 5.8 7.421000 -1.620999348
## 881 7.0 7.764189 -0.764188953
## 882 5.9 7.906645 -2.006644795
## 883 9.8 6.516620 3.283380231
## 884 7.0 10.147089 -3.147088598
## 885 8.4 6.905136 1.494863388
## 886 4.4 8.981540 -4.581539741
## 887 22.0 15.204276 6.795724408
## 888 5.2 12.057294 -6.857293921
## 889 6.8 7.904486 -1.104486226
## 890 4.2 7.921754 -3.721753981
## 891 9.0 9.354947 -0.354947258
## 892 8.2 12.702662 -4.502662583ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "blue",se = F)+
theme_bw()
Regresi Plinomial Derajat 2 (ordo 2)
mod_polinomial2 = lm(triceps ~ poly(age,2,raw = T),
data=triceps)
summary(mod_polinomial2)##
## Call:
## lm(formula = triceps ~ poly(age, 2, raw = T), data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.5677 -2.4401 -0.4587 1.6368 24.9961
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0229191 0.3063806 19.658 < 2e-16 ***
## poly(age, 2, raw = T)1 0.2434733 0.0364403 6.681 4.17e-11 ***
## poly(age, 2, raw = T)2 -0.0006257 0.0007926 -0.789 0.43
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.008 on 889 degrees of freedom
## Multiple R-squared: 0.3377, Adjusted R-squared: 0.3362
## F-statistic: 226.6 on 2 and 889 DF, p-value: < 2.2e-16AIC(mod_polinomial2)## [1] 5012.89ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = F)+
theme_bw()
ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = T)+
theme_bw()
Regresi Polinomial Derajat 3 (ordo 3)
mod_polinomial3 = lm(triceps ~ poly(age,3,raw = T),data=triceps)
summary(mod_polinomial3)##
## Call:
## lm(formula = triceps ~ poly(age, 3, raw = T), data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.5832 -1.9284 -0.5415 1.3283 24.4440
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.004e+00 3.831e-01 20.893 < 2e-16 ***
## poly(age, 3, raw = T)1 -3.157e-01 7.721e-02 -4.089 4.73e-05 ***
## poly(age, 3, raw = T)2 3.101e-02 3.964e-03 7.824 1.45e-14 ***
## poly(age, 3, raw = T)3 -4.566e-04 5.612e-05 -8.135 1.38e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.868 on 888 degrees of freedom
## Multiple R-squared: 0.3836, Adjusted R-squared: 0.3815
## F-statistic: 184.2 on 3 and 888 DF, p-value: < 2.2e-16ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,3,raw=T),
lty = 1, col = "blue",se = T)+
theme_bw()
AIC <- cbind(AIC(mod_linear),AIC(mod_polinomial2),AIC(mod_polinomial3))MSE(predict(mod_linear),triceps$triceps)## [1] 16.01758MSE(predict(mod_polinomial2),triceps$triceps)## [1] 16.00636MSE(predict(mod_polinomial3),triceps$triceps)## [1] 14.89621Regresi Fungsi Tangga (5)
mod_tangga5 = lm(triceps ~ cut(age,5),data=triceps)
summary(mod_tangga5)##
## Call:
## lm(formula = triceps ~ cut(age, 5), data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5474 -2.0318 -0.4465 1.3682 23.3759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.2318 0.1994 36.260 < 2e-16 ***
## cut(age, 5)(10.6,20.9] 1.6294 0.3244 5.023 6.16e-07 ***
## cut(age, 5)(20.9,31.2] 5.9923 0.4222 14.192 < 2e-16 ***
## cut(age, 5)(31.2,41.5] 7.5156 0.4506 16.678 < 2e-16 ***
## cut(age, 5)(41.5,51.8] 7.4561 0.5543 13.452 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.939 on 887 degrees of freedom
## Multiple R-squared: 0.3617, Adjusted R-squared: 0.3588
## F-statistic: 125.7 on 4 and 887 DF, p-value: < 2.2e-16ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "blue",se = F)+
theme_bw()
Regresi FUngsi Tangga (7)
mod_tangga7 = lm(triceps ~ cut(age,7),data=triceps)
summary(mod_tangga7)##
## Call:
## lm(formula = triceps ~ cut(age, 7), data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8063 -1.7592 -0.4366 1.2894 23.1461
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.5592 0.2219 34.060 < 2e-16 ***
## cut(age, 7)(7.62,15] -0.6486 0.3326 -1.950 0.0515 .
## cut(age, 7)(15,22.3] 3.4534 0.4239 8.146 1.27e-15 ***
## cut(age, 7)(22.3,29.7] 5.8947 0.4604 12.804 < 2e-16 ***
## cut(age, 7)(29.7,37] 7.8471 0.5249 14.949 < 2e-16 ***
## cut(age, 7)(37,44.4] 6.9191 0.5391 12.835 < 2e-16 ***
## cut(age, 7)(44.4,51.8] 6.3013 0.6560 9.606 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.805 on 885 degrees of freedom
## Multiple R-squared: 0.4055, Adjusted R-squared: 0.4014
## F-statistic: 100.6 on 6 and 885 DF, p-value: < 2.2e-16ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,7),
lty = 1, col = "blue",se = F)+
theme_bw()
Perbandingan Model
MSE = function(pred,actual){
mean((pred-actual)^2)
}Membandingkan model-modelnya
nilai_MSE <- rbind(MSE(predict(mod_linear),triceps$triceps),
MSE(predict(mod_polinomial2),triceps$triceps),
MSE(predict(mod_polinomial3),triceps$triceps),
MSE(predict(mod_tangga5),triceps$triceps),
MSE(predict(mod_tangga7),triceps$triceps))
nama_model <- c("Linear","Poly (ordo=2)", "Poly (ordo=3)",
"Tangga (breaks=5)", "Tangga (breaks=7)")
data.frame(nama_model,nilai_MSE)## nama_model nilai_MSE
## 1 Linear 16.01758
## 2 Poly (ordo=2) 16.00636
## 3 Poly (ordo=3) 14.89621
## 4 Tangga (breaks=5) 15.42602
## 5 Tangga (breaks=7) 14.36779Evaluasi Model Menggunakan Cross Validation
Regresi Linear
set.seed(1401201089)
cross_val <- vfold_cv(triceps,v=10,strata = "triceps")
metric_linear <- map_dfr(cross_val$splits,
function(x){
mod <- lm(triceps ~ age,data=triceps[x$in_id,])
pred <- predict(mod,newdata=triceps[-x$in_id,])
truth <- triceps[-x$in_id,]$triceps
rmse <- mlr3measures::rmse(truth = truth,response = pred)
mae <- mlr3measures::mae(truth = truth,response = pred)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_linear## # A tibble: 10 × 2
## rmse mae
## <dbl> <dbl>
## 1 3.69 2.82
## 2 4.14 2.99
## 3 4.45 3.11
## 4 3.16 2.29
## 5 4.15 2.85
## 6 3.63 2.63
## 7 4.35 2.96
## 8 3.92 2.70
## 9 3.86 2.88
## 10 4.55 3.16# menghitung rata-rata untuk 10 folds
mean_metric_linear <- colMeans(metric_linear)
mean_metric_linear## rmse mae
## 3.989941 2.838418Polynomial Derajat 2 (ordo 2)
set.seed(1401201089)
cross_val <- vfold_cv(triceps,v=10,strata = "triceps")
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(triceps ~ poly(age,2,raw = T),data=triceps[x$in_id,])
pred <- predict(mod,newdata=triceps[-x$in_id,])
truth <- triceps[-x$in_id,]$triceps
rmse <- mlr3measures::rmse(truth = truth,response = pred)
mae <- mlr3measures::mae(truth = truth,response = pred)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_poly2## # A tibble: 10 × 2
## rmse mae
## <dbl> <dbl>
## 1 3.68 2.83
## 2 4.20 3.07
## 3 4.45 3.11
## 4 3.18 2.30
## 5 4.14 2.86
## 6 3.62 2.64
## 7 4.35 2.97
## 8 3.92 2.71
## 9 3.85 2.90
## 10 4.55 3.17# menghitung rata-rata untuk 10 folds
mean_metric_poly2 <- colMeans(metric_poly2)
mean_metric_poly2## rmse mae
## 3.994084 2.856381Polynomial Derajat 3 (ordo 3)
set.seed(1401201089)
cross_val <- vfold_cv(triceps,v=10,strata = "triceps")
metric_poly3 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(triceps ~ poly(age,3,raw = T),data=triceps[x$in_id,])
pred <- predict(mod,newdata=triceps[-x$in_id,])
truth <- triceps[-x$in_id,]$triceps
rmse <- mlr3measures::rmse(truth = truth,response = pred)
mae <- mlr3measures::mae(truth = truth,response = pred)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_poly3## # A tibble: 10 × 2
## rmse mae
## <dbl> <dbl>
## 1 3.39 2.58
## 2 4.09 2.71
## 3 4.30 2.84
## 4 3.19 2.18
## 5 4.12 2.72
## 6 3.50 2.52
## 7 4.22 2.78
## 8 3.62 2.44
## 9 3.70 2.55
## 10 4.42 3.03# menghitung rata-rata untuk 10 folds
mean_metric_poly3 <- colMeans(metric_poly3)
mean_metric_poly3## rmse mae
## 3.855217 2.634753Fungsi Tangga
set.seed(1401201089)
cross_val <- vfold_cv(triceps,v=10,strata = "triceps")
breaks <- 3:10
best_tangga <- map_dfr(breaks, function(i){
metric_tangga <- map_dfr(cross_val$splits,
function(x){
training <- triceps[x$in_id,]
training$age <- cut(training$age,i)
mod <- lm(triceps ~ age,data=training)
labs_x <- levels(mod$model[,2])
labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))
testing <- triceps[-x$in_id,]
age_new <- cut(testing$age,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
pred <- predict(mod,newdata=list(age=age_new))
truth <- testing$triceps
data_eval <- na.omit(data.frame(truth,pred))
rmse <- mlr3measures::rmse(truth = data_eval$truth,response = data_eval$pred)
mae <- mlr3measures::mae(truth = data_eval$truth,response = data_eval$pred)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_tangga
# menghitung rata-rata untuk 10 folds
mean_metric_tangga <- colMeans(metric_tangga)
mean_metric_tangga
}
)
best_tangga <- cbind(breaks=breaks,best_tangga)
# menampilkan hasil all breaks
best_tangga## breaks rmse mae
## 1 3 3.837104 2.613577
## 2 4 3.905216 2.660754
## 3 5 3.938642 2.729484
## 4 6 3.835198 2.622106
## 5 7 3.801116 2.560375
## 6 8 3.820517 2.559213
## 7 9 3.792659 2.514480
## 8 10 3.812491 2.532252#berdasarkan rmse
best_tangga %>% slice_min(rmse)## breaks rmse mae
## 1 9 3.792659 2.51448#berdasarkan mae
best_tangga %>% slice_min(mae)## breaks rmse mae
## 1 9 3.792659 2.51448Perbandingan Model
nilai_metric <- rbind(mean_metric_linear,
mean_metric_poly2,
mean_metric_poly3,
best_tangga %>% select(-1) %>% slice_min(mae))
nama_model <- c("Linear","Poly2","Poly3","Tangga (breaks=9)")
data.frame(nama_model,nilai_metric)## nama_model rmse mae
## 1 Linear 3.989941 2.838418
## 2 Poly2 3.994084 2.856381
## 3 Poly3 3.855217 2.634753
## 4 Tangga (breaks=9) 3.792659 2.514480PERTEMUAN 9
Untuk tugas bagian pertemuan 9 ini saya akan melakukan interpretasi pada data bangkitan set.seed 1401201089 saja, sedangkan pada data triceps tidak saya interpretasikan karena isinya kurang lebih sama. ## Library
library(tidyverse)
library(ggplot2)
library(dplyr)
library(purrr)
library(rsample)
library(splines)#Data Data yang digunakan yaitu data yang dibangkitkan dari fungsi set.seed 1401201089 (NIM saya). Dibangkitkan 1000 data.x yang menyebar normal dengan rata-rata dan ragam sama dengan 1. Model yang digunakan pada data tersebut yaitu Model Polynomial Ordo 2 dengan persamaan sebagai berikut: \[y=5+2(data.x)+3(data.x)^2+error\]
set.seed(1401201089)
data.x <- rnorm(1000,1,1)
err <- rnorm(1000,0,5)
y <- 5+2*data.x+3*data.x^2+err
plot(data.x,y,xlim=c(-3,5), ylim=c(-10,70), pch=16, col="orange")
abline(v=1, col="red", lty=2)
Pada plot di atas, kita memberi titik potong (knot) di x=1.
Piecewise Cubic Polynomial
dt.all <- cbind(y,data.x)
##knots=1
dt1 <- dt.all[data.x <=1,]
dim(dt1)## [1] 503 2dt2 <- dt.all[data.x >1,]
dim(dt2)## [1] 497 2maksud dari knots = 1 yaitu titik potongnya berada pada titik 1. Dimensi data yang berada di samping kiri titik potong (<1) yaitu sebanyak 503 data, sedangkan dimensi data yang berada di kanan titik potong (>1) yaitu sebanyak 947 data.
plot(data.x,y,xlim=c(-3,5), ylim=c(-10,70), pch=16, col="orange")
cub.mod1 <- lm(dt1[,1]~dt1[,2]+I(dt1[,2]^2)+I(dt1[,2]^3))
ix <- sort(dt1[,2], index.return=T)$ix
lines(dt1[ix,2],cub.mod1$fitted.values[ix],col="blue", lwd=2)
cub.mod2 <- lm(dt2[,1]~dt2[,2]+I(dt2[,2]^2)+I(dt2[,2]^3))
ix <- sort(dt2[,2], index.return=T)$ix
lines(dt2[ix,2],cub.mod2$fitted.values[ix],col="blue", lwd=2)
garis biru di atas merupakan plot Piecewise Cubic Polynomial antara dt1 dan dt2. Terlihat patahan (cutoff) pada titik knots = 1.
Truncated Power Basis
#dengan menggunakan truncated power basis
plot(data.x,y,xlim=c( -3,5), ylim=c( -10,70), pch=16,col="orange")
abline(v=1,col="red", lty=2)
hx <- ifelse( data.x>1,(data.x-1)^3,0)
cubspline.mod <- lm( y~data.x+I(data.x^2)+I(data.x^3)+hx)
ix <- sort( data.x,index.return=T)$ix
lines( data.x[ix], cubspline.mod$fitted.values[ix],col="green", lwd=2)
pada plot Truncated Power Basis, plot yang pada awalnya memiliki cutoff
pada titik 1 dihubungkan menjadi satu kesatuan.
Perbandingan dengan k-fold CV
plot( data.x,y,xlim=c(-3,5), ylim=c( -10,70), pch=16,col="orange")
abline(v=0,col="red", lty=2)
abline(v=2,col="red", lty=2)
hx1 <- ifelse( data.x>0,(data.x-0)^3,0)
hx2 <- ifelse( data.x>2,(data.x-2)^3,0)
cubspline.mod2 <- lm( y~data.x+I(data.x^2)+I(data.x^3)+hx1+hx2)
ix <- sort( data.x,index.return=T)$ix
lines( data.x[ix],cubspline.mod2$fitted.values[ix],col="green", lwd=2)
Sekarang kita membuat knot pada 2 titik, yaitu pada titik 0 dan titik 2. Maka kita memiliki 3 fungsi, yaitu kurang dari 0, dari 0 sampai 2, dan lebih dari 2.
Perbandingan 1 knot dan 2 knot dengan 5-fold CV
Kita bandingkan nilai residunya
##1 knot (knot=1)
data.all <- cbind( y,data.x,hx)
set.seed(456)
ind <- sample(1:5,length( data.x),replace=T)
res <- c()
for( i in 1:5){
dt.train <- data.all[ ind!=i,]
x.test <- data.all[ ind==i, -1]
y.test <- data.all[ ind==i,1]
mod1 <- lm( dt.train[,1]~dt.train[,2]+I( dt.train[,2]^2)+I( dt.train[,2]^3)+dt.train[,3])
x.test.olah <- cbind(1,x.test[,1], x.test[,1]^2,x.test[,1]^3,x.test[,2])
beta <- coefficients(mod1)
prediksi <- x.test.olah%*%beta
res <- c( res,mean(( y.test-prediksi)^2))
}
res## [1] 27.23126 20.16745 24.23764 21.42984 23.41967mean(res)## [1] 23.29717##2 knot (knot = 0, 2)
data.all2 <- cbind(y,data.x,hx1,hx2)
set.seed(456)
ind2 <- sample(1:5,length( data.x),replace=T)
res2 <- c()
for( i in 1:5){
dt.train2 <- data.all2[ind2!=i,]
x.test2 <- data.all2[ind2==i,-1]
y.test2 <- data.all2[ind2==i,1]
mod2 <- lm(dt.train2[,1]~dt.train2[,2]+I(dt.train2[,2]^2)+I(dt.train2[,2]^3)+dt.train2[,3]+dt.train2[,4])
x.test.olah2 <- cbind(1,x.test2[,1],x.test2[,1]^2,x.test2[,1]^3,x.test2[,2],x.test2[,3])
beta2 <- coefficients(mod2)
prediksi2 <- x.test.olah2%*%beta2
res2 <- c(res2,mean((y.test2-prediksi2)^2))
}
res2## [1] 27.24392 20.19229 24.28500 21.49835 23.42328mean(res2)## [1] 23.32857Dari perbandingan residu perlipatan, 1 knot lebih bagus daripada 2 knot karena nilai residunya lebih kecil.
Smoothing Spline
Fungsi smooth.spline secara otomatis memilih paramter lambda terbaik dengan menggunakan Cross-Validation (CV) dan Generalized Cross-Validation (GCV). Perbedaan utama antara keduanya adalah pada ukuran kebaikan model yang digunakan. Kalau CV menggunakan MSE sebagai ukuran kebaikan model, sedangkan GCV menggunakan Weighted-MSE. Banyaknya folds (lipatan) pada CV dan GCV ini adalah sejumlah obserbasinya.
ss1 <- smooth.spline(data.x,y,all.knots = T)
plot(data.x,y,xlim=c(-3,5), ylim=c(-10,70), pch=16,col="orange")
lines(ss1,col="purple", lwd=2)
Contoh Data TRICEPS
Data yang digunakan untuk ilustrasi berasal dari studi antropometri terhadap 892 perempuan di bawah 50 tahun di tiga desa Gambia di Afrika Barat. Data terdiri dari 3 Kolom yaitu Age, Intriceps dan tricepts. Berikut adalah penjelasannya pada masing-masing kolom:
- age : umur responden
- Intriceps : logaritma dari ketebalan lipatan kulit triceps
- triceps: ketebalan lipatan kulit triceps
Lipatan kulit trisep diperlukan untuk menghitung lingkar otot lengan atas. Ketebalannya memberikan informasi tentang cadangan lemak tubuh, sedangkan massa otot yang dihitung memberikan informasi tentang cadangan protein
data("triceps", package="MultiKink")kita gambarkan dalam bentuk scatterplot
ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
theme_bw()
Regresi SPline
#Menentukan banyaknya fungsi basis dan knots
dim(bs(triceps$age, knots = c(10, 20,40)))## [1] 892 6 #atau
dim(bs(triceps$age, df=6))## [1] 892 6#nilai knots yang ditentukan oleh komputer
attr(bs(triceps$age, df=6),"knots")## 25% 50% 75%
## 5.5475 12.2100 24.7275b-spline
Menggunakan knot yang ditentukan manual
knot_value_manual_3 = c(10, 20,40)
mod_spline_1 = lm(triceps ~ bs(age, knots =knot_value_manual_3 ),data=triceps)
summary(mod_spline_1)##
## Call:
## lm(formula = triceps ~ bs(age, knots = knot_value_manual_3),
## data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.385 -1.682 -0.393 1.165 22.855
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.22533 0.61873 13.294 < 2e-16 ***
## bs(age, knots = knot_value_manual_3)1 -0.06551 1.14972 -0.057 0.955
## bs(age, knots = knot_value_manual_3)2 -4.30051 0.72301 -5.948 3.90e-09 ***
## bs(age, knots = knot_value_manual_3)3 7.80435 1.17793 6.625 6.00e-11 ***
## bs(age, knots = knot_value_manual_3)4 6.14890 1.27439 4.825 1.65e-06 ***
## bs(age, knots = knot_value_manual_3)5 5.56640 1.42225 3.914 9.78e-05 ***
## bs(age, knots = knot_value_manual_3)6 7.90178 1.54514 5.114 3.87e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.76 on 885 degrees of freedom
## Multiple R-squared: 0.4195, Adjusted R-squared: 0.4155
## F-statistic: 106.6 on 6 and 885 DF, p-value: < 2.2e-16ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~bs(x, knots = knot_value_manual_3),
lty = 1,se = F)
Menggunakan knot yang ditentukan fungsi komputer
knot_value_pc_df_6 = attr(bs(triceps$age, df=6),"knots")
mod_spline_1 = lm(triceps ~ bs(age, knots =knot_value_pc_df_6 ),data=triceps)
summary(mod_spline_1)##
## Call:
## lm(formula = triceps ~ bs(age, knots = knot_value_pc_df_6), data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.0056 -1.7556 -0.2944 1.2008 23.0695
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.5925 0.8770 7.517 1.37e-13 ***
## bs(age, knots = knot_value_pc_df_6)1 3.7961 1.5015 2.528 0.01164 *
## bs(age, knots = knot_value_pc_df_6)2 -2.0749 0.8884 -2.335 0.01974 *
## bs(age, knots = knot_value_pc_df_6)3 1.5139 1.1645 1.300 0.19391
## bs(age, knots = knot_value_pc_df_6)4 11.6394 1.3144 8.855 < 2e-16 ***
## bs(age, knots = knot_value_pc_df_6)5 5.9680 1.5602 3.825 0.00014 ***
## bs(age, knots = knot_value_pc_df_6)6 8.9127 1.4053 6.342 3.60e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.757 on 885 degrees of freedom
## Multiple R-squared: 0.4206, Adjusted R-squared: 0.4167
## F-statistic: 107.1 on 6 and 885 DF, p-value: < 2.2e-16ggplot(triceps,aes(x=age, y=triceps)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~bs(x, knots = knot_value_pc_df_6),
lty = 1,se = F)
Naturan Spline
knot_value_manual_3 = c(10, 20,40)
mod_spline3ns = lm(triceps ~ ns(age, knots = knot_value_manual_3),data=triceps)
summary(mod_spline3ns)##
## Call:
## lm(formula = triceps ~ ns(age, knots = knot_value_manual_3),
## data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7220 -1.7640 -0.3985 1.1908 22.9684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.8748 0.3742 23.714 < 2e-16 ***
## ns(age, knots = knot_value_manual_3)1 7.0119 0.6728 10.422 < 2e-16 ***
## ns(age, knots = knot_value_manual_3)2 6.0762 0.8625 7.045 3.72e-12 ***
## ns(age, knots = knot_value_manual_3)3 2.0780 1.0632 1.954 0.051 .
## ns(age, knots = knot_value_manual_3)4 8.8616 0.9930 8.924 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.762 on 887 degrees of freedom
## Multiple R-squared: 0.4176, Adjusted R-squared: 0.415
## F-statistic: 159 on 4 and 887 DF, p-value: < 2.2e-16knot_value_manual_3 = c(10, 20,40)
mod_spline3ns = lm(triceps ~ ns(age, knots = knot_value_manual_3),data=triceps)
summary(mod_spline3ns)##
## Call:
## lm(formula = triceps ~ ns(age, knots = knot_value_manual_3),
## data = triceps)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7220 -1.7640 -0.3985 1.1908 22.9684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.8748 0.3742 23.714 < 2e-16 ***
## ns(age, knots = knot_value_manual_3)1 7.0119 0.6728 10.422 < 2e-16 ***
## ns(age, knots = knot_value_manual_3)2 6.0762 0.8625 7.045 3.72e-12 ***
## ns(age, knots = knot_value_manual_3)3 2.0780 1.0632 1.954 0.051 .
## ns(age, knots = knot_value_manual_3)4 8.8616 0.9930 8.924 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.762 on 887 degrees of freedom
## Multiple R-squared: 0.4176, Adjusted R-squared: 0.415
## F-statistic: 159 on 4 and 887 DF, p-value: < 2.2e-16Smoothing Splin
Fungsi smooth.spline() dapat digunakan untuk melakukan smoothing spline pada R
model_sms <- with(data = triceps,smooth.spline(age,triceps))
model_sms ## Call:
## smooth.spline(x = age, y = triceps)
##
## Smoothing Parameter spar= 0.5162704 lambda= 1.232259e-06 (13 iterations)
## Equivalent Degrees of Freedom (Df): 50.00639
## Penalized Criterion (RSS): 8591.581
## GCV: 13.77835Fungsi smooth.spline secara otomatis memilih paramter lambda terbaik dengan menggunakan Cross-Validation (CV) dan Generalized Cross-Validation (GCV). Perbedaan utama antara keduanya adalah pada ukuran kebaikan model yang digunakan. Kalau CV menggunakan MSE sebagai ukuran kebaikan model, sedangkan GCV menggunakan Weighted-MSE. Banyaknya folds (lipatan) pada CV dan GCV ini adalah sejumlah obserbasinya.
pred_data <- broom::augment(model_sms)
ggplot(pred_data,aes(x=x,y=y))+
geom_point(alpha=0.55, color="black")+
geom_line(aes(y=.fitted),col="blue",
lty=1)+
xlab("age")+
ylab("triceps")+
theme_bw()
Pengaruh parameter lambda terhadap tingkat smoothness kurva bisa dilihat dari ilustrasi berikut:
model_sms_lambda <- data.frame(lambda=seq(0,5,by=0.5)) %>%
group_by(lambda) %>%
do(broom::augment(with(data = triceps,smooth.spline(age,triceps,lambda = .$lambda))))
p <- ggplot(model_sms_lambda,
aes(x=x,y=y))+
geom_line(aes(y=.fitted),
col="blue",
lty=1
)+
facet_wrap(~lambda)
p
Jika kita tentukan df=7, maka hasil kurva model smooth.spline akan lebih merepresentasikan data
model_sms <- with(data = triceps,smooth.spline(age,triceps,df=7))
model_sms ## Call:
## smooth.spline(x = age, y = triceps, df = 7)
##
## Smoothing Parameter spar= 1.049874 lambda= 0.008827582 (11 iterations)
## Equivalent Degrees of Freedom (Df): 7.000747
## Penalized Criterion (RSS): 10112.16
## GCV: 14.20355pred_data <- broom::augment(model_sms)
ggplot(pred_data,aes(x=x,y=y))+
geom_point(alpha=0.55, color="black")+
geom_line(aes(y=.fitted),col="blue",
lty=1)+
xlab("age")+
ylab("triceps")+
theme_bw()
LOESS
Masih menggunakan data triceps seperti pada ilustrasi sebelumnya, kali ini kita akan mencoba melakukan pendekatan local regression dengan fungsi loess().
model_loess <- loess(triceps~age,
data = triceps)
summary(model_loess)## Call:
## loess(formula = triceps ~ age, data = triceps)
##
## Number of Observations: 892
## Equivalent Number of Parameters: 4.6
## Residual Standard Error: 3.777
## Trace of smoother matrix: 5.01 (exact)
##
## Control settings:
## span : 0.75
## degree : 2
## family : gaussian
## surface : interpolate cell = 0.2
## normalize: TRUE
## parametric: FALSE
## drop.square: FALSEPengaruh parameter span terhadap tingkat smoothness kurva bisa dilihat dari ilustrasi berikut:
model_loess_span <- data.frame(span=seq(0.1,5,by=0.5)) %>%
group_by(span) %>%
do(broom::augment(loess(triceps~age,
data = triceps,span=.$span)))
p2 <- ggplot(model_loess_span,
aes(x=age,y=triceps))+
geom_line(aes(y=.fitted),
col="blue",
lty=1
)+
facet_wrap(~span)
p2
Pendekatan ini juga dapat dilakukan dengan fungsi stat_smooth() pada package ggplot2.
library(ggplot2)
ggplot(triceps, aes(age,triceps)) +
geom_point(alpha=0.5,color="black") +
stat_smooth(method='loess',
formula=y~x,
span = 0.75,
col="blue",
lty=1,
se=F)
Tuning span dapat dilakukan dengan menggunakan cross-validation:
set.seed(123)
cross_val <- vfold_cv(triceps,v=10,strata = "triceps")
span <- seq(0.1,1,length.out=50)
best_loess <- map_dfr(span, function(i){
metric_loess <- map_dfr(cross_val$splits,
function(x){
mod <- loess(triceps ~ age,span = i,
data=triceps[x$in_id,])
pred <- predict(mod,
newdata=triceps[-x$in_id,])
truth <- triceps[-x$in_id,]$triceps
data_eval <- na.omit(data.frame(pred=pred,
truth=truth))
rmse <- mlr3measures::rmse(truth = data_eval$truth,
response = data_eval$pred
)
mae <- mlr3measures::mae(truth = data_eval$truth,
response = data_eval$pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
head(metric_loess,20)
# menghitung rata-rata untuk 10 folds
mean_metric_loess <- colMeans(metric_loess)
mean_metric_loess
}
)
best_loess <- cbind(span=span,best_loess)
# menampilkan hasil all breaks
best_loess## span rmse mae
## 1 0.1000000 3.773822 2.495752
## 2 0.1183673 3.771342 2.493506
## 3 0.1367347 3.766450 2.487557
## 4 0.1551020 3.757385 2.480571
## 5 0.1734694 3.749539 2.473494
## 6 0.1918367 3.746288 2.470915
## 7 0.2102041 3.742741 2.467990
## 8 0.2285714 3.740397 2.465388
## 9 0.2469388 3.737724 2.463311
## 10 0.2653061 3.737105 2.462252
## 11 0.2836735 3.735645 2.460514
## 12 0.3020408 3.734522 2.459225
## 13 0.3204082 3.733298 2.457694
## 14 0.3387755 3.732538 2.456280
## 15 0.3571429 3.732082 2.455452
## 16 0.3755102 3.731261 2.454642
## 17 0.3938776 3.730882 2.454293
## 18 0.4122449 3.730585 2.454291
## 19 0.4306122 3.730351 2.454869
## 20 0.4489796 3.730717 2.456222
## 21 0.4673469 3.730831 2.456768
## 22 0.4857143 3.731387 2.458307
## 23 0.5040816 3.732091 2.460189
## 24 0.5224490 3.732392 2.461811
## 25 0.5408163 3.733609 2.464543
## 26 0.5591837 3.734393 2.467180
## 27 0.5775510 3.735554 2.470517
## 28 0.5959184 3.736679 2.473012
## 29 0.6142857 3.738227 2.476519
## 30 0.6326531 3.741401 2.476887
## 31 0.6510204 3.742660 2.479842
## 32 0.6693878 3.744336 2.483581
## 33 0.6877551 3.746651 2.488792
## 34 0.7061224 3.748708 2.492180
## 35 0.7244898 3.750734 2.495218
## 36 0.7428571 3.752815 2.498359
## 37 0.7612245 3.753998 2.504175
## 38 0.7795918 3.755981 2.511421
## 39 0.7979592 3.756794 2.514375
## 40 0.8163265 3.759456 2.522994
## 41 0.8346939 3.767892 2.543250
## 42 0.8530612 3.777168 2.557863
## 43 0.8714286 3.782407 2.566972
## 44 0.8897959 3.788428 2.578042
## 45 0.9081633 3.797002 2.595698
## 46 0.9265306 3.803259 2.609852
## 47 0.9448980 3.809478 2.623386
## 48 0.9632653 3.824192 2.655691
## 49 0.9816327 3.835596 2.679155
## 50 1.0000000 3.860514 2.718911#berdasarkan rmse
best_loess %>% slice_min(rmse)## span rmse mae
## 1 0.4306122 3.730351 2.454869#berdasarkan mae
best_loess %>% slice_min(mae)## span rmse mae
## 1 0.4122449 3.730585 2.454291library(ggplot2)
ggplot(triceps, aes(age,triceps)) +
geom_point(alpha=0.5,color="black") +
stat_smooth(method='loess',
formula=y~x,
span = 0.4306122,
col="blue",
lty=1,
se=F)
LATIHAN
library(ISLR)## Warning: package 'ISLR' was built under R version 4.2.2AutoData = Auto %>% select(mpg, horsepower, origin)
tibble(AutoData)## # A tibble: 392 × 3
## mpg horsepower origin
## <dbl> <dbl> <dbl>
## 1 18 130 1
## 2 15 165 1
## 3 18 150 1
## 4 16 150 1
## 5 17 140 1
## 6 15 198 1
## 7 14 220 1
## 8 14 215 1
## 9 14 225 1
## 10 15 190 1
## # … with 382 more rowsplot(AutoData$mpg, AutoData$horsepower,col="black")
Dari grafik di atas, secara ekploratif, data terlihat sedikit melengkung sehingga tidak membentuk pola linier.
Regresi Linear
regresi linear merupakan suatu pendekatan untuk mengetahui hubungan antara satu atau lebih variabel dependen dan juga variabel independen.
lin.mod2 <-lm(AutoData$horsepower~AutoData$mpg)
plot(AutoData$mpg,AutoData$horsepower)
lines(AutoData$mpg,lin.mod2$fitted.values,col="blue")
Dari grafik di atas, dapat diketahui bahwa garis regresi linier tidak mendekati pola data, sehingga tidak akan tepat apabila digunakan fungsi liner biasa pada data ini.
summary(lin.mod2)##
## Call:
## lm(formula = AutoData$horsepower ~ AutoData$mpg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -64.892 -15.716 -2.094 13.108 96.947
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 194.4756 3.8732 50.21 <2e-16 ***
## AutoData$mpg -3.8389 0.1568 -24.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 24.19 on 390 degrees of freedom
## Multiple R-squared: 0.6059, Adjusted R-squared: 0.6049
## F-statistic: 599.7 on 1 and 390 DF, p-value: < 2.2e-16Polynomial
pol.mod2 <- lm(AutoData$horsepower~AutoData$mpg+I(AutoData$mpg^2)) #ordo 2
ix2 <- sort(AutoData$mpg,index.return=T)$ix
plot(AutoData$mpg,AutoData$horsepower)
lines(AutoData$mpg[ix2], pol.mod2$fitted.values[ix2],col="blue", cex=2)
Garis biru berupakan hasil fungsi polynomial orde 2. Dapat dilihat bahwa garis tersebut mendekati observasi yang ada, sehingga errornya menjadi lebih kecil daripada fungsi regresi linier sebelumnya.
summary(pol.mod2)##
## Call:
## lm(formula = AutoData$horsepower ~ AutoData$mpg + I(AutoData$mpg^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -72.52 -11.24 -0.11 9.47 92.98
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 302.06824 9.25689 32.63 <2e-16 ***
## AutoData$mpg -13.32406 0.77456 -17.20 <2e-16 ***
## I(AutoData$mpg^2) 0.18804 0.01513 12.43 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.49 on 389 degrees of freedom
## Multiple R-squared: 0.718, Adjusted R-squared: 0.7165
## F-statistic: 495.1 on 2 and 389 DF, p-value: < 2.2e-16Fungsi Tangga
range(AutoData$mpg)## [1] 9.0 46.6data tersebut memiliki nilai x minimum 9 dan nilai maksimum 46.6
d1 <- as.factor(ifelse(AutoData$mpg<=20,1,0))
d2 <- as.factor(ifelse(AutoData$mpg<=30 & AutoData$mpg>20,1,0))
d3 <- as.factor(ifelse(AutoData$mpg>30,1,0))
data.frame(AutoData$horsepower,d1,d2,d3)## AutoData.horsepower d1 d2 d3
## 1 130 1 0 0
## 2 165 1 0 0
## 3 150 1 0 0
## 4 150 1 0 0
## 5 140 1 0 0
## 6 198 1 0 0
## 7 220 1 0 0
## 8 215 1 0 0
## 9 225 1 0 0
## 10 190 1 0 0
## 11 170 1 0 0
## 12 160 1 0 0
## 13 150 1 0 0
## 14 225 1 0 0
## 15 95 0 1 0
## 16 95 0 1 0
## 17 97 1 0 0
## 18 85 0 1 0
## 19 88 0 1 0
## 20 46 0 1 0
## 21 87 0 1 0
## 22 90 0 1 0
## 23 95 0 1 0
## 24 113 0 1 0
## 25 90 0 1 0
## 26 215 1 0 0
## 27 200 1 0 0
## 28 210 1 0 0
## 29 193 1 0 0
## 30 88 0 1 0
## 31 90 0 1 0
## 32 95 0 1 0
## 33 100 1 0 0
## 34 105 1 0 0
## 35 100 1 0 0
## 36 88 1 0 0
## 37 100 1 0 0
## 38 165 1 0 0
## 39 175 1 0 0
## 40 153 1 0 0
## 41 150 1 0 0
## 42 180 1 0 0
## 43 170 1 0 0
## 44 175 1 0 0
## 45 110 1 0 0
## 46 72 0 1 0
## 47 100 1 0 0
## 48 88 1 0 0
## 49 86 0 1 0
## 50 90 0 1 0
## 51 70 0 1 0
## 52 76 0 1 0
## 53 65 0 0 1
## 54 69 0 0 1
## 55 60 0 1 0
## 56 70 0 1 0
## 57 95 0 1 0
## 58 80 0 1 0
## 59 54 0 1 0
## 60 90 1 0 0
## 61 86 0 1 0
## 62 165 1 0 0
## 63 175 1 0 0
## 64 150 1 0 0
## 65 153 1 0 0
## 66 150 1 0 0
## 67 208 1 0 0
## 68 155 1 0 0
## 69 160 1 0 0
## 70 190 1 0 0
## 71 97 1 0 0
## 72 150 1 0 0
## 73 130 1 0 0
## 74 140 1 0 0
## 75 150 1 0 0
## 76 112 1 0 0
## 77 76 0 1 0
## 78 87 0 1 0
## 79 69 0 1 0
## 80 86 0 1 0
## 81 92 0 1 0
## 82 97 0 1 0
## 83 80 0 1 0
## 84 88 0 1 0
## 85 175 1 0 0
## 86 150 1 0 0
## 87 145 1 0 0
## 88 137 1 0 0
## 89 150 1 0 0
## 90 198 1 0 0
## 91 150 1 0 0
## 92 158 1 0 0
## 93 150 1 0 0
## 94 215 1 0 0
## 95 225 1 0 0
## 96 175 1 0 0
## 97 105 1 0 0
## 98 100 1 0 0
## 99 100 1 0 0
## 100 88 1 0 0
## 101 95 0 1 0
## 102 46 0 1 0
## 103 150 1 0 0
## 104 167 1 0 0
## 105 170 1 0 0
## 106 180 1 0 0
## 107 100 1 0 0
## 108 88 1 0 0
## 109 72 0 1 0
## 110 94 0 1 0
## 111 90 1 0 0
## 112 85 1 0 0
## 113 107 0 1 0
## 114 90 0 1 0
## 115 145 1 0 0
## 116 230 1 0 0
## 117 49 0 1 0
## 118 75 0 1 0
## 119 91 1 0 0
## 120 112 1 0 0
## 121 150 1 0 0
## 122 110 0 1 0
## 123 122 1 0 0
## 124 180 1 0 0
## 125 95 1 0 0
## 126 100 1 0 0
## 127 100 1 0 0
## 128 67 0 0 1
## 129 80 0 1 0
## 130 65 0 0 1
## 131 75 0 1 0
## 132 100 1 0 0
## 133 110 1 0 0
## 134 105 1 0 0
## 135 140 1 0 0
## 136 150 1 0 0
## 137 150 1 0 0
## 138 140 1 0 0
## 139 150 1 0 0
## 140 83 0 1 0
## 141 67 0 1 0
## 142 78 0 1 0
## 143 52 0 0 1
## 144 61 0 0 1
## 145 75 0 1 0
## 146 75 0 1 0
## 147 75 0 1 0
## 148 97 0 1 0
## 149 93 0 1 0
## 150 67 0 0 1
## 151 95 1 0 0
## 152 105 1 0 0
## 153 72 1 0 0
## 154 72 1 0 0
## 155 170 1 0 0
## 156 145 1 0 0
## 157 150 1 0 0
## 158 148 1 0 0
## 159 110 1 0 0
## 160 105 1 0 0
## 161 110 1 0 0
## 162 95 1 0 0
## 163 110 0 1 0
## 164 110 1 0 0
## 165 129 1 0 0
## 166 75 0 1 0
## 167 83 0 1 0
## 168 100 1 0 0
## 169 78 0 1 0
## 170 96 0 1 0
## 171 71 0 1 0
## 172 97 0 1 0
## 173 97 1 0 0
## 174 70 0 1 0
## 175 90 1 0 0
## 176 95 0 1 0
## 177 88 0 1 0
## 178 98 0 1 0
## 179 115 0 1 0
## 180 53 0 0 1
## 181 86 0 1 0
## 182 81 0 1 0
## 183 92 0 1 0
## 184 79 0 1 0
## 185 83 0 1 0
## 186 140 1 0 0
## 187 150 1 0 0
## 188 120 1 0 0
## 189 152 1 0 0
## 190 100 0 1 0
## 191 105 0 1 0
## 192 81 0 1 0
## 193 90 0 1 0
## 194 52 0 1 0
## 195 60 0 1 0
## 196 70 0 1 0
## 197 53 0 0 1
## 198 100 1 0 0
## 199 78 1 0 0
## 200 110 1 0 0
## 201 95 1 0 0
## 202 71 0 1 0
## 203 70 0 0 1
## 204 75 0 1 0
## 205 72 0 1 0
## 206 102 1 0 0
## 207 150 1 0 0
## 208 88 1 0 0
## 209 108 1 0 0
## 210 120 1 0 0
## 211 180 1 0 0
## 212 145 1 0 0
## 213 130 1 0 0
## 214 150 1 0 0
## 215 68 0 0 1
## 216 80 0 1 0
## 217 58 0 0 1
## 218 96 0 1 0
## 219 70 0 0 1
## 220 145 1 0 0
## 221 110 1 0 0
## 222 145 1 0 0
## 223 130 1 0 0
## 224 110 1 0 0
## 225 105 0 1 0
## 226 100 1 0 0
## 227 98 1 0 0
## 228 180 1 0 0
## 229 170 1 0 0
## 230 190 1 0 0
## 231 149 1 0 0
## 232 78 0 1 0
## 233 88 0 1 0
## 234 75 0 1 0
## 235 89 0 1 0
## 236 63 0 0 1
## 237 83 0 0 1
## 238 67 0 1 0
## 239 78 0 0 1
## 240 97 0 1 0
## 241 110 0 1 0
## 242 110 0 1 0
## 243 48 0 0 1
## 244 66 0 0 1
## 245 52 0 0 1
## 246 70 0 0 1
## 247 60 0 0 1
## 248 110 1 0 0
## 249 140 1 0 0
## 250 139 0 1 0
## 251 105 1 0 0
## 252 95 0 1 0
## 253 85 0 1 0
## 254 88 0 1 0
## 255 100 0 1 0
## 256 90 1 0 0
## 257 105 0 1 0
## 258 85 0 1 0
## 259 110 1 0 0
## 260 120 1 0 0
## 261 145 1 0 0
## 262 165 1 0 0
## 263 139 1 0 0
## 264 140 1 0 0
## 265 68 0 1 0
## 266 95 0 1 0
## 267 97 0 1 0
## 268 75 0 0 1
## 269 95 0 1 0
## 270 105 0 1 0
## 271 85 0 1 0
## 272 97 0 1 0
## 273 103 0 1 0
## 274 125 1 0 0
## 275 115 0 1 0
## 276 133 1 0 0
## 277 71 0 0 1
## 278 68 0 1 0
## 279 115 0 1 0
## 280 85 1 0 0
## 281 88 0 1 0
## 282 90 0 1 0
## 283 110 0 1 0
## 284 130 1 0 0
## 285 129 1 0 0
## 286 138 1 0 0
## 287 135 1 0 0
## 288 155 1 0 0
## 289 142 1 0 0
## 290 125 1 0 0
## 291 150 1 0 0
## 292 71 0 0 1
## 293 65 0 0 1
## 294 80 0 0 1
## 295 80 0 1 0
## 296 77 0 1 0
## 297 125 0 1 0
## 298 71 0 1 0
## 299 90 0 1 0
## 300 70 0 0 1
## 301 70 0 0 1
## 302 65 0 0 1
## 303 69 0 0 1
## 304 90 0 1 0
## 305 115 0 1 0
## 306 115 0 1 0
## 307 90 0 0 1
## 308 76 0 0 1
## 309 60 0 0 1
## 310 70 0 0 1
## 311 65 0 0 1
## 312 90 0 1 0
## 313 88 0 1 0
## 314 90 0 1 0
## 315 90 1 0 0
## 316 78 0 0 1
## 317 90 0 1 0
## 318 75 0 0 1
## 319 92 0 0 1
## 320 75 0 0 1
## 321 65 0 0 1
## 322 105 0 1 0
## 323 65 0 0 1
## 324 48 0 0 1
## 325 48 0 0 1
## 326 67 0 0 1
## 327 67 0 1 0
## 328 67 0 0 1
## 329 67 0 0 1
## 330 62 0 1 0
## 331 132 0 0 1
## 332 100 0 1 0
## 333 88 0 0 1
## 334 72 0 0 1
## 335 84 0 1 0
## 336 84 0 1 0
## 337 92 0 1 0
## 338 110 0 1 0
## 339 84 0 1 0
## 340 58 0 0 1
## 341 64 0 0 1
## 342 60 0 0 1
## 343 67 0 0 1
## 344 65 0 0 1
## 345 62 0 0 1
## 346 68 0 0 1
## 347 63 0 0 1
## 348 65 0 0 1
## 349 65 0 1 0
## 350 74 0 0 1
## 351 75 0 0 1
## 352 75 0 0 1
## 353 100 0 0 1
## 354 74 0 0 1
## 355 80 0 1 0
## 356 76 0 0 1
## 357 116 0 1 0
## 358 120 0 1 0
## 359 110 0 1 0
## 360 105 0 1 0
## 361 88 0 1 0
## 362 85 1 0 0
## 363 88 0 1 0
## 364 88 0 1 0
## 365 88 0 0 1
## 366 85 0 0 1
## 367 84 0 1 0
## 368 90 0 1 0
## 369 92 0 1 0
## 370 74 0 0 1
## 371 68 0 0 1
## 372 68 0 0 1
## 373 63 0 0 1
## 374 70 0 0 1
## 375 88 0 0 1
## 376 75 0 0 1
## 377 70 0 0 1
## 378 67 0 0 1
## 379 67 0 0 1
## 380 67 0 0 1
## 381 110 0 1 0
## 382 85 0 0 1
## 383 92 0 1 0
## 384 112 0 1 0
## 385 96 0 0 1
## 386 84 0 0 1
## 387 90 0 1 0
## 388 86 0 1 0
## 389 52 0 0 1
## 390 84 0 0 1
## 391 79 0 1 0
## 392 82 0 0 1d1 = data x yang lebih kecil dari sama dengan 20 diberi nilai 1, sisanya diberi nilai 0 d2 = data x yang nilainya di antara 20 hingga 30 diberi nilai 1, sisanya diberi nilai 0 d3 = data x yang lebih besar dari 30 diberi nilai 1, sisanya diberi nilai 0
step.mod2 <- lm(AutoData$horsepower~d1+d2+d3)
plot(AutoData$mpg,AutoData$horsepower)
lines(AutoData$mpg,lin.mod2$fitted.values,col="blue")
lines(AutoData$mpg[ix2], step.mod2$fitted.values[ix2],col="green")
Garis hijau merupakan fungsi tangga dengan 3 selang, x<=20 (d1), 20<x<=30 (d2), dan x>30 (d3), sedangkan garis biru merupakan fungsi regresi linier. Apabila garis hijau dan garis biru dibandingkan secara eksploratif, maka dapat dilihat bahwa garis fungsi tangga lebih mengikuti pola dibandingkan garis fungsi regresi linier.
summary(step.mod2)##
## Call:
## lm(formula = AutoData$horsepower ~ d1 + d2 + d3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.450 -12.966 0.034 12.550 92.550
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 70.518 2.886 24.431 < 2e-16 ***
## d11 66.932 3.557 18.816 < 2e-16 ***
## d21 17.448 3.602 4.844 1.84e-06 ***
## d31 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 26.3 on 389 degrees of freedom
## Multiple R-squared: 0.5356, Adjusted R-squared: 0.5332
## F-statistic: 224.3 on 2 and 389 DF, p-value: < 2.2e-16Komparasi Model
nilai_AIC2 <- rbind(AIC(lin.mod2),
AIC(pol.mod2),
AIC(step.mod2))
nama_model <- c("Linear","Poly (ordo=2)","Tangga (breaks=3)")
data.frame(nama_model,nilai_AIC2)## nama_model nilai_AIC2
## 1 Linear 3614.324
## 2 Poly (ordo=2) 3485.217
## 3 Tangga (breaks=3) 3680.688Komparasi model dilakukan berdasarkan nilai AIC-nya. Semakin kecil nilai AIC, maka semakin baik model. Dari ketiga nilai AIC di atas, model polynomial ordo 2 memiliki nilai terkecil.
MSE2 = function(pred,actual){
mean((pred-actual)^2)
}
nilai_MSE2 <- rbind(MSE2(predict(lin.mod2),AutoData$horsepower),
MSE2(predict(pol.mod2),AutoData$horsepower),
MSE2(predict(step.mod2),AutoData$horsepower))
nama_model <- c("Linear","Poly (ordo=2)","Tangga (breaks=3)")
data.frame(nama_model,nilai_MSE2)## nama_model nilai_MSE2
## 1 Linear 582.3257
## 2 Poly (ordo=2) 416.7871
## 3 Tangga (breaks=3) 686.2376Selain menggunakan AIC, komparasi model juga dapat dilakukan berdasarkan nilai MSE-nya.Semakin kecil nilai MSE, semakin bagus modelnya. Sesuai dengan hasil AIC di atas, model polynomial ordo 2 juga memiliki nilai MSE terkecil.
Maka dari itu, baik secara eksploratif maupun berdasarkan nilai AIC dan MSE-nya, model plynomial ordo 2 merupakan model terbaik bagi data tersebut.
plot(AutoData$mpg, AutoData$horsepower,xlim=c( 7,50), ylim=c( 40,250), pch=16,col="orange")
abline(v=25, col="red", lty=2)
Pada plot di atas, kita memberi titik potong (knot) di x=25 . ## Piecewise Cubic Polynomial
dt.all <- cbind(AutoData$horsepower, AutoData$mpg)
##knots=1
dt11 <- dt.all[AutoData$mpg <=25,]
dim(dt11)## [1] 236 2dt22 <- dt.all[AutoData$mpg >25,]
dim(dt22)## [1] 156 2maksud dari knots = 25 yaitu titik potongnya berada pada titik 25. Dimensi data yang berada di samping kiri titik potong (<25) yaitu sebanyak 236 data, sedangkan dimensi data yang berada di kanan titik potong (>25) yaitu sebanyak 156 data.
plot(AutoData$mpg, AutoData$horsepower,xlim=c( 7,50), ylim=c( 40,250), pch=16,col="orange")
cub.mod1 <- lm(dt11[,1]~dt11[,2]+I(dt11[,2]^2)+I(dt11[,2]^3))
ix <- sort(dt11[,2], index.return=T)$ix
lines(dt11[ix,2],cub.mod1$fitted.values[ix],col="blue", lwd=2)
cub.mod2 <- lm(dt22[,1]~dt22[,2]+I(dt22[,2]^2)+I(dt22[,2]^3))
ix <- sort(dt22[,2], index.return=T)$ix
lines(dt22[ix,2],cub.mod2$fitted.values[ix],col="blue", lwd=2)
garis biru di atas merupakan plot Piecewise Cubic Polynomial antara dt1 dan dt2. Terlihat patahan (cutoff) pada titik knots = 25.
Truncated Power Basis
#dengan menggunakan truncated power basis
plot(AutoData$mpg, AutoData$horsepower,xlim=c( 7,50), ylim=c( 40,250), pch=16,col="orange")
abline(v=25, col="red", lty=2)
hx <- ifelse( AutoData$mpg>25,(AutoData$mpg-1)^3,0)
cubspline.mod <- lm( AutoData$horsepower~AutoData$mpg+I(AutoData$mpg^2)+I(AutoData$mpg^3)+hx)
ix <- sort( AutoData$mpg,index.return=T)$ix
lines( AutoData$mpg[ix], cubspline.mod$fitted.values[ix],col="green", lwd=2)
pada plot Truncated Power Basis, plot yang pada awalnya memiliki cutoff
pada titik 25 dihubungkan menjadi satu kesatuan.
Perbandingan dengan k-fold CV
plot(AutoData$mpg, AutoData$horsepower,xlim=c( 7,50), ylim=c( 40,250), pch=16,col="orange")
abline(v=20,col="red", lty=2)
abline(v=40,col="red", lty=2)
hx11 <- ifelse( AutoData$mpg>20,(AutoData$mpg-20)^3,0)
hx22 <- ifelse( AutoData$mpg>40,(AutoData$mpg-40)^3,0)
cubspline.mod22 <- lm( AutoData$horsepower~AutoData$mpg+I(AutoData$mpg^2)+I(AutoData$mpg^3)+hx11+hx22)
ix <- sort( AutoData$mpg,index.return=T)$ix
lines( AutoData$mpg[ix], cubspline.mod$fitted.values[ix],col="green", lwd=2)
Sekarang kita membuat knot pada 2 titik, yaitu pada titik 20 dan titik 40. Maka kita memiliki 3 fungsi, yaitu kurang dari 20, dari 20 sampai 40, dan lebih dari 40.
Perbandingan 1 knot dan 2 knot dengan 5-fold CV
Kita bandingkan nilai residunya
##1 knot (knot=1)
data.all <- cbind( AutoData$horsepower,AutoData$mpg, hx)
set.seed(456)
ind <- sample(1:5,length( AutoData$mpg),replace=T)
res <- c()
for( i in 1:5){
dt.train <- data.all[ ind!=i,]
x.test <- data.all[ ind==i, -1]
y.test <- data.all[ ind==i,1]
mod1 <- lm( dt.train[,1]~dt.train[,2]+I( dt.train[,2]^2)+I( dt.train[,2]^3)+dt.train[,3])
x.test.olah <- cbind(1,x.test[,1], x.test[,1]^2,x.test[,1]^3,x.test[,2])
beta <- coefficients(mod1)
prediksi <- x.test.olah%*%beta
res <- c( res,mean(( y.test-prediksi)^2))
}
res## [1] 337.5968 376.6104 332.1561 465.5146 443.7971mean(res)## [1] 391.135##2 knot (knot = 0, 2)
data.all2 <- cbind(AutoData$horsepower,AutoData$mpg,hx1,hx2)## Warning in cbind(AutoData$horsepower, AutoData$mpg, hx1, hx2): number of rows of
## result is not a multiple of vector length (arg 1)set.seed(456)
ind2 <- sample(1:5,length( AutoData$mpg),replace=T)
res2 <- c()
for( i in 1:5){
dt.train2 <- data.all2[ind2!=i,]
x.test2 <- data.all2[ind2==i,-1]
y.test2 <- data.all2[ind2==i,1]
mod2 <- lm(dt.train2[,1]~dt.train2[,2]+I(dt.train2[,2]^2)+I(dt.train2[,2]^3)+dt.train2[,3]+dt.train2[,4])
x.test.olah2 <- cbind(1,x.test2[,1],x.test2[,1]^2,x.test2[,1]^3,x.test2[,2],x.test2[,3])
beta2 <- coefficients(mod2)
prediksi2 <- x.test.olah2%*%beta2
res2 <- c(res2,mean((y.test2-prediksi2)^2))
}
res2## [1] 362.3182 408.8973 339.1532 479.3794 484.5075mean(res2)## [1] 414.8511Dari perbandingan residu perlipatan, 1 knot lebih bagus daripada 2 knot karena nilai residunya lebih kecil.
Smoothing Spline
Fungsi smooth.spline secara otomatis memilih paramter lambda terbaik dengan menggunakan Cross-Validation (CV) dan Generalized Cross-Validation (GCV). Perbedaan utama antara keduanya adalah pada ukuran kebaikan model yang digunakan. Kalau CV menggunakan MSE sebagai ukuran kebaikan model, sedangkan GCV menggunakan Weighted-MSE. Banyaknya folds (lipatan) pada CV dan GCV ini adalah sejumlah obserbasinya.
ss1 <- smooth.spline(AutoData$mpg,AutoData$horsepower,all.knots = T)
plot(AutoData$mpg, AutoData$horsepower,xlim=c( 7,50), ylim=c( 40,250), pch=16,col="orange")
lines(ss1,col="purple", lwd=2)